【实验记录】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