|
今天继续进行下一步，也是序列文件的获取，有了这些数据，我们才可以进行下一步的工作，才能够绘制一些图片。
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
 文件
 CDS原始序列
 处理后的文件，只有单行序列
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数据库，这里就不赘述了！
|
