序列的比对及提取已知ID的序列

小白我最近做的事有些杂乱,是时候整理一波了。

序列的比对及提取已知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看啦,虽然结果是 并没有什么卵用 ! 但还是学到了点东西的。

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 194,457评论 5 459
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 81,837评论 2 371
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 141,696评论 0 319
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 52,183评论 1 263
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 61,057评论 4 355
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 46,105评论 1 272
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 36,520评论 3 381
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 35,211评论 0 253
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 39,482评论 1 290
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 34,574评论 2 309
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 36,353评论 1 326
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 32,213评论 3 312
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 37,576评论 3 298
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 28,897评论 0 17
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 30,174评论 1 250
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 41,489评论 2 341
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 40,683评论 2 335

推荐阅读更多精彩内容