如何使用perl实现ncbi基因组数据格式转换

小编给大家分享一下如何使用perl实现ncbi基因组数据格式转换，希望大家阅读完这篇文章之后都有所收获，下面让我们一起去探讨吧！

ncbi基因组数据格式转换

大家都知道从ncbi下载的基因组数据与正常的基因组数据有所不同，所有的基因转录本CDS的ID都是重新起的，用起来非常不方便，为此写了一perl脚本将其ID转换回来。

用法如下：

Usage
Forced parameter:
-gff        genoma gff file           <infile>   must be given
-fa         genoma fasta file         <infile>   must be given
-o          output gff file           <outfile>   must be given
-out        output fasta file         <outfile>   must be given
Other parameter:
-h           Help document

代码如下：

use Getopt::Long;
use strict;
use Bio::SeqIO;
use Bio::Seq;
#get opts
my %opts;
GetOptions(\%opts, "gff=s", "o=s", "fa=s", "out=s","h");
if(!defined($opts{gff}) || !defined($opts{o}) || !defined($opts{fa}) || !defined($opts{out}) || defined($opts{h}))
{
print <<"Usage End.";
Usage
Forced parameter:
-gff        genoma gff file           <infile>must be given
-fa         genoma fasta file  <infile>must be given
-o          output gff file  <outfile>must be given
-out        output fasta file  <outfile>must be given
Other parameter:
-h           Help document
Usage End.
exit;
}

open(IN,"$opts{gff}") || die "open $opts{gff} failed\n";
open(OUT,">$opts{o}") || die "open $opts{o} failed\n";
my %chr;
my %gene;
my %gene_mrna_num;
my %mrna;
my %mrna_exon_num;
while(<IN>){
if(/^#/){
print OUT $_;
next;
}
chomp;

my @line = split("\t");

###########################################################################
if($line[2] eq "region"){
if($line[8] =~/chromosome/){
$line[8] =~ /chromosome=([^;]*)/;
my $chromosome = $1;
$chr{$line[0]} = $chromosome;
}else{
$chr{$line[0]} = $line[0];
}

}

###########################################################################
if($line[2] eq "gene"){
$line[8] =~ /ID=([^;]*);Name=([^;]*)/;
my $geneid = $1;
my $genename = $2;
$line[8] =~ s/$geneid/$genename/;
$gene{$geneid} = $genename;
$gene_mrna_num{$geneid} = 0;
}

###########################################################################
if($line[2] eq "mRNA"){
$line[8] =~ /ID=([^;]*);Parent=([^;]*)/;
my $mrnaid = $1;
my $geneid = $2;
$line[8] =~ s/$geneid/$gene{$geneid}/;
$gene_mrna_num{$geneid}++;
$line[8] =~ s/$mrnaid/$gene{$geneid}\.$gene_mrna_num{$geneid}/;

$mrna{$mrnaid} = "$gene{$geneid}.$gene_mrna_num{$geneid}";

$mrna_exon_num{$mrnaid} = 0;
}

###########################################################################
if($line[2] eq "exon"){
$line[8] =~ /ID=([^;]*);Parent=([^;]*)/;
my $exonid = $1;
my $mrnaid = $2;
$mrna_exon_num{$mrnaid}++;

$line[8] =~ s/$mrnaid/$mrna{$mrnaid};Name=$mrna{$mrnaid}\.exon$mrna_exon_num{$mrnaid}/;
$line[8] =~ s/ID=$exonid;//;
}

###########################################################################

if($line[2] eq "CDS"){
$line[8] =~ /ID=([^;]*);Parent=([^;]*)/;
my $cdsid = $1;
my $mrnaid = $2;
$line[8] =~ s/$mrnaid/$mrna{$mrnaid}/;
$line[8] =~ s/$cdsid/$mrna{$mrnaid}/;
}

$line[0] = $chr{$line[0]};

print OUT join("\t",@line)."\n";

}
close(IN);
close(OUT);

my $in  = Bio::SeqIO->new(-file => "$opts{fa}" , -format => 'Fasta');
my $out = Bio::SeqIO->new(-file => ">$opts{out}" , -format => 'fasta');
while ( my $seq = $in->next_seq() ) {
#my($id,$sequence)=($seq->id,$seq->seq);
my($id,$sequence,$desc)=($seq->id,$seq->seq,$seq->desc);

my $newSeqobj = Bio::Seq->new(-seq => $sequence,
       -desc => $desc,
       -id => $chr{$id},
 );
 $out->write_seq($newSeqobj);
}

看完了这篇文章，相信你对“如何使用perl实现ncbi基因组数据格式转换”有了一定的了解，如果想了解更多相关知识，欢迎关注亿速云行业资讯频道，感谢各位的阅读！