作者：钱胜

原载于本实验室微信公众号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跑上回宿舍睡觉吧！