||
作者:钱胜
原载于本实验室微信公众号EvoDevo好好玩(http://mp.weixin.qq.com/s/pGZdmSkwt9Rb9NeCEg-C3Q )
先给大家安利一个网站,生信技能树,里面有几千个帖子,有问题的时候不妨在里面搜索,或许就有你想要解决的问题。技能树里面有一个类似的问题:在里面所搜根据ID提取序列,出来的结果是这样的
这个问题很经典,被当做一个主题,里面有很多用shell,Python,Perl等语言写的脚本来实现相关的功能,原本我也是想偷懒,直接使用论坛的代码先run,在看看这些代码的意思,如果你点进去,你就知道这个主题和我想实现的任务是不一样的。那么,应该如何从给定的gene ID 从Reference genome中提取出相应的序列呢?
这里就需要基因组的注释文件作为桥梁,先给大家看一下(依旧是斑马鱼):
前面5行是注释信息,可以略过,后面5行则告诉我们某个基因,如ENSDARG00000104632 基因的起始终止位置(6732-52059),转录本起始终止位置,外显子的位置和CDS的位置,那么我们的思路就有了,根据gtf file提取出所有基因的起始终止位置,再从中获得我需要的gene list中基因的位置,最后根据bedtools的getfasta提取出目标基因的序列,最后进行blast。OK 思路有了,看代码实现:
awk -F ';''{print $1}'Danio_rerio.GRCz10.87.gtf| awk -F "\t"'{if($3~/gene/) print $1"\t"$4"\t"$5"\t"$9 }'|sed 's/gene_id//g'|sed 's/"//g'>Danio_rerio.GRCz10.geneID.bed
#获得所有基因在染色体上的位置及起始终止位点
下面是我的目标gene ID,这些就是我想要的fasta序列的基因
下面在R语言上进行关联合并,获得目标基因在染色体上的信息
Danio.bed <- read.delim("/home/qians/WGD/zebrafish/Danio_rerio.GRCz10.geneID.txt",sep =" ",header = FALSE)
rownames(Danio.bed)<-Danio.bed$V5
zebrafish_one2one.txt <- read.delim("/home/qians/WGD/zebrafish/zebrafish_one2one.txt", header = FALSE)
rownames(zebrafish_one2one.txt)<- zebrafish_one2one.txt$V1
zebrafish_one2one.bed <-Danio.bed[match(rownames(zebrafish_one2one.txt),rownames(Danio.bed)),c(1:4)]
zebrafish_one2one.bed$V4 <- rownames(zebrafish_one2one.bed)
write.table(zebrafish_one2one.bed,file ="/home/qians/WGD/zebrafish/zebrafish_one2one.bed",sep ="\t", quote = FALSE, col.names = NA)
上面我得到了目标基因在染色体上的位置及起始终止位点,但这不是一个标准的bed文件,我需要把第一列去掉,然后用tab分割,作为bedtools的标准输入
awk -F "\t"'{print $2"\t"$3"\t"$4"\t"$5}' zebrafish_one2one.bed > zebrafish_one2one_final.bed
bedtools getfasta -name -fi Danio_rerio.GRCz10.dna_sm.toplevel.fa -bed zebrafish_one2one_final.bed -fo zebrafish_one2one.fa
可以看到我得到了目标基因的fasta文件,最后运行blast进行比对
makeblastdb -in Danio_rerio.GRCz10.dna_sm.toplevel.fa -dbtype nucl
blastn -db Danio_rerio.GRCz10.dna_sm.toplevel.fa -query zebrafish_one2one.fa -out zebrafish_one2one_blast
可以看到ENSEMBL87版本注释的基因共32266个,blast的基因超过四分之一,所以需要的时间还是挺长的,用screen跑上回宿舍睡觉吧!
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-9-24 03:49
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社