如何使用perl合并和处理文件中的多行以生成报告

Question

我是一个新的perl，只是尝试了一些杂乱的代码。

cat input1.txt

  ##gff-version 2
  ##source-version geneious 5.6.4
  Xm_ABL1 Geneious  CDS 1   168 .   +   .   Name=Xm_ABL1;created by=User;modified by=User;ID=w0IVHutPuN4H4FVDCg4sFVRaJjQ.1340919460469.4
  Xm_ABL1 Geneious  CDS 169 334 .   +   .   Name=Xm_ABL1;created by=User;modified by=User;ID=w0IVHutPuN4H4FVDCg4sFVRaJjQ.1340919460469.4
  Xm_ABL1 Geneious  CDS 335 628 .   +   .   Name=Xm_ABL1;created by=User;modified by=User;ID=w0IVHutPuN4H4FVDCg4sFVRaJjQ.1340919460469.4
  Xm_ABL1 Geneious  CDS 629 901 .   +   .   Name=Xm_ABL1;created by=User;modified by=User;ID=w0IVHutPuN4H4FVDCg4sFVRaJjQ.1340919460469.4
  Xm_ABL1 Geneious  CDS 902 985 .   +   .   Name=Xm_ABL1;created by=User;modified by=User;ID=w0IVHutPuN4H4FVDCg4sFVRaJjQ.1340919460469.4
  Xm_ABL1 Geneious  CDS 986 1165    .   +   .   Name=Xm_ABL1;created by=User;modified by=User;ID=w0IVHutPuN4H4FVDCg4sFVRaJjQ.1340919460469.4
  Xm_ABL1 Geneious  CDS 1166    1350    .   +   .   Name=Xm_ABL1;created by=User;modified by=User;ID=w0IVHutPuN4H4FVDCg4sFVRaJjQ.1340919460469.4
  Xm_ABL1 Geneious  CDS 1351    1504    .   +   .   Name=Xm_ABL1;created by=User;modified by=User;ID=w0IVHutPuN4H4FVDCg4sFVRaJjQ.1340919460469.4
  Xm_ABL1 Geneious  BLAST Hit   169 334 .   +   .   
  Xm_ABL1 Geneious  extracted region    1   168 .   +   .   Name=Extracted region from gi|371443098|gb|JH556762.1|;Extracted interval="351297 -> 351464"
  Xm_ABL1 Geneious  extracted region    169 334 .   +   .   Name=Extracted region from gi|371443098|gb|JH556762.1|;Extracted interval="371785 -> 371950"
  Xm_ABL1 Geneious  extracted region    335 628 .   +   .   Name=Extracted region from gi|371443098|gb|JH556762.1|;Extracted interval="372554 -> 372847"
  Xm_ABL1 Geneious  extracted region    629 901 .   +   .   Name=Extracted region from gi|371443098|gb|JH556762.1|;Extracted interval="374760 -> 375032"
  Xm_ABL1 Geneious  extracted region    902 985 .   +   .   Name=Extracted region from gi|371443098|gb|JH556762.1|;Extracted interval="375230 -> 375313"
  Xm_ABL1 Geneious  extracted region    986 1165    .   +   .   Name=Extracted region from gi|371443098|gb|JH556762.1|;Extracted interval="375992 -> 376171"
  Xm_ABL1 Geneious  extracted region    1166    1350    .   +   .   Name=Extracted region from gi|371443098|gb|JH556762.1|;Extracted interval="376575 -> 376759"
  Xm_ABL1 Geneious  extracted region    1351    1504    .   +   .   Name=Extracted region from gi|371443098|gb|JH556762.1|;Extracted interval="376914 -> 377067"

