很多基因组数据库中只有基因组，CDS和蛋白质序列，而没有基因的序列，这就需要我们自己去提取
1. 输入文件
首先我们需要知道所有基因在染色体上的位置信息，这个可以在基因组注释gff文件中找到，gff文件格式如下：
image.png
我们需要将gene行里的染色体/基因位置/正负链和基因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
# 第一个管道符：匹配第三列为gene的行
# 第二个管道符：打印第1,4,5,7,9列
# 第三个管道符：去掉；以后的内容
# 第四个管道符：去掉ID=gene-
结果如下：
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是输出的基因序列文件
这样我们就将所有基因序列都提取出来了。
