基于全基因组的基因家族分析（3）：SlNRAMP家族基因CDS和Genomic DNA序列获取
大家好，我是技能树的老朋友啦，三年前在群主的第一波RNA-seq入门8步活动中因为表现优异获得群主青睐成为技能树VIP一员，也开启了自己的学习经验分享人生！
考虑到技能树过于偏重于肿瘤等疾病领域经验分享，我有必要自告奋勇推荐一下自己的我们植物学领域的生物信息学应用心得体会，会以4个头条的形式发布，也欢迎大家点击原文直达我的博客！
今天继续进行下一步，也是序列文件的获取，有了这些数据，我们才可以进行下一步的工作，才能够绘制一些图片。
1. CDS序列获取
SlNRAMP家族成员鉴定这篇文章已经讲解了如何获得蛋白质序列，那么同样的方法，我们也就可以得到CDS序列，变化的只不过是数据集，将之前的protein.fa改成CDS.fa即可，而且这两个数据我在第一篇文章数据准备中已经下载，不知道的小伙伴可以回看一下。CDS序列获取代码如下：
# Nramp_num文件含有id号
xargs samtools faidx CDS.fa < id > Nramp_cds_seq
less Nramp_cds_seq
# 发现CDS序列有回车，分为很多行，用Perl单行命令将其变成一行
perl -pe '/^>/ ? print "\n" : chomp' Nramp_cds_seq | tail -n +2 > Nramp_cds.txt
less Nramp_cds.txt
2. Genomic DNA序列
Gene序列比起前面的protein和CDS序列要稍微复杂一点，因为我们要获得的ID号不仅仅是geneid，而且还需要在染色体上的位置信息——起始和终止位置，以及染色体编号。
这里就需要用到其他两个数据文件了，就是基因组序列和基因组注释文件，思路——首先根据已经获得的ID号从gff文件中获取染色体位置信息，然后再用bedtools工具根据得到的染色体位置信息来获取基因的序列，最终得到基因集。代码如下：
# 根据geneid批量获取基因的染色体位置信息文件
grep -f geneid.file ITAG2.4_gene_models.gff3 | cut -f 1,4,5,9 | grep "ID=mRNA" | sed 's/;Name.*//g' | sed 's/ID.*://g' > gene.bed
# 利用bedtools工具，将基因序列提取出来并重定向到新文件
bedtools getfasta -fi genome.fa -bed gene.bed -name > gene_seq
代码1解释：首先查看gff3文件的格式，在转录组入门中的了解参考基因组及注释已经做了相应的解释，可以参考那篇文章。然后利用grep -f 将id文件中的geneid号逐行匹配（第一个管道之前的代码）。接着将匹配到的内容用cut -f 进行列截取，包括染色体名，起始位置，终止位置以及基因注释说明（第一个和第二个管道符之间的代码）。因为我们需要的是gene序列也就是编码mRNA的那部分而不是exon或者intron，但是在ID部分都属于该基因。因此还要将切割完的序列匹配ID=mRNA的行，这才是我们真正需要的序列。用grep匹配字符（第二和第三管道符之间的代码）。最后用sed将多余的信息替换成空白，将文本最简化（染色体位置信息图所示）。
代码2解释：使用bedtools工具的getfasta选项，用来提取序列，-fi是用来指定基因组fasta格式文件，-bed用来指定染色体位置信息文件，-name用来确定基因的名字，最后重重定向。可以看到>后的名称是基因ID，染色体及起始和终止位置。
结束了今天的任务，CDS和gene序列全部获取得到，接下来可以进行一些图的绘制。
当然，上面提到的CDS和gene序列你也可以直接在数据库里面下载，比如ENSEMBL数据库，这里就不赘述了！