基于全基因组的基因家族分析(3):SlNRAMP家族基因CDS和
今天继续进行下一步,也是序列文件的获取,有了这些数据,我们才可以进行下一步的工作,才能够绘制一些图片。
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
文件
CDS原始序列
处理后的文件,只有单行序列
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序列全部获取得到,接下来可以进行一些图的绘制。
染色体位置信息基因序列