很多基因組數(shù)據(jù)庫(kù)中只有基因組，CDS和蛋白質(zhì)序列，而沒有基因的序列，這就需要我們自己去提取

1. 輸入文件

首先我們需要知道所有基因在染色體上的位置信息，這個(gè)可以在基因組注釋gff文件中找到，gff文件格式如下：

image.png

我們需要將gene行里的染色體/基因位置/正負(fù)鏈和基因ID提取出來（也就是上圖的第1,4,5，7,9行）。直接使用Linux命令即可：

less GCA_002102615.1_NepCla1.0_genomic.gff.gz| awk '{if($3~/^gene$/)print}' |cut -f 1,4,5,7,9|cut -f 1 -d ";"|sed 's/ID=gene-//'>geneid.txt
# 第一個(gè)管道符：匹配第三列為gene的行
# 第二個(gè)管道符：打印第1,4,5,7,9列
# 第三個(gè)管道符：去掉；以后的內(nèi)容
# 第四個(gè)管道符：去掉ID=gene-

結(jié)果如下：

image.png

2.提取基因序列

有了基因的位置信息就可以到基因組序列中截取基因序列。這里需要用到Perl腳本。腳本如下：

die "perl $0 <genome.fa> <weizhi.txt> <OUT> " unless(@ARGV==3 );
use Math::BigFloat;
use Bio::SeqIO;
use Bio::Seq;
$in = Bio::SeqIO -> new(-file => "$ARGV[0]",
                                  -format => 'Fasta');
$out = Bio::SeqIO -> new(-file => ">$ARGV[2]",
                                  -format => 'Fasta');
my %keep=() ;
open IN,"$ARGV[0]" or die "$!";
my%ref=();
while ( my $seq = $in->next_seq() ) {
     my($id,$sequence,$desc)=($seq->id,$seq->seq,$seq->desc);

         $ref{$id}=$seq;

}

$in->close();
open IN,"$ARGV[1]" or die "$!";
while (<IN>) {
     chomp;
     next if /^#/;
     my @a= split /\t/;
     my$seq=$ref{$a[0]};

     my$seq_string=$seq->subseq($a[1],$a[2]);
     my$newseqobj1=Bio::Seq -> new(-seq => $seq_string,
                     -id => "$a[4]"
     ) ;
     if( $a[3]  eq "-" ){
               my$reseq = $newseqobj1 ->revcom();
               $out->write_seq($reseq);
     }elsif ( $a[3]  eq "+" ){

              $out->write_seq($newseqobj1);
     }

}

$perl  gene.fa.jiequ.pl  genome.fa  geneid.txt  gene.fa
#gene.fa.jiequ.pl為腳本名稱;genome.fa為基因組序列文件;geneid.txt是基因位置信息;gene.fa是輸出的基因序列文件

這樣我們就將所有基因序列都提取出來了。