извлечение последовательностей fasta на основе позиции

В одну сторону. Содержание script.pl:

#!/usr/bin/env perl

use warnings;
use strict;

my ($adn, $l, $header);

while ( <> ) { 
    chomp;

    ## First line is known, a header, so print it and process next one.
    if ( $. == 1 ) { 
        printf qq|%s_%s\n|, $_, q|200-300|;
        next;
    }   

    ## Concat adn while not found a header.
    if ( '>' ne substr $_, 0, 1 ) { 
        if ( ! $l ) { $l = length }
        $adn .= $_; 
        if ( ! eof ) { next }
    }   
    else {
        $header = sprintf qq|%s_%s\n|, $_, q|200-300|;
    }   

    ## Extract range 200-300 and insert newlines to set same length of 
    ## line than before.
    my $s = substr $adn, 199, 100;
    $s =~ s/(.{$l})/$1\n/g;
    printf qq|%s\n|, $s; 
    undef $adn;

    ## If not end of file, print the header of the following adn.
    if ( ! eof ) { print $header }
}

Запустите это так:

perl script.pl infile

Это дает:

>Contig[0001]_200-300
AGAAATCGACGTTTTAAGTTTCTGTTTCCAAATTCAAACGGATGTATCTTCGCCAATAAT
TGTCAGAAGTTAGAATTTCTTTCAACATTATGAAGCCCTT

Я предполагаю, что вы знаете, как включить это в свою программу. Взгляните на функцию substr. Он делает то, что вы хотите.

substr EXPR,OFFSET,LENGTH 

Итак, что вам нужно сделать, это:

my $snippet = substr($sequence, 200, 100);

Обрати внимание на CPAN: есть модуль под названием Bio :: SeqReader :: Fasta, который можно использовать для чтения файла и получения последовательностей. На мой взгляд, это довольно хорошо задокументировано, и я заинтригован этим. Итак, вот решение.

use strict;
use warnings;
use feature qw(say);
use Bio::SeqReader::Fasta;

# Put your positions here (start, end)
my @positions = (
  [ 200, 300 ], 
  [ 50, 180 ],
);

open my $fh, '<', '/path/to/file.fasta' or die $!;
my $in = new Bio::SeqReader::Fasta( fh => $fh ); # create the B::SR::F object

my $i = 0;
while ( my $so = $in->next() ) {
  # $so represents one dataset and is a Bio::SeqReader::FastaRecord

  say substr($so->seq(),               # we want a part of the sequence string
    $positions[$i]->[0],               # starting position
    $positions[$i]->[1] - $positions[$i]->[0]); # end - start --> length
}

close $fh;

Один из способов использования Perl и модуля Bio :: SeqIO. Беги как:

./process_fasta.pl file.fa 200 300

Содержание ./process_fasta.pl:

#!/usr/bin/perl 

use strict;
use warnings;
use Bio::SeqIO;

my $in_file = $ARGV[0];
my $start_pos = $ARGV[1];
my $end_pos = $ARGV[2];

my $in = Bio::SeqIO->new ( -file => $in_file, -format => 'fasta');
my $out = Bio::SeqIO->new( -file => ">$in_file.out", -format => 'fasta');


while (my $seq = $in->next_seq() ) {

    $seq->display_id( $seq->display_id() . "_$start_pos-$end_pos" );
    $out->write_seq( $seq->trunc($start_pos, $end_pos) );
}

Полученные результаты:

>Contig[0001]_200-300
AGAAATCGACGTTTTAAGTTTCTGTTTCCAAATTCAAACGGATGTATCTTCGCCAATAAT
TGTCAGAAGTTAGAATTTCTTTCAACATTATGAAGCCCTTT

Я аплодирую тем, кто использовал Bio:: модули, потому что я предпочитаю их писать что-то новое. Тем не менее, возможно, будет полезно следующее:

use strict;
use warnings;

my $end   = pop;
my $start = pop;
local $/ = '>';

while (<>) {
    chomp;
    next unless /(Contig.+)/;
    my ($header) = "$/$1_$start-$end\n";
    my $seq = ${^POSTMATCH};
    $seq =~ s/\s//g;
    print $header;
    print +( substr $seq, $start - 1, $end - $start ) . "\n";
}

Использование: perl script.pl inFile start end [>outFile]

Последний необязательный параметр направляет вывод в файл.

Пример: perl script.pl data.fasta 200 300

Вывод в вашем наборе данных:

>Contig[0001]_200-300
AGAAATCGACGTTTTAAGTTTCTGTTTCCAAATTCAAACGGATGTATCTTCGCCAATAATTGTCAGAAGTTAGAATTTCTTTCAACATTATGAAGCCCTT

Параметры начала и конца pop отключены @ARGV, а затем разделитель записей устанавливается на ">". По мере чтения файла - записи fasta за раз - информация заголовка захватывается, оставляя последовательность в ${^POSTMATCH}. Все пробелы удаляются из последовательности. Наконец, переформатированный заголовок printed, как и диапазон символов в последовательности.

Надеюсь это поможет!

Полный рабочий скрипт наилегчайшего веса:

#!/usr/bin/env perl

use strict;
use warnings;

my $first=1;
if ( $ARGV[0] eq '-0' )
{
    shift @ARGV;
    $first=0;
}

die("Usage: cat <coords> | get_subseqs.pl (-0) <fasta files> > out; where coords is id, start, end. Use -0 if coords start with 0 instead of 1.\n") unless @ARGV;

# read coords
my %contigs = (); # id => [ start, end ]
while (<STDIN>) {
    chomp;
    my @row = split(/\t/);
    die("Require tab-delim: id, start, end\n") unless @row == 3;
    $contigs{$row[0]} = [ $row[1], $row[2] ];
}

# get subseqs
my ($id,$seq,$start,$end);
foreach my $infile (@ARGV) {
    open(IN, '<', $infile) or die($!);
    while (<IN>) {
        if (/>(\S+)/) {
            my $id2 = $1;
            print ">${id}_$start-$end\n", reformat(substr($seq, $start-$first, $end-$start+1)) if $id;
            if ( exists($contigs{$id2}) ) {
                ($id,$seq,$start,$end) = ($id2,'',@{$contigs{$id2}});
                delete($contigs{$id2});
            } else { $id=undef }
        } elsif($id) { $seq .= $_ }
    }
    close(IN);
    print ">${id}_$start-$end\n", reformat(substr($seq, $start-$first, $end-$start+1)) if $id;
    $id = undef;
}

sub reformat { # add newline every 60 bases
    my $seq = shift;
    my $seq2 = '';
    while ( length($seq) > 60 ) {
        $seq2 .= substr($seq,0,60)."\n";
        $seq = substr($seq,60);
    }
    $seq2 .= $seq."\n";
    return $seq2;
}