其他分享
首页 > 其他分享> > 【实验记录】9月9日-9月11日

【实验记录】9月9日-9月11日

作者:互联网

今天的计划就是说,我要把这块自己给自己弄出来点东西来。
希望在汇报前能够完成的事情:
(1) 我要弄明白为什么在基因和表观组层面存在那么多的不一样。
(2) 希望把gwas的结果弄出来。
(3) 找到一些候选的感兴趣的片段,并被各方面的数据认真地证实。

你要认真,你要靠自己,你要专注,你要自己探索出一条可靠的路来~

(一)在基因范围内2kb,以及不在2kb范围内的序列的提取。

Q:现在遇到的问题,就是说两个文件内容之间的比较,在比较之前,需要先对其按照染色体号进行sort,但是这一步经常会出错?
A:
(1)需求:先按照染色体进行排序,接着按照基因组上数字的大小进行排序。
(2)需求:提取两个文件中A文件中有,而B文件中没有的部分。
这个事情,之前有做过类似的。但是总是忘记了代码的备份,所以总是忘记怎么做。而且另一方面,我觉得自己做一些事情,就是“不求甚解”,没有做到理解到位,就很容易一些基本的问题出错。

现在开始继续探索和执行:
http://bbs.chinaunix.net/thread-4298085-1-1.html
这个文件下有若干选择,我们来一个个的试试看。

awk '{match($1,"chr(([0-9]|[A-Z])+)",a);b=a[1]~/[0-9]/?sprintf("%02d",a[1]):a[1];c[b" "$2]=$0}END{PROCINFO["sorted_in"]="@ind_str_asc";for(i in c){print c[i]}}' repeat_2k.bed

这个代码的运行结果是,第一列确实是按照我们所说的chr1-chr22的顺序排列的。但是问题是,第二列反而不是按照数字的大小来排列的,而是按照字典顺序的。所以,还是有问题。以及awk和shell的代码我一直以来都弄不明白。

sort -k1.4,1n -k2,2n repeat_2k.bed >repeat_2k_sort.bed
sort -k1.4,1n -k2,2n repeat_interest.bed >repeat_interest_sort.bed

这个代码倒是实现了我们的需求。
首先,其第一列是按照chrX,chr1-chr22的顺序进行排列的。
其次,其第二列是确乎是按照数字从小到大的顺序排列的。

接下来,我们来研究其究竟是如何实现的。
-k1.4,1n 指的是第一列,从第四个字符开始排序。逗号后面的1n规定结束的位置,指的是只对第一列进行排序。并且按照数值的大小升序排列。
-k2,2n 指的是对第二列,从第二列的第1个字符开始,按顺序只对第二列的数据进行排序,并且按照数值的大小升序排列。

现在的话就是顺利的将两个文件进行了排序,接下来就是取两个文件之间的区别的部分。

diff repeat_2k_sort.bed repeat_interest_sort.bed -y

通过这种side by side的可视化的查看的话,我们发现A文件确乎为B文件的子集。所以,证明我们前面排序确乎是排对了。如果不对的话,容易产生的问题是对错位。
那么,接下来我们想办法是否可以输出两个文件的差级。
参考链接:https://blog.csdn.net/sonia_liss/article/details/102914830

sort repeat_interest_sort.bed repeat_2k_sort.bed repeat_2k_sort.bed | uniq -u >repeat_2k_intergenic.bed
sort -k1.4,1n -k2,2n repeat_2k_intergenic.bed >repeat_2k_intergenic_sort.bed #重新进行排序

这个方法实现了我们想要的东西。

好的,接下来我们要做的事情就是我们想看,在genic和不在genic的序列在Roadmap中的富集结果会有怎样大的差别?

(二)将两个文件与Roadmap索引进行富集分析。

预期的结果:在genic区域之外的序列应该不太可能再富集在Twk的区间内了,或者更有可能在增强子和启动子之类的位置。
我们现在先记录好文件的位置:

/home/xxzhang/workplace/project/geneRegion/transcript_genic_2k_5K_10k

然后找到之前Roadmap的代码,换用新的输入文件重新运行。

bgzip -c repeat_2k_intergenic_sort.bed >repeat_2k_intergenic_sort.bed.gz
bgzip -c repeat_2k_sort.bed >repeat_2k_sort.bed.gz
sed 's/\s\+/\t/g' repeat_2k_intergenic_sort.bed >repeat_2k_intergenic_sort_v2.bed
time giggle search -i split_sort_b -q repeat_2k_sort.bed.gz -s >repeat_2k_sort.bed.giggle.result

real    0m0.659s
user    0m0.620s
sys     0m0.010s

time giggle search -i split_sort_b -q repeat_2k_intergenic_sort.bed.gz -s >repeat_2k_intergenic_sort.bed.giggle.result

real    0m1.107s
user    0m1.049s
sys     0m0.016s

