2016-07-08 5 views
0

나는이 문제를 놓치고 있습니다. 디렉토리에 (매우 지저분하고 일관성이)의 GenBank 파일 (예 : GBK 파일을 여기에 : ftp://ftp.ncbi.nlm.nih.gov/genomes/archive/old_genbank/Bacteria/Acetobacter_pasteurianus_IFO_3283_01_uid31129/AP011121.gbk)을
1. 프로세스 :
2. 분할 각 파일을 유전자에 의해
3. foreach 문에 나는 펄 스크립트가 루프는 각 유전자에 대한 관련 정보를 얻습니다.
4. 각 루프가 끝날 때 해당 유전자에 대한 정보를 출력합니다.

문제 : 특정 파일의 중간에 무작위로 매달려 있지만 문제에 대해서는 분명히 다릅니다. 그들이 멈추는 유전자는 72K가 넘는 전체 파일에 흩어져 있습니다. 명령 행에 출력 된 결과는 파일에 인쇄 된 출력보다 몇 루프 (유전자)에 있습니다 (그림 참조).

문제 해결 : 다른 파일에 대해 멈추는 변수가 다르다. 변수 인쇄 중 중간에 멈추는 경우가 있는데, RAM/CPU 사용량이 낮다. 메모리/CPU를 계속 사용하면 멈추고 Windows에서는 둘 다 멈 춥니 다. 최신 버전의 딸기 펄)과 리눅스 (플럭스 HPC) 시스템에서는 출력에 아무런 문제가 없다. (커맨드 라인 출력이 파일 출력보다 빠르기 때문에, 그 동안 끊어져있는 유전자를 처리 할 수 ​​있다는 것을 알 수있다. 파일로 인쇄). 그것이 내가 유래에 표시되는 코드의 일부처럼 영리/속기가되지 않도록 코드의 길이 죄송합니다
펄이 매달려 있습니다. 아마도 인쇄와 관련이 있습니다.

, 나는 미생물 학자이다 (도의 GenBank 파일 형식으로 매우 일관성이/언어는 그래서 필요 그 주위의 프로그램). 다른 코드 제안을 구현하게되어 기쁩니다.

#!/usr/bin/perl 
use warnings; 

##### open output files ##### 
$GBKINFO = 'L:\NCBI_DAT\GBFF\Bacteria_tablist.txt'; 
open(GBKINFO,'>', $GBKINFO)||die "unable to open $GBKINFO:$!\n"; 
$debug  = 'L:\NCBI_DAT\GBFF\Bacteria_debug.txt'; 
open(DEBUG,'>', $debug)||die "unable to open $debug:$!\n"; 
$GBKGENOMES = 'L:\NCBI_DAT\GBFF\Bacteria_taxonomy.txt'; 

##### load taxonomy hash ##### 
$taxons = 'R:\1_Downloads\taxa.Bacteria.dat'; 
open(TAXONS, $taxons) || die "unable to open $taxons: $!\n"; 
my %TAXhash; 
while(<TAXONS>){ 
(my $orgID, my $phylog)=split('\t',$_); $TAXhash{$orgID}=$phylog;} 
close(TAXONS); 

##### load protIDs hash ##### 
$protid = 'R:\1_Downloads\geneinfo.Bacteria.dat'; 
open(PROTID, $protid) || die "unable to open $taxons: $!\n"; 
my %IDhash; 
while(<PROTID>){ 
(my $prot, my $ID)=split('\t',$_); 
$prot =~ s/\s//g; $ID =~ s/\n//g; 
$IDhash{$prot}=$ID;} 
close(PROTID); 

