今天繼續(xù)進(jìn)行下一步，也是序列文件的獲取，有了這些數(shù)據(jù)，我們才可以進(jìn)行下一步的工作，才能夠繪制一些圖片。
1. CDS序列獲取
SlNRAMP家族成員鑒定這篇文章已經(jīng)講解了如何獲得蛋白質(zhì)序列，那么同樣的方法，我們也就可以得到CDS序列，變化的只不過是數(shù)據(jù)集，將之前的protein.fa改成CDS.fa即可，而且這兩個(gè)數(shù)據(jù)我在第一篇文章數(shù)據(jù)準(zhǔn)備中已經(jīng)下載，不知道的小伙伴可以回看一下。CDS序列獲取代碼如下：
# Nramp_num文件含有id號(hào)
xargs samtools faidx CDS.fa < id > Nramp_cds_seq
less Nramp_cds_seq
# 發(fā)現(xiàn)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序列要稍微復(fù)雜一點(diǎn)，因?yàn)槲覀円@得的ID號(hào)不僅僅是geneid，而且還需要在染色體上的位置信息——起始和終止位置，以及染色體編號(hào)。
這里就需要用到其他兩個(gè)數(shù)據(jù)文件了，就是基因組序列和基因組注釋文件，思路——首先根據(jù)已經(jīng)獲得的ID號(hào)從gff文件中獲取染色體位置信息，然后再用bedtools工具根據(jù)得到的染色體位置信息來獲取基因的序列，最終得到基因集。代碼如下：
# 根據(jù)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文件的格式，在轉(zhuǎn)錄組入門中的了解參考基因組及注釋已經(jīng)做了相應(yīng)的解釋，可以參考那篇文章。然后利用grep -f 將id文件中的geneid號(hào)逐行匹配（第一個(gè)管道之前的代碼）。接著將匹配到的內(nèi)容用cut -f 進(jìn)行列截取，包括染色體名，起始位置，終止位置以及基因注釋說明（第一個(gè)和第二個(gè)管道符之間的代碼）。因?yàn)槲覀冃枰氖莋ene序列也就是編碼mRNA的那部分而不是exon或者intron，但是在ID部分都屬于該基因。因此還要將切割完的序列匹配ID=mRNA的行，這才是我們真正需要的序列。用grep匹配字符（第二和第三管道符之間的代碼）。最后用sed將多余的信息替換成空白，將文本最簡(jiǎn)化（染色體位置信息圖所示）。
代碼2解釋：使用bedtools工具的getfasta選項(xiàng)，用來提取序列，-fi是用來指定基因組fasta格式文件，-bed用來指定染色體位置信息文件，-name用來確定基因的名字，最后重重定向。可以看到>后的名稱是基因ID，染色體及起始和終止位置。
結(jié)束了今天的任務(wù)，CDS和gene序列全部獲取得到，接下來可以進(jìn)行一些圖的繪制。
