☰

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과 그 후속 행을 두 번 인쇄합니다.

**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";
}   
} 
}

관련 정보