从NCBI基因组数据中获得cds，pep和geneID对应表

在做基因组相关分析时，我们常常需要从基因组中提取cds，并翻译成相应的pep序列。此脚本，以NCBI数据库中标准的基因组序列文件和对应的gff文件为输入文件，快速获得cds序列，pep序列，RNA,Protein和gene的对应关系表等相关文件。

A perl script which deals  with ncbi raw data,and from which get cds ,pep and gene,mRNA and protein ID list.

使用方法如下：

perl    $0    gff_filegenomes_fa_file        file_prefix

#! /usr/bin/perl -w

=head1 #######################################
	
=head1 name
	grep_cds_pep_from_ncbi_genomes_datas.pl

=head1 description
	deal with ncbi raw data,and from which get cds ,pep and gene,mRNA and protein ID list.

=head1 example
	perl  $0 ref_ConCri1.0_top_level.gff3.gz ccr_ref_ConCri1.0_chrUn.fa.gz  mole
	perl  $0 gff_file genomes_fa_file  file_prefix

=head1 author
	original from Xiangfeng Li,xflee0608@163.com
		##2014-4-19/21##

=head1 #######################################
=cut

use strict;
die `pod2text $0` unless @ARGV==3;
my ($gff,$fa,$prefix)=@ARGV;

##deal gff file##
$gff=~/gz$/ ? (open GFF,"gzip -cd $gff|"||die):(open GFF,$gff||die);
my (%mrna,%cds,%pep);
while(<GFF>){
	chomp;
	next if(/^#/);
	my @p=split/\t/,$_;
	my @q=split/;/,$p[8];
	my ($rna,$pep,$nt,$gene);
	my $chr=$p[0];
	if($p[2] eq "mRNA"){
		($rna=$q[0])=~s/ID=//;
		($nt=$q[1])=~s/Name=//;
		($gene=$q[-3])=~s/gene=//;
		$mrna{$chr}{$rna}{strand}=$p[6];
		$cds{$rna}=[$gene,$nt];
	}
        if($p[2] eq "CDS"){
                ($pep=$q[1])=~s/Name=//;
                ($rna=$q[2])=~s/Parent=//;
                push @{$mrna{$chr}{$rna}{nt}},[$p[3],$p[4]];
                $pep{$rna}=$pep;
        }
}
close GFF;

##get id list##
my %anno;
my $id_gene_cds_pep=$prefix."_id_gene_cds_pep.lst";
open ID, ">",$id_gene_cds_pep||die;
foreach my $i(sort keys %pep){
	if($cds{$i}){
		my $out=join "\t",$i,$cds{$i}[0],$cds{$i}[1],$pep{$i};
		print ID $out,"\n";
		($anno{$i}=$out)=~s/\t/\|/g;
	}
}
warn "create file:$id_gene_cds_pep\n";
close ID;

##deal fa file##
my %max;
$fa=~/gz$/ ? (open FA,"gzip -cd $fa|"||die):(open FA,$fa||die);
my $raw_cds=$prefix."_raw_cds";
open CDS1,">$raw_cds"||die;
my ($start,$end);
$/=">";<FA>;$/="\n";
while(<FA>){
        my $name=$1 if(/(\S+)/);
	my $info=(split/\|/,$name)[-1];
        $/=">";
        my $seq=<FA>;
        $/="\n";
        $seq=~s/>|\n+//g;
	my $scaf=$mrna{$info};
	foreach my $k(sort keys %$scaf){
		next if(! exists $scaf->{$k}{nt});
		my @p=@{$scaf->{$k}{nt}};
		my $strand=$$scaf{$k}{strand};
		my $get;
		@p=sort{$a->[0]<=>$b->[0]}@p;
		my $loc1=$p[0][0];
		my $loc2=$p[-1][1];
		my ($get_len,$add,$out,$gene);
		if(exists $anno{$k}){
			$add=$anno{$k}; 
			$gene=(split/\|/,$add)[1];
		}else{next;}
		foreach(@p){
			($start,$end)=@$_[0,1];
			$get.=uc(substr($seq,$start-1,$end-$start+1));
		}
		if($strand eq "+"){
			$get_len=length$get;
			$get=~s/([A-Z]{50})/$1\n/g;
			chop($get) unless($get_len%50);
			$out=">$add LOC=$info :$loc1:$loc2:+ length=$get_len\n$get\n";
			push @{$max{$gene}},[$get_len,$out];
			print CDS1 $out;
		}
		if($strand eq "-"){
			$get=&reverse_complement($get);
			$get_len=length$get;
			$get=~s/([A-Z]{50})/$1\n/g;
			chop($get) unless($get_len%50);
			$out=">$add LOC=$info :$loc1:$loc2:- length=$get_len\n$get\n";
			push @{$max{$gene}},[$get_len,$out];
			print CDS1 $out;
		}
	}
}
warn "create file:$raw_cds\n";
close FA;
close CDS1;

##get max transcript##
my $filter_cds=$prefix."_filter_cds";
open CDS2,">$filter_cds"||die;
my @a;
foreach my $j(keys %max){
	my @trans=@{$max{$j}};
	@trans=sort{$a->[0] cmp $b->[0]}@trans;
	push @a,$trans[-1][1];
}
my @a_new;
foreach(@a){
	my $r=$1 if(/^>rna(\d+)/);
	push @a_new,[$r,$_];
}
my @cds_sort=map{$_->[1]}
		sort{$a->[0] <=> $b->[0]}@a_new;
print CDS2 $_ for@cds_sort;
close CDS2;
warn "create file:$filter_cds\n";

##get raw pep sequences##
my $raw_pep=$prefix."_raw_pep";
open PEP1,">",$raw_pep||die;
my @raw_pep=&cds2pep($raw_cds);
print PEP1 $_ for@raw_pep;
close PEP1;
warn "create file:$raw_pep\n";

##get filter pep sequences##
my $filter_pep=$prefix."_filter_pep";
open PEP2,">$filter_pep"||die;
my @filter_pep=&cds2pep($filter_cds);
print PEP2 $_ for@filter_pep;
close PEP2;
warn "create file:$filter_pep\n";

##add label for cds and pep of filter##
my $label=uc($prefix);
open IN1,$filter_cds||die;
my $filter_cds_label=$prefix."_filter_cds_label";
open OUT1,">",$filter_cds_label||die;
while(<IN1>){
	chomp;
	if(/^>/){
		my @a=split/\|/,$_,2;
		my $name=$a[0]."_$label";
		print OUT1"$name |$a[1]\n";
	}else{print OUT1 "$_\n";}
}
close IN1;
close OUT1;
warn "create file:$filter_cds_label\n";

open IN2,$filter_pep||die;
my $filter_pep_label=$prefix."_filter_pep_label";
open OUT2,">",$filter_pep_label||die;
while(<IN2>){
        chomp;
        if(/^>/){
                my @a=split/\|/,$_,2;
                my $name=$a[0]."_$label";
                print OUT2 "$name |$a[1]\n";
        }else{print OUT2 "$_\n";}
}
close IN2;
close OUT2;
warn "create file:$filter_pep_label\n";

##timing##
my $time=times;
my $time_out=sprintf "%.2f",$time/60;
print "##########\nElapsed Time :$time_out mins\n##########\n";

##subroutine##
sub reverse_complement{
	my ($seq)=shift;
	$seq=reverse$seq;
	$seq=~tr/AaGgCcTt/TtCcGgAa/;
	return $seq;
}

##subroutine##
sub cds2pep{
	my $file=shift;
	my %code = (
                        "standard" =>
                                {
                                'GCA' => 'A', 'GCC' => 'A', 'GCG' => 'A', 'GCT' => 'A',                               # Alanine
                                'TGC' => 'C', 'TGT' => 'C',                                                           # Cysteine
                                'GAC' => 'D', 'GAT' => 'D',                                                           # Aspartic Aci
                                'GAA' => 'E', 'GAG' => 'E',                                                           # Glutamic Aci
                                'TTC' => 'F', 'TTT' => 'F',                                                           # Phenylalanin
                                'GGA' => 'G', 'GGC' => 'G', 'GGG' => 'G', 'GGT' => 'G',                               # Glycine
                                'CAC' => 'H', 'CAT' => 'H',                                                           # Histidine
                                'ATA' => 'I', 'ATC' => 'I', 'ATT' => 'I',                                             # Isoleucine
                                'AAA' => 'K', 'AAG' => 'K',                                                           # Lysine
                                'CTA' => 'L', 'CTC' => 'L', 'CTG' => 'L', 'CTT' => 'L', 'TTA' => 'L', 'TTG' => 'L',   # Leucine
                                'ATG' => 'M',                                                                         # Methionine
                                'AAC' => 'N', 'AAT' => 'N',                                                           # Asparagine
                                'CCA' => 'P', 'CCC' => 'P', 'CCG' => 'P', 'CCT' => 'P',                               # Proline
                                'CAA' => 'Q', 'CAG' => 'Q',                                                           # Glutamine
                                'CGA' => 'R', 'CGC' => 'R', 'CGG' => 'R', 'CGT' => 'R', 'AGA' => 'R', 'AGG' => 'R',   # Arginine
                                'TCA' => 'S', 'TCC' => 'S', 'TCG' => 'S', 'TCT' => 'S', 'AGC' => 'S', 'AGT' => 'S',   # Serine
                                'ACA' => 'T', 'ACC' => 'T', 'ACG' => 'T', 'ACT' => 'T',                               # Threonine
                                'GTA' => 'V', 'GTC' => 'V', 'GTG' => 'V', 'GTT' => 'V',                               # Valine
                                'TGG' => 'W',                                                                         # Tryptophan
                                'TAC' => 'Y', 'TAT' => 'Y',                                                           # Tyrosine
                                'TAA' => 'U', 'TAG' => 'U', 'TGA' => 'U'                                              # Stop
                                }
                        ## more translate table could be added here in future
                        ## more translate table could be added here in future
                        ## more translate table could be added here in future
        );
	open IN,$file||die;
	$/=">";<IN>;$/="\n";
	my @results_set;
	while(<IN>){
		my $info=$_;
        	my @a=split/\s+/,$info;
        	$/=">";my $seq=<IN>;$/="\n";
        	$seq=~s/\n|>//g;
        	my $len=length$seq;
        	my $info_out=join " ",@a[0..($#a-1)];
        	my ($pep_out,$triplet);
        	for(my $i=0;$i<$len;$i+=3){
                	$triplet=substr($seq,$i,3);
                	next if(length$triplet!=3);
                	if(exists $code{standard}{$triplet}){
                        	$pep_out.=$code{standard}{$triplet};
                	}else{$pep_out.="X"}
        	}
        	$pep_out=~s/U$// if($pep_out=~/U$/);
        	my $pep_len=length$pep_out;
        	$pep_out=~s/([A-Z]{50})/$1\n/g;
        	chop($pep_out) unless($pep_len%50);
        	my $results= ">$info_out length=$pep_len\n$pep_out\n";
		push @results_set,$results;
	}
	return @results_set;
}

__END__