perlを使ってgff3ファイルの重複領域を削除するには？

Question

9列のgff3ファイルから重複領域を削除しようとしています。

**Input file:** 
scaffold591 Source gene 3322458 3376057 0.41 - . ID=g24007 
scaffold591 Source transcript 3322458 3376057 0.41 - . ID=g24007.t1;Parent=g24007 
scaffold591 Source transcription_end_site 3322458 3322458 . - . Parent=g24007.t1 
scaffold591 Source gene 3322500 3346055 0.41 - . ID=g24007 
scaffold591 Source transcript 3322500 3346055 0.41 - . ID=g24007.t1;Parent=g24007 
scaffold591 Source transcription_end_site 3322500 3322500 . - . Parent=g24007.t1 
scaffold591 Source gene 3377307 3513095 0.46 + . ID=g24008 
scaffold591 Source transcript 3377307 3513095 0.41 + . ID=g24008.t1;Parent=g24008 
scaffold591 Source transcription_end_site 3377307 3377307 . + . Parent=g24008.t1 

ここでは、同じ鎖の「遺伝子」、すなわち「 - 」または「+」（第7列）を有する行のみを比較しようとしています。 「 - 」鎖（第7列）の例行1と行4

scaffold591 Source gene 3322458 3376057 0.41 - . ID=g24007 
scaffold591 Source gene 3322500 3346055 0.41 - . ID=g24007 

それらが同じ足場と同じから「遺伝子」であるため

。 row4の座標（列4と5）は、行1の座標の範囲内にあります。そのような場合、私のコードは、重複行4を削除し、より大きな範囲を持つrow1を保持する必要があります。二回

My code: 
#!/usr/bin/perl 
use warnings; 
use strict; 
use List::Util qw{ max }; 

open (IN, "<scaffold_sample.txt"); 

my $previous_seqid =""; 
my $previous_strand; 
my $previous_start; 
my $previous_end; 
my @gff; 
my @tmp; 
while (<IN>) 
{ 
    chomp; 
    my ($seqid,$source, $region, $start, $end, $score, $strand, $frame, $attribute) = split ("\t",$_); 
    @gff = ($seqid,$source, $region, $start, $end, $score, $strand, $frame, $attribute); 

    if ($seqid eq $previous_seqid && $strand eq $previous_strand && $region eq 'gene') 
    { 
     if($start < $previous_end && $end < $previous_end) 
     { 
      @gff = @tmp; 
      $previous_seqid = $gff[0]; 
      $previous_strand = $gff[6]; 
      $previous_start = $gff[3]; 
      $previous_end = $gff[4]; 
      print join "\t",@gff; 
      print "\n"; 
     } 
     else 
     { 
      @tmp = @gff; 
     } 

    } 
    else 
    { 
     @tmp = ($seqid,$source, $region, $start, $end, $score, $strand, $frame, $attribute); 
     $previous_seqid = $seqid; 
     $previous_strand = $strand; 
     $previous_start = $start; 
     $previous_end = $end; 
     print join "\t",@tmp; 
     print "\n"; 
    } 

} 

助けてくださいROW1とその次の行

**My expected output:** 
scaffold591 Source gene 3322458 3376057 0.41 - . ID=g24007 
scaffold591 Source transcript 3322458 3376057 0.41 - . ID=g24007.t1;Parent=g24007 
scaffold591 Source transcription_end_site 3322458 3322458 . - . Parent=g24007.t1 
scaffold591 Source gene 3377307 3513095 0.46 + . ID=g24008 
scaffold591 Source transcript 3377307 3513095 0.41 + . ID=g24007.t1;Parent=g24008 
scaffold591 Source transcription_end_site 3377307 3377307 . + . Parent=g24008.t1 

私のコードを印刷します。

Answer 1

これは興味深い問題であることが判明しました。ローをデデップしたいのですが、ファイルの後半で大きな範囲が見つかった場合は、元の小さな範囲が見つかった位置にこの大きな範囲を出力したいと考えています。

正直言って、私はあなたのソリューションを見ていませんでしたが、最初から始めました。

2つのデータ構造を使用しました。 %line_dataには、処理した行の詳細が含まれています。これはマルチレベルのハッシュであり、seqid、strand、およびregionにキーイングされています。新しいレコードがハッシュの値と一致しない場合、我々は最初にseqid、strandおよびregionの組み合わせである。新しいレコードが一致する場合は、前にこの組み合わせを見て、2つのうちどれが最大の範囲を持ち、必要に応じて上書きします。

次に、私たちが出力しようとしているデータを含む@linesがあります。これには、%line_dataのハッシュへの参照が含まれています。より広い範囲が見つかったときに最新の状態に保つためには、ちょっとしたハウスキーピングが必要です。

これが私の結論です。それはあなたの入力に対して正しい出力を与えますが、より多様な入力を壊すかどうかは分かりません。

#!/usr/bin/perl 

use strict; 
use warnings; 
use feature 'say'; 

my @lines; 
my %line_data; 

# Column names (for use as hash keys)  
my @cols = qw[seqid source region start end score strand frame attribute]; 

# Store the input data in DATA for easier testing 
while (<DATA>) { 
    my %record; 
    # Split a record into a hash 
    @record{@cols} = split; 

    # If this key combination exists... 
    if (exists $line_data{$record{seqid}}{$record{strand}}{$record{region}}) { 
    # Get the previous record with these keys... 
    my $prev = $line_data{$record{seqid}}{$record{strand}}{$record{region}}; 
    # See if the new range is larger... 
    if ($record{start} > $prev->{start} and $record{end} > $prev->{end}) { 
     # If so, overwrite it. 
     $line_data{$record{seqid}}{$record{strand}}{$record{region}} = \%record; 
     $lines[$prev->{pos}] = \%record; 
     $record{post} = $prev->{pos}; 
    } 
    } else { 
    # We haven't seen this key combination before. 
    # So just store it. 
    $line_data{$record{seqid}}{$record{strand}}{$record{region}} = \%record; 
    push @lines, \%record; 
    $record{pos} = $#lines; 
    } 
} 

# Having processed the data, we walk the @lines array, 
# de-referencing the hash and joining the values with a space. 
foreach (@lines) { 
    say join ' ', @$_{@cols}; 
} 

__DATA__ 
scaffold591 Source gene 3322458 3376057 0.41 - . ID=g24007 
scaffold591 Source transcript 3322458 3376057 0.41 - . ID=g24007.t1;Parent=g24007 
scaffold591 Source transcription_end_site 3322458 3322458 . - . Parent=g24007.t1 
scaffold591 Source gene 3322500 3346055 0.41 - . ID=g24007 
scaffold591 Source transcript 3322500 3346055 0.41 - . ID=g24007.t1;Parent=g24007 
scaffold591 Source transcription_end_site 3322500 3322500 . - . Parent=g24007.t1 
scaffold591 Source gene 3377307 3513095 0.46 + . ID=g24008 
scaffold591 Source transcript 3377307 3513095 0.41 + . ID=g24008.t1;Parent=g24008 
scaffold591 Source transcription_end_site 3377307 3377307 . + . Parent=g24008.t1