☰

Perlを使用してgff3ファイルから重複する領域を削除する方法は？

tag-icon • tag-icon perl
Perlを使用してgff3ファイルから重複する領域を削除する方法は？

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番目の列）を持つ行のみを比較したいと思います。

たとえば、行1と行4です。

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

これらは、同じ足場と同じ「-」鎖（7列目）からの「遺伝子」です。行4（列4と5）の座標は、行1の座標範囲内にあります。この場合、私のコードは重複する行4を削除し、行1をより広い範囲に保つ必要があります。

**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

私のコードはrow1とそれに続く行を2回印刷します。

**My code:**
#!/usr/bin/perl
use warnings;
use strict;

open (IN, "<scaffold_sample.txt");
#open (OUT, ">output.txt");
my $previous_seqid = "";
my $previous_strand;
my $previous_start;
my $previous_end;
my @gff;
my @tmp;
my @tmp2;
my @transcripts;
while (<IN>)
 { 
    chomp;
    @gff = split ("\t",$_);

   if ($gff[2] eq "gene")
    {
        #print "yes"."\n";
        if($gff[0] eq $previous_seqid && $gff[6] eq $previous_strand)
        {
            if($gff[3] < $previous_end && $gff[4] < $previous_end)
            {
                @tmp2 = @tmp;
                $previous_seqid = $tmp2[0];
                $previous_strand = $tmp2[6];
                $previous_start = $tmp2[3];
                $previous_end = $tmp2[4];
            }
            else
            {
                @gff=@tmp;
                print join "\t",@gff;
                print "\n";
                $previous_seqid = $gff[0];
                $previous_strand = $gff[6];
                $previous_start = $gff[3];
                $previous_end = $gff[4];

            }
        }
        else
        {
            @tmp = @gff;
            $previous_seqid = $tmp[0];
            $previous_strand = $tmp[6];
            $previous_start = $tmp[3];
            $previous_end = $tmp[4];
        }
        print join "\t",@tmp2;
        print "\n"; 
    }
    else
    {
        print join "\t",@gff;
        print "\n";
    }

}
close (IN);

答え1

@Jesvinこれを試してみてください。お役に立てば幸いです。

#!/usr/local/perl

use strict;
use warnings;
use Data::Dumper;

open( IN, "<gene_sorted.txt" );
open(OUT, ">genes_out.txt");
my $prev_scaffold = "";
my $prev_strand   = "";
my @data;

while (<IN>) {
$_ =~ s/\n|\r//g;

my @tmp = split( /\t/, $_ );
if ( $prev_scaffold . "#" . $prev_strand ne $tmp[0] . "#" . $tmp[6] ) {
$prev_scaffold = $tmp[0];
$prev_strand   = $tmp[6];
check( \@data );
undef @data;
push( @data, $_ );
}
else 
{
push( @data, $_ );
}

}
close IN;

check( \@data );

sub check {
my $prev_start = 0;
my $prev_end   = 0;
my @array      = @{ $_[0] };

foreach (@array) {

#print $_."\n";
(
my $scaffold,
my $source,
my $annotation,
my $start,
my $end,
my $score,
my $strand,
my $column,
my $id
) = split( /\t/, $_ );


if ( $start > $prev_start && $end > $prev_end ) {
$prev_start = $start;
$prev_end   = $end;
#print $start. " > " . $end . "\n";
print OUT $_."\n";
}   
} 
}

関連情報