##### get genbank files ##### 
my $dir = 'L:\NCBI_DAT\GBFF\Bacteria'; 
die unless opendir DIR, $dir; 
foreach my $file (readdir DIR) { 
    $/="\n//\n"; 
    next if $file eq '.' or $file eq '..'; 
    $gbk = $dir.'/'.$file; 
    $gbk =~ s/\//\\/g; 
    $gbk =~ /GCA_(\d+)/; 

####### open .gbk file and split contigs ####### 
if(-f $gbk && $gbk =~ /\.(gbk|gbff)/ && $gbk !~ /gz$/){ 
    open(GBK, $gbk) || die "unable to open $gbk: $!\n"; 
    $i=1; $count =0; 

    while(<GBK>){ 
     $BIG=$_; 
     $BIG =~ s/[\@<>\%\n]//g; 
     if($BIG =~ /(\s+CDS\s{5}|\s+\S*RNA\s{5})/){ 

      # Get Phylogeny 
      $BIG =~ ~ /db_xref\=\"taxon:(\d+)/; 
      $taxID=$1; $count = 1; 
      $BIG =~ /ORGANISM\s+(.*?)\s+(\w+\;.*?)\./; 
      $SPECIES = $1; $PHYLOGENY = $2; 
      $PHYLOGENY =~ s/\(.*?\)//g; $PHYLOGENY =~ s/[^\w\;]//g; 
      $taxonomy = $TAXhash{$taxID}; 
      if($taxonomy =~ /\w/){ 
       $taxonomy =~ s/(\s+\;$|\s+$|\s+\;\s+$)//g; 
       $PHYLOGENY = $taxonomy;} 
      else{$PHYLOGENY=$PHYLOGENY."\;".$SPECIES;} 
      $PHYLOGENY =~ s/[\t\n]//g; 

      if($PHYLOGENY =~ /Bacteria/i){$org="B";} 
      elsif($PHYLOGENY =~ /Virus/i){$org="V";} 
      elsif($PHYLOGENY =~ /Fungi/i){$org="F";} 
      elsif($PHYLOGENY =~ /Archaea/i){$org="A";} 
      elsif($PHYLOGENY =~ /Chordata/i){$org="C";} 
      elsif($PHYLOGENY =~ /(Viridiplantae|Stramenopiles|Rhodophyta)/i){$org="P";} 
      elsif($PHYLOGENY =~ /Eukaryota/i && $org !~ /[ABCFHPRV]/){$org="I";} 
      else{$org="U";} 

      # Print Phylogeny 
      open(GENO, '>>', $GBKGENOMES)||die "unable to open $GBKGENOMES:$!\n"; 
      print GENO "$taxID\t$PHYLOGENY\n"; close(GENO); 

      # get genome seq ### 
      $BIG =~ /ORIGIN(.+)/; 
      $GenomeSeq=$1; 
      $GenomeSeq =~ s/[^a-z]//ig; 
      $GenomeSeq = uc($GenomeSeq); 
      if(length($GenomeSeq)<100){next;} # eg Bos Taurus genome had no gene seqs 

      # split file by genes 
      $BIG =~ /VERSION\s{5,}(\w.*?)\s/; $Accession = $1; 
      $BIG =~ s/(\s+gene\s{5})/\%$1/g; 
      $BIG =~ s/(\s+[a-z]RNA\s{5})/\%$1/g; 
      $BIG =~ s/(\s+CDS\s{5})/\%$1/g; 
      $BIG =~ s/order\((\d+)\W*.*?\W(\d+)\)+/$1\.\.$2/g; 
      $BIG =~ s/join\((\d+)\W*.*?\W(\d+)\)+/$1\.\.$2/g; 
      @genes = split("\%",$BIG); $junk=shift(@genes); 

      # get gene info 
      foreach(@genes){ 
       $gline = $_; 
       if($gline =~ /^\s+gene\s+/){next;} 

       # get gene type 
       if($gline =~ /\s+CDS\s+[\dc]/)  {$type = "Protein";} 
       elsif($gline =~ /\s\/pseudo/)   {$type = "Pseudo";} 
       elsif($gline =~ /\s+\S*[^mr]RNA\s{5}/){$type = "ncRNA"; 
        $gline =~ /\/note\=\".*\;*(.*)\"/; $LOC=$1; 
        if($LOCUS !~ /\w/){$LOCUS=$LOC;}} 
       elsif($gline =~ /\s+rRNA\s{5}/)  {$type = "rRNA";} 
       elsif($gline =~ /\s+tRNA\s{5}/)  {$type = "tRNA";} 
       else{next;} 

       # get gene names and ids 
       if($gline =~ /\/note\=\".*(COG\d\d\d\d)/)  {$COG = $1; $COG =~ s/\s//g;} else{$COG ='';} 
       if($gline =~ /\/note\=\".*:(K\d\d\d\d\d)/) {$KO = $1; $KO =~ s/\s//g;} else{$KO ='';} 
       if($gline =~ /\/locus_tag\=\"(.*?)\"/)  {$LOCUS = $1; $LOCUS =~ s/\s//g;} else{$LOCUS ='';} 
       if($gline =~ /\/protein_id\=\"(.*?)\"/)  {$ProtID = $1; $ProtID =~ s/\s//g;} else{$ProtID ='';} 
       if($gline =~ /\/product\=\"(.*?)\"/)   {$Product = $1;} else{$Product ='';} 
       if($gline =~ /\/gene\=\"(.*?)\"/)    {$GName = $1;} else{$GName ='';} 
       if($gline =~ /\/inference\=\".*(RF\d+)\"/) {$Rfam = $1;} else{$Rfam ='';} 
       if($gline =~ /\/translation\=\"([\w\s]+)\"/) {$AAseq = $1; $AAseq =~ s/\s//g;} else{$AAseq ='';} 

       # get gene seq and coords     
       if($gline =~/(RNA|CDS)\s+(\d+)\D*\.\.\D*(\d+)/){ 
        if($2>$3){$start = $3; $end = $2;} 
        else{$start = $2; $end = $3;} 
        $strand = "\+"; $seq= substr $GenomeSeq, $start-1, $end-$start+1;} 
       elsif($gline =~/(RNA|CDS)\s+compl\S*?(\d+)\D*\.\.\D*(\d+)/){ 
        if($2>$3){$start = $3; $end = $2;} 
        else{$start = $2; $end = $3;} 
        $strand = "\-"; $seq= substr $GenomeSeq, $start-1, $end-$start+1; 
        $seq =~ tr/atgcrykmbvhdATGCRYKMBVHD/tacgyrmkvbdhTACGYRMKVBDH/; 
        $rseq=reverse($seq); $seq=$rseq;} 
       else{print DEBUG "no coords $gline\t$gbk\n"; next;} 

       $seq=uc($seq); $Glen = length($seq); $coords = "$start\.\.$end"; 
       if($Glen < 5){print DEBUG "gene length issue $gline\t$gbk\n"; next;} 

       # get gene IDs 
       print "prot id $ProtID and $gbk\n"; 
       $IDS = $IDhash{$ProtID}; $IDS =~ s/\n//g; $IDS =~ s/.*\&//; $Func = ''; $DATname = ''; 
       if($IDS =~ /\#/){($DATname, $Func) = split("\#", $IDS); $Func =~ s/(\s+$|^\s+)//;} 
       if($Func !~ /$COG/ && $COG =~ /COG\d\d\d\d/){$Func = $COG."\@".$Func;} 
       if($Func !~ /$KO/ && $KO =~ /K\d\d\d\d\d/){$Func = $KO."\@".$Func;} 
       $Func =~ s/(\@$|^\@)//g; 

       # fix gene name issues 
       if($GName =~ /((hypothetical|uncharacterized|conserved|predicted)\s+protein|unknown function|scaffold|contig)/i || length($GName)<3 || $GName !~ /\w/){ 
         if(length($Product)>length($GName) && $Product !~ /((hypothetical|uncharacterized|conserved|predicted)\s+protein|unknown function|scaffold|contig)/i){$GName=$Product;} 
        elsif(length($DATname)>length($GName) && $DATname !~ /((hypothetical|uncharacterized|conserved|predicted)\s+protein|unknown function|scaffold|contig)/i){$GName=$DATname;} 
       else{ if($type =~ /(protein|pseudo)/i){$GName = "Uncharacterized protein";} else{$GName = "Uncharacterized gene";}}}     
       $GName =~ s/([\;\,\.\@\<\>\%\|]|\(.*\))//g; $GName =~ s/(\s$|^\s)//g; $GName =~ s/\s+/_/g; 
       if($LOCUS !~ /\w/){$LOCUS = $Accession."&".$coords;} 

       $FINAL = "$LOCUS\t$ProtID\t$GName\t$type\t$Glen\t$strand\t$taxID\t$org\t$Func\t$AAseq\t$seq"; $FINAL =~ s/\n//g; 
       if($LOCUS =~ /\w/){print GBKINFO "$FINAL\n";} 

      } # foreach gene 
       last; 
      } # if big matches protein 
      else{$i++; print "no protein $i\n"; next;} 
     } # close while gbk 
     close(GBK)||die "unable to close GBK:$!\n"; #just added to check it is closing 

     if($count==0){ print "no genes unlinked $gbk\n"; unlink $gbk or warn "Could not unlink $gbk: $!"; next;} 
    } # closes 1st if for getting genomes 
} # closes 1st foreach file 

close(GBKINFO); 
close(DEBUG); 

Image showing command line print is ahead of print to file and memory use is fine

+1

인쇄물이 얼어 있다고 불평하면 대개 실수로 버퍼링으로 속일 수 있습니다. '$ fh-> autoflush'를 사용하여 버퍼링을 비활성화 할 수 있습니다. – ikegami

+0

정말 인쇄가 멈추는 것은 파이프에 쓰고 버퍼가 가득차 있기 때문입니다. 쓰기가 완료되기 전에 파이프 판독기가 일부 데이터를 읽어야합니다. – ikegami

+0

출력이 파일로 전송되면 변경됩니까? 또는 명령 프롬프트의 설정이이 [답변] (http://superuser.com/a/1053104/544510) –

답변

0

감사합니다! 나는 일반적인 $ | = 1; 중첩 된 모든 FH와 작동하지 않는 FH 특정 다음을 사용합니다.
select ((select (FH), $ | = 1) [0]);
그것은 매달려있는 곳을 찾는데 도움이되었습니다 ... 하나의 정규식은 일부 gbk 파일의 엉망으로 잘 맞지 않습니다. 나쁜 정규 표현식 -> $ gline = ~ // note \ = \ "\; (. *) \"/;

관련 문제