python /home/xxzhang/workplace/software/giggle/scripts/giggle_heat_map.py -s /home/xxzhang/workplace/software/giggle/examples/rme/states.txt -c /home/xxzhang/workplace/software/giggle/examples/rme/EDACC_NAME.txt -i repeat_2k_sort.bed.giggle.result -o repeat_2k_sort.bed.3x11.pdf -n /home/xxzhang/workplace/software/giggle/examples/rme/new_groups.txt --x_size 3 --y_size 11 --stat combo --ablines 15,26,31,43,52,60,72,82,87,89,93,101,103,116,120,122,127 --t "genic 2k" --state_names /home/xxzhang/workplace/software/giggle/examples/rme/short_states.txt --group_names /home/xxzhang/workplace/software/giggle/examples/rme/new_groups_names.txt -159.19499252187038 120.78549515881221  

最后画到了想要的图,这样的结果其实是符合我们的预期的。
也就是说,这些片段,位于基因中的大多数富集在Twk位置,其它的不位于基因区2kb范围内的则大部分是转录不激活的。
接下来我们就是说,通过什么来提现这种关系。或许接下来我们考虑的范围会是,插入到基因区的这部分亚家族中又会有怎样的功能。插入到非基因区的这部分亚家族是否有一些我们比较意外的功能。是我下一步比较感兴趣的问题。
可以联系到咱们现在找到的几个比较显著的家族。
希望Part1和Part2的结果有一个连贯性。

(三)分家族去看序列家族在基因区/非基因区的富集情况

结合part1和通路富集的结果,看看是否有一些有意思的发现。


2022-09-11 我们前面已经知道了,在gene区的片段会更多的位于转录区,而位于intergenic区域的,会更倾向于被异染色质修饰或没有功能。
所以接下来,我们就想分家族看看。
“我”享受双脚踏踏实实的踩在地上的感觉。

(1)genic基因内片段的富集

awk '{print $4}' repeat_2k_sort_v2.bed |sort |uniq >Hs_genic_subfamily.txt

首先是Hs_genic_subfamily.txt文件的构造,这个文件存放的是存在片段在基因间的重复序列的subfamily的名字。然后我们根据这个名字可以提取相应的片段。我们对信息提取主要根据的文件是repeat_2k_sort.bed文件。该文件存放的是位于已知基因间2kb范围内的片段。

其次是构造一些文件夹。

mkdir genic_subfamily
mkdir genic_result
mkdir genic_pdf_result
#!perl
open (MARK, "< Hs_genic_subfamily.txt") or die "can not open it!";
while ($line = <MARK>){
		print($line);
		chomp($line);
		print($line);
		system_call("awk \'\$4==\"".$line."\"\'   /home/xxzhang/data/Epigenome/Roadmap/repeat_2k_sort_v2.bed |sed 's\/\\s\\+\/\\t\/g' |bgzip -c >./genic_subfamily/genic_".$line.".bed.gz");
		system_call("time giggle search -i split_sort_b -q ./genic_subfamily/genic_".$line.".bed.gz -s >./genic_result/genic_".$line.".bed.gz.giggle.result");
		system_call("python /home/xxzhang/workplace/software/giggle/scripts/giggle_heat_map.py  -t \"genic ".$line."\" -s /home/xxzhang/workplace/software/giggle/examples/rme/states.txt -c /home/xxzhang/workplace/software/giggle/examples/rme/EDACC_NAME.txt -i ./genic_result/genic_".$line.".bed.gz.giggle.result -o ./genic_pdf_result/genic_".$line.".bed.gz.3x11.pdf -n /home/xxzhang/workplace/software/giggle/examples/rme/new_groups.txt --x_size 3 --y_size 11 --stat combo --ablines 15,26,31,43,52,60,72,82,87,89,93,101,103,116,120,122,127 --state_names /home/xxzhang/workplace/software/giggle/examples/rme/short_states.txt --group_names /home/xxzhang/workplace/software/giggle/examples/rme/new_groups_names.txt");
} 
close(MARK);

sub system_call
{
  my $command=$_[0];
  print "\n\n".$command."\n\n";
  system($command);
}

这里有一个之前没有做过的尝试,就是如何将许多的pdf文件合并为一个pdf文件,这样老师就不用一个个打开了,只需要拖动即可。
通过查询,发现一些平台有合并的pdf文件的限制,于是就很坑。不推荐。
后来使用了一个cleverPDF的工具,没有合并的文件的限制,使用方便,满足需求。
https://www.cleverpdf.com/cn/combine-pdf

(2)intergenic区域

过程与genic区域同理,略。

在程序运行的过程中,请允许我对结果有一些基本的猜测。我觉得这里应该不太可能在转录区,应该会在enhancer和promoter区域有一定富集的情况。

结果和我们预期的一致,今天的任务完成。

回去看文献和珊珊玩游戏去。

标签:11,sort,repeat,记录,bed,实验,giggle,genic,2k
来源: https://www.cnblogs.com/zjuer/p/16672780.html