如果输入文件包含( -> )向前箭头，我希望输出如下if($array7=~/.*interval=\"\d+ -> \d+\"$/gm){ $array5="+"；}

cat output1.txt

gi_371443098_gb_JH556762.1  gene    351297  377067  .   +   .   Name=Xm_ABL1;created by=User;modified by=User;ID=w0IVHutPuN4H4FVDCg4sFVRaJjQ.1340919460469.4
gi_371443098_gb_JH556762.1  CDS 351297  351464  .   +   .   Name=Xm_ABL1;created by=User;modified by=User;ID=w0IVHutPuN4H4FVDCg4sFVRaJjQ.1340919460469.4
gi_371443098_gb_JH556762.1  CDS 371785  371950  .   +   .   Name=Xm_ABL1;created by=User;modified by=User;ID=w0IVHutPuN4H4FVDCg4sFVRaJjQ.1340919460469.4
gi_371443098_gb_JH556762.1  CDS 372554  372847  .   +   .   Name=Xm_ABL1;created by=User;modified by=User;ID=w0IVHutPuN4H4FVDCg4sFVRaJjQ.1340919460469.4
gi_371443098_gb_JH556762.1  CDS 374760  375032  .   +   .   Name=Xm_ABL1;created by=User;modified by=User;ID=w0IVHutPuN4H4FVDCg4sFVRaJjQ.1340919460469.4
gi_371443098_gb_JH556762.1  CDS 375230  375313  .   +   .   Name=Xm_ABL1;created by=User;modified by=User;ID=w0IVHutPuN4H4FVDCg4sFVRaJjQ.1340919460469.4
gi_371443098_gb_JH556762.1  CDS 375992  376171  .   +   .   Name=Xm_ABL1;created by=User;modified by=User;ID=w0IVHutPuN4H4FVDCg4sFVRaJjQ.1340919460469.4
gi_371443098_gb_JH556762.1  CDS 376575  376759  .   +   .   Name=Xm_ABL1;created by=User;modified by=User;ID=w0IVHutPuN4H4FVDCg4sFVRaJjQ.1340919460469.4
gi_371443098_gb_JH556762.1  CDS 376914  377067  .   +   .   Name=Xm_ABL1;created by=User;modified by=User;ID=w0IVHutPuN4H4FVDCg4sFVRaJjQ.1340919460469.4
###

如果输入文件包含(<-)反向箭头，则cat output1.txt。如果($array7=~/.*interval=\“\d+ <- \d+\"$/gm){ $array5="-"；}

gi_371443098_gb_JH556762.1  gene    351297  377067  .   -   .   Name=Xm_ABL1;created by=User;modified by=User;ID=w0IVHutPuN4H4FVDCg4sFVRaJjQ.1340919460469.4
gi_371443098_gb_JH556762.1  CDS 351297  351464  .   -   .   Name=Xm_ABL1;created by=User;modified by=User;ID=w0IVHutPuN4H4FVDCg4sFVRaJjQ.1340919460469.4
gi_371443098_gb_JH556762.1  CDS 371785  371950  .   -   .   Name=Xm_ABL1;created by=User;modified by=User;ID=w0IVHutPuN4H4FVDCg4sFVRaJjQ.1340919460469.4
gi_371443098_gb_JH556762.1  CDS 372554  372847  .   -   .   Name=Xm_ABL1;created by=User;modified by=User;ID=w0IVHutPuN4H4FVDCg4sFVRaJjQ.1340919460469.4
gi_371443098_gb_JH556762.1  CDS 374760  375032  .   -   .   Name=Xm_ABL1;created by=User;modified by=User;ID=w0IVHutPuN4H4FVDCg4sFVRaJjQ.1340919460469.4
gi_371443098_gb_JH556762.1  CDS 375230  375313  .   -   .   Name=Xm_ABL1;created by=User;modified by=User;ID=w0IVHutPuN4H4FVDCg4sFVRaJjQ.1340919460469.4
gi_371443098_gb_JH556762.1  CDS 375992  376171  .   -   .   Name=Xm_ABL1;created by=User;modified by=User;ID=w0IVHutPuN4H4FVDCg4sFVRaJjQ.1340919460469.4
gi_371443098_gb_JH556762.1  CDS 376575  376759  .   -   .   Name=Xm_ABL1;created by=User;modified by=User;ID=w0IVHutPuN4H4FVDCg4sFVRaJjQ.1340919460469.4
gi_371443098_gb_JH556762.1  CDS 376914  377067  .   -   .   Name=Xm_ABL1;created by=User;modified by=User;ID=w0IVHutPuN4H4FVDCg4sFVRaJjQ.1340919460469.4
###

因为我是初学者，所以我已经尝试了一些杂乱的代码。

#usr/bin/perl;
use strict;
open(FH,"$ARGV[0]");

while(<FH>){
chomp $_; 
my @array=split("\t");
my $key="$array[2]-$array[0]-$array[1]-$array[2]-$array[3]"; 
if($array[1] eq "CDS"){
$cds_cnt{$key}++;
$cds{$key}="$array[4]\t$array[5]\t$array[6]\t$array[7]";
    }
if($array[1] eq "extracted region"){ 
    (my $pos1,my $pos2)=($array[7]=~/.*interval=\"(\d+) -> (\d+)\"$/gm);
$extract_cnt{$key}++;
$extract{$key}="$pos1\t$pos2";
            }
}

 foreach $i (  sort {$a<=>$b} keys %cds){ 
 my $a=$i; #print "$i\n";
 $a=~s/CDS/extracted region/g;
 if($cds_cnt{$i} == $extract_cnt{$a}){
 #print "$i\t$cds{$i}\n$a\t$extract{$a}\n"; 
 my @array=split /\-/,$i;
 my @pos=split "\t",$extract{$a};
 print "$array[1]\t$array[2]\t$pos[0]\t$pos[1]\t$cds{$i}\n";
   }
  }
 print "###";

更新

我需要在代码中修改的内容

1.要从提取的区域行中获取值(例如array7=/gi|371443098|gb|JH556762.1|/)，它可以是任何值，在其上添加下划线(即gi_371443098_gb_JH556762.1)并在output1.txt中以数组形式打印，如下所示。

2.打印时添加新行作为第一行(gi_371443098_gb_JH556762.1基因)，第三列取CDS的起始值(即351297)，第四列取CDS的结束值(即377067)，第一行打印，如output1.txt所示

3. in /extracted region/ block-all rows for.e.g.Extracted interval="351297 -> 351464“(即前进箭头)将array5打印为"+”符号，包括输出中的基因标题。如果e.g.Extracted interval="351297 <- 351464“(反向箭头)将array5打印为输出中包含基因标题的"-”符号。

Answer 1

看起来您要实现的目标是将标记为CDS的行的详细信息与标记为extracted region的匹配行合并，然后根据一些最小值和最大值，按Name分组，使用前导摘要标题打印合并的结果。对吗？

我假设你所说的$array0和$array2足以将他们结合在一起，但在你的例子中这一点并不是很清楚。

您的第一个问题只是一个正则表达式，我认为您已经掌握了它的一般技巧。我认为问题在于你是如何获取数据的。

要做你要求的第二件事，在你的第一遍中捕获hi和lo值，并存储它们。

我并不打算写一个完整的解决方案，但它就在这里……

use strict;
use warnings;

my $metadata = {}; # hashref to store CDS info in..
my $group = {}; # hashref to store summary/detail in..
my $arrow = { "->" => '+', "<-" => '-' }; # decode arrow to pos/neg

open(FH,"$ARGV[0]");
while(<FH>){
    chomp;
    next if /^#/;
    my @array=split("\t");
    my $key = join(":", $array[0], $array[2]);
    if ($array[1] =~ /CDS/){
        $metadata->{$key} = $array[7];
    }
    if ($array[1] =~ /extracted region/){
        #assert CDS already processed..
        die "No CDS record for $key!\n" unless $metadata->{$key};
        (my $label = $array[7]) =~ s/.*region from (.*)\|;.*/$1/;
        $label =~ s/\|/_/g;
        $group->{$label} ||= { #seed summary if not exists
                pos1 => 1e10,
                pos2 => 0,
                metadata => $metadata->{$key},
                sequences => [],
        };
        (my $pos1, my $arr, my $pos2) = ($array[7]=~/.*interval=\"(\d+) (<?->?) (\d+)\"$/gm);
        # capture hi/lo values for group
        $group->{$label}->{pos1} = $pos1 if $pos1 < $group->{$label}->{pos1};
        $group->{$label}->{pos2} = $pos2 if $pos2 > $group->{$label}->{pos2};
        # push this sequence onto the group's array
        push(@{ $group->{$label}->{sequences} }, [ $pos1, $pos2, $arrow->{$arr} ]);
    }
}
for my $gene (sort keys %{ $group }){
    #write out header
    printf "%s\t%s\t%d\t%d\t.\t%s\t.\t%s\n",
        $gene, 'gene',
        $group->{$gene}->{pos1}, $group->{$gene}->{pos2},
        $group->{$gene}->{sequences}->[0]->[2],
        $group->{$gene}->{metadata};
    foreach my $sequence ( @{ $group->{$gene}->{sequences} } ){
        # write out details
        printf "%s\t%s\t%d\t%d\t.\t%s\t.\t%s\n",
            $gene, 'CDS',
            $sequence->[0], $sequence->[1], $sequence->[2],
            $group->{$gene}->{metadata};
    }
}
print "###\n";

我希望这篇文章已经有了足够的注释，让它有意义。像这样写的，如果你必须在六个月后回来修改它，它将是一个更容易维护的lot。

更新

在第四条注释之后，我修改了这段代码。arrayref块现在捕获$array5的值并将其存储在序列arrayref中。

更新#2

• 请注意，使用$arrow hashref将->和<-解码为所需的符号。(第6行)
• 基因标题基于第一个序列的第3个字段中的值显示+或-。有一个假设，一个基因的所有序列都有相同的方向。(从末尾算起11行)

我认为我们现在已经偏离了问答的快车道，进入了自由软件开发局。我写的不是复杂的代码，它有注释和结构。现在是时候让您掌握它的逻辑了。给我的答案加票。