資源描述：
《grep_cds_pep_from_ncbi_genomes_datas.pl 英文文獻(xiàn)資料》由會(huì)員上傳分享，免費(fèi)在線閱讀，更多相關(guān)內(nèi)容在學(xué)術(shù)論文-天天文庫。
1、#!/usr/bin/perl-w=head1#######################################=head1namegrep_cds_pep_from_ncbi_genomes_datas.pl=head1descriptiondealwithncbirawdata,andfromwhichgetcds,pepandgene,mRNAandproteinIDlist.=head1exampleperlref_ConCri1.0_top_level.gff3.gzccr_ref_ConCri1.0_chrUn.fa.gzmoleperlgff_filege
2、momes_fa_filefile_prefix=head1authororiginalfromXiangfengLi,xflee@163.com##2014-4-19/21##=head1#######################################=cutusestrict;die`pod2text$0`unless@ARGV==3;my($gff,$fa,$prefix)=@ARGV;##dealgfffile##$gff=~/gz$/?(openGFF,"gzip-cd$gff
3、"
4、
5、die):(openGFF,$gff
6、
7、die);my(%mrna,%cd
8、s,%pep);while(){chomp;nextif(/^#/);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/N
9、ame=//;($rna=$q[2])=~s/Parent=//;push@{$mrna{$chr}{$rna}{nt}},[$p[3],$p[4]];$pep{$rna}=$pep;}}closeGFF;##getidlist##my%anno;openID,">",$prefix."_id_gene_cds_pep.lst"
10、
11、die;foreachmy$i(sortkeys%pep){if($cds{$i}){my$out=join"t",$i,$cds{$i}[0],$cds{$i}[1],$pep{$i};printID$out,"";($anno{$i}=$out
12、)=~s/t/
13、/g;}}closeID;##dealfafile##my%max;$fa=~/gz$/?(openFA,"gzip-cd$fa
14、"
15、
16、die):(openFA,$fa
17、
18、die);my$raw_cds=$prefix."_raw_cds";openCDS1,">$raw_cds"
19、
20、die;my($start,$end);$/=">";;$/="";while(){my$name=$1if(/(S+)/);my$info=(split/
21、/,$name)[-1];$/=">";my$seq=;$/="";$seq=~s/>
22、
23、+//g;my$scaf=$mrna{$info};foreachmy$k(sortkeys%$scaf){nextif(!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/
24、/,$a
25、dd)[1];}else{next;}foreach(@p){($start,$end)=@$_[0,1];$get.=uc(substr($seq,$start-1,$end-$start+1));}if($strandeq"+"){$get_len=length$get;$get=~s/([A-Z]{50})/$1/g;chop($get)unless($get_len%50);$out=">$addLOC=$info:$loc1:$loc2:+length=
