从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,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_pepn";
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})/$1n/g;
			chop($get) unless($get_len%50);
			$out=">$add LOC=$info :$loc1:$loc2:+ length=$get_lenn$getn";
			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})/$1n/g;
			chop($get) unless($get_len%50);
			$out=">$add LOC=$info :$loc1:$loc2:- length=$get_lenn$getn";
			push @{$max{$gene}},$out];
			print CDS1 $out;
		}
	}
}
warn "create file:$raw_cdsn";
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_cdsn";

##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_pepn";

##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_pepn";

##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_labeln";

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_labeln";

##timing##
my $time=times;
my $time_out=sprintf "%.2f",$time/60;
print "##########nElapsed Time :$time_out minsn##########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,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})/$1n/g;
        	chop($pep_out) unless($pep_len%50);
        	my $results= ">$info_out length=$pep_lenn$pep_outn";
		push @results_set,$results;
	}
	return @results_set;
}

__END__