나는이 문제를 놓치고 있습니다. 디렉토리에 (매우 지저분하고 일관성이)의 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
인쇄물이 얼어 있다고 불평하면 대개 실수로 버퍼링으로 속일 수 있습니다. '$ fh-> autoflush'를 사용하여 버퍼링을 비활성화 할 수 있습니다. – ikegami
정말 인쇄가 멈추는 것은 파이프에 쓰고 버퍼가 가득차 있기 때문입니다. 쓰기가 완료되기 전에 파이프 판독기가 일부 데이터를 읽어야합니다. – ikegami
출력이 파일로 전송되면 변경됩니까? 또는 명령 프롬프트의 설정이이 [답변] (http://superuser.com/a/1053104/544510) –