Cómo obtener toda función en un intervalo desde un archivo GFF3 en Perl?

Question

Me gustaría escribir una función en Perl que recibe un nombre de archivo y una gama GFF3 (es decir, 100.000 .. 2000000). y devuelve una referencia a una matriz que contiene todos los nombres / adhesiones de genes que se encuentran en este intervalo.

supongo usando bioperl tendrá sentido, pero tengo muy poca experiencia con ella. Puedo escribir un script que analiza un GFF3 por mi mismo, pero si se utiliza bioperl (u otra packagae) no es demasiado complicado -. Yo prefiero volver a utilizar su código

Answer 1

use Bio::Tools::GFF;

my $range_start = 100000;
my $range_end   = 200000;

my @features_in_range = ( );


my $gffio = Bio::Tools::GFF->new(-file => $gff_file, -gff_version => 3);

while (my $feature = $gffio->next_feature()) {

    ## What about features that are not contained within the coordinate range but
    ## do overlap it?  Such features won't be caught by this check.            
    if (
        ($feature->start() >= $range_start)
        &&
        ($feature->end()   <= $range_end)
       ) {

        push @features_in_range, $feature;

    }

}

$gffio->close();

AVISO: aplicación Naive. Acabo un golpe de eso, ha tenido ninguna prueba. Ni siquiera voy a garantizar que se compila.

Answer 2

Usted desea utilizar BioPerl para esto, usando posiblemente la Bio::Tools::GFF módulo.

Usted realmente debe preguntar en la BioPerl lista de correo . Es muy amable y los suscriptores están muy bien - que va a duda será capaz de ayudarle. Y una vez que usted consigue una respuesta (y si no se consigue uno aquí en primer lugar), sugiero responder a su propia pregunta aquí con la respuesta de modo que todos puedan beneficiarse!

Answer 3

La siguiente función toma un hash de los objetivos y rangos y devuelve una función que iterar sobre todos los objetivos que se superponen cualquiera de los intervalos. Los objetivos deberían ser una referencia a una serie de referencias:

my $targets =    
[
  [
    $start,
    $end,
  ],
  ...,
]

Los rangos debe ser una referencia a una matriz de hashes:

my $ranges =
[
  {
    seqname   => $seqname,
    source    => $source,
    feature   => $feature,
    start     => $start,
    end       => $end,
    score     => $score,
    strand    => $strand,
    frame     => $frame,
    attribute => $attribute,
  },
  ...,
]

Puede, por supuesto, sólo tiene que pasar un solo objetivo.

my $brs_iterator
= binary_range_search( targets => $targets, ranges => $ranges );

while ( my $gff_line = $brs_iterator->() ) {
   # do stuff
}

sub binary_range_search {
    my %options = @_;

    my $targets = $options{targets}  || croak 'Need a targets parameter';
    my $ranges  = $options{ranges} || croak 'Need a ranges parameter';

    my ( $low, $high ) = ( 0, $#{$ranges} );
    my @iterators = ();

TARGET:
    for my $range (@$targets) {

    RANGE_CHECK:
        while ( $low <= $high ) {

            my $try = int( ( $low + $high ) / 2 );

            $low = $try + 1, next RANGE_CHECK
                if $ranges->[$try]{end} < $range->[0];
            $high = $try - 1, next RANGE_CHECK
                if $ranges->[$try]{start} > $range->[1];

            my ( $down, $up ) = ($try) x 2;
            my %seen = ();

            my $brs_iterator = sub {

                if (    $ranges->[ $up + 1 ]{end} >= $range->[0]
                    and $ranges->[ $up + 1 ]{start} <= $range->[1]
                    and !exists $seen{ $up + 1 } )
                {
                    $seen{ $up + 1 } = undef;
                    return $ranges->[ ++$up ];
                }
                elsif ( $ranges->[ $down - 1 ]{end} >= $range->[0]
                    and $ranges->[ $down - 1 ]{start} <= $range->[1]
                    and !exists $seen{ $down - 1 }
                    and $down > 0 )
                {
                    $seen{ $down - 1 } = undef;
                    return $ranges->[ --$down ];
                }
                elsif ( !exists $seen{$try} ) {
                    $seen{$try} = undef;
                    return $ranges->[$try];
                }
                else {
                    return;
                }
            };
            push @iterators, $brs_iterator;
            next TARGET;
        }
    }

# In scalar context return master iterator that iterates over the list of range iterators.
# In list context returns a list of range iterators.
    return wantarray
        ? @iterators
        : sub {
        while (@iterators) {
            if ( my $range = $iterators[0]->() ) {
                return $range;
            }
            shift @iterators;
        }
        return;
        };
}