小白我最近做的事有些杂乱，是时候整理一波了。
序列的比对及提取已知ID的序列
故事是这样开始的：
故事的主要内容
如图所示，我们遇到了三个问题，当然，这么difficult的问题我是搞不定的啦，在谢大佬的帮助下这三个问题是已经解决了，但我们不满足于此，我们想要把和这些探针序列匹配上的转座子序列提取出来，以便我们进一步设计更多的探针。谢大佬说剩下的工作比较简单，就把他前面的笔记给我们，让我们自己把转座子序列提取出来。于是乎，照猫画虎模式开启。
1st Step
去NCBI把blast+下载到服务器账号，ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/，在这个链接里选择合适的版本，然后
wget -c 网址 —— 解压 —— 添加环境变量 —— 搞定
注意：-c 断点续传 顾名思义，划重点哦
此外还要在citrus.hzau.edu.cn下载甜橙csi.chromosome.fa序列, 并准备探针序列的fasta格式文件。
2nd Step
将探针序列blast到甜橙基因组上，生成格式7，以名为“CL”的探针为例，格式7包含以下信息
格式7包含的信息
blast过程也很简单，可以写个脚本，提交到队列中运行
blast 脚本第一行，打点当前目录下文件，并将其作为运行参数
第二行，建立基因组的索引文件，数据为Nt即DNA
这个 “解析序列ID” 我还没懂是神马意思第三行，将query序列blast到基因组上，输出.out7文件
num_threads参数表明运行blast所用CPU数，应该和向服务器申请的CPU数一致
3rd Step
从fish_csi.out7文件中提取出各个探针Identity>80%, 且coverage>80%的blast结果
提取blast结果，写入*_Ident0.8Cover0.8文件
这是以CL探针为例，我们就是这么一条探针一条探针地提取的，写到脚本文件中运行，也可能有其他更快捷的办法，但因为探针数量不多，为了节约时间成本就没再深入追究。
知识点：对于“&&”的解释可以参看下面的链接，意为 &&左边的command1执行成功(返回0表示成功)后，&&右边的command2才能被执行。与之对应的还有“||”和“()”。
https://www.cnblogs.com/chenggang816/p/10303508.html
4th Step
将上一步生成的每条序列的筛选后的结果文件转换成bed格式，用于下步分析, 通过vim TransToBed.sh新建脚本文件，并在TransToBed.sh脚本文件中写入
TransToBed.sh脚本文件中写入的内容
知识点：①Shell脚本for循环结构如图所示；②“$”为申请参数，形成队列；③“echo”为在Shell中打印，运行脚本后可以看到终端打印出一行行*.Ident0.8Cover0.8的文件，起到监视识别文件是否正确的作用，此外“echo”输出的结尾自带换行符，所以该命令结束后的 [账户名 目录名]$ 是新开一行的，而如果用“cat”命令显示一个结尾无换行符的文本文件后 [账户名 目录名]$ 是紧跟在文档最后一个字符后面的，而不是新开一行，这在有利于在合并FASTA文件之判断合并前的FASTA文件末尾是否有换行符；④bed文件的分隔符为“\t”；⑤awk工具的If,else语句如图；⑥图中提供了将目录下后缀相同的文件全部执行操作后分别输出到加了新后缀名的文件中；⑦在目录下对文件进行批量操作时同一批操作的文件使用相同的后缀名，方便进行批处理（我算是明白为何在NCBI的Gbrowse上检索下载的序列文件有那么长而且整齐的文件名了）
5th Step
将转座子的.gff3文件转为bed格式。谢大佬不知从哪里拿来了citrus中的转座子的注释文件，先看看长什么样子吧
转座子.gff3文件格式
由此看在citrus基因组不同位置上分布的转座子是有不同的ID号的。
知道长啥样，转起格式来岂不是小菜一碟？
.gff3转为.bed文件
6th Step
利用bedtools将每条序列的bed格式文件跟csi的TE的gff注释文件的bed格式进行intersect，并判断每条序列完全被某个TE包含(包含度为80%，填入的数值为已经通过探针长度*0.8计算过的)
bedTools intersect是求两个文件所表示的染色体区域的交集bedtools 更多使用方式可参看以下大佬的链接
https://www.jianshu.com/p/6c3b87301491
提取求交集后得到的转座子的ID
7th Step
我们有了转座子的ID名字下一步就是提取该转座子在染色体上的位置信息了（虽然后来我发现其实用bedtools intersect 的 -wo选项输出的内容中有转座子的位置信息，唉，玩游戏还是要看攻略呀！）康康我是咋弄的吧。
刚开始，感觉自己还很机智
刚开始的错误脚本
结果啥也没输出来，sort排过续后的文件根本没有重复的，后来在uniq后面加了 -f 选项，忽略几行也不奏效，感觉可能是文件太大，于是乎将文件拆分了，把各个探针对应的转座子ID分别提出来
把探针CL对应的转座子ID提出来
这里在第一列的名称行加了“z”，是因为有的探针名字的ASCII码是小于chr的，在sort的时候会被排在chr前面，提取重复的时候就不能提取到位置信息的行
再把所有能和探针有交集的转座子的注释信息提出来，我借鉴了该链接下的回答，得到all_te.bed文件
https://segmentfault.com/q/1010000010766902
这个链接让我学到了新的东西：如何用awk工具对两个文件进行不同却有关联的操作，感谢大佬。结果出来后发现得到的行数比repeat_te.bed的行数少一百多行，我认为应该是有的转座子中包含了两个探针，在awk双文件操作时重复的下标识别的是同一个参数，相当于把ID一样的转座子又合并回同一个转座子了。
之后就是重复第一步了
得到CL探针所对应的位置信息
得到这些我就可以去基因组中提序列了
8th Step
首先，将.loca文件转为了.bed文件，之后用bedtools getfasta工具提取序列
提取序列
然后就可以用BioEdit看啦，虽然结果是 并没有什么卵用 ！ 但还是学到了点东西的。
网友评论