Como remover regiões sobrepostas no arquivo gff3 usando perl?

0

Estou tentando remover regiões sobrepostas de um arquivo gff3 de 9 colunas.

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

Aqui estou tentando comparar apenas as linhas com "gene" da mesma cadeia, ou seja, "-" ou "+" (7ª coluna).

Por exemplo, linha 1 e linha 4.

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

Eles são "gene" do mesmo andaime e da mesma linha "-" (7ª coluna). As coordenadas row4 (colunas 4 e 5) estão dentro do intervalo de coordenadas da linha 1. Nesses casos, meu código deve remover a linha 4 sobreposta e reter a linha1 com o intervalo maior.

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

Meu código imprime a linha 1 e suas linhas a seguir duas vezes

**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);
    
por Jesvin Joy 07.11.2017 / 15:34

1 resposta

0

@Jesvin Tente com isso. Espero que seja útil.

#!/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";
}   
} 
}
    
por 07.11.2017 / 17:56

Tags