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

Question

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

猫输入1.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($array[7]=~/.*interval=\"\d+ -> \d+\"$/gm){ $array[5]= 这样的输出"+"; }

猫输出1.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 如果输入文件包含 (<-) 反向箭头。if($array[7]=~/.*interval=\"\d+ <- \d+\"$/gm){ $array[5]="-"; }

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.要从提取区域的行中获取值，（即array[7]=/gi|371443098|gb|JH556762.1|/）它可以是任何值，添加下划线（即gi_371443098_gb_JH556762.1）并打印如图所示，在 output1.txt 中的数组 [0] 中。

2.在打印时添加新行作为第一行（gi_371443098_gb_JH556762.1基因），在第3列获取CDS的起始值（即351297）并在第4列（即377067）获取CDS的结束值，并在第一行打印为显示在 ouput1.txt

3.如果/extracted region/ block-all rows for.egExtracted interval="351297 -> 351464"（即向前箭头）将array[5]打印为“+”符号，包括输出中的基因标题。如果 egExtracted interval="351297 <- 351464"（反向箭头）将数组 [5] 打印为“-”符号，包括输出中的基因标题。

score 2 · Accepted Answer

看起来您想要实现的是将标记为CDS的行中的详细信息与标记为提取区域的匹配行合并，然后根据一些最小值和最大值，使用前导摘要标题打印合并的结果，分组为名称。那是对的吗？

我将假设您所说的 $array[0] (Xm_ABL1 Geneious) 和 $array[2] (169, 335 等) 足以将它们结合在一起，但这在您的示例中并不是很清楚。

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

要执行您要求的第二件事，请在第一遍中捕获 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";

我希望这有足够的评论才能说得通。这么写，如果过了六个月还得回来修改，维护起来会容易很多。

更新

我在第四条评论之后修改了这段代码。$array[7] 正则表达式块现在捕获 $array[5] 的值并将其存储在序列 arrayref 中。

更新#2

注意使用 $arrow hashref 来解码你需要的符号->。<-（第 6 行）
基因标头根据第一个序列的第三个字段中的值显示 + 或 -。假设一个基因的所有序列都具有相同的方向。（从末尾开始 11 行）

我想我们现在已经偏离了问答高速公路，进入了自由软件开发局。我写的不是复杂的代码，它有注释和结构。是时候让你掌握它的逻辑了。并投票赞成我的答案。

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

1 回答 1

Related

Reference