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

Question

Я новичок в perl. Все еще учусь.

У меня есть файл в формате fasta. Я хотел бы извлечь последовательности, охватывающие определенную позицию. Так, например, от позиции 200 до 300

>Contig[0001] 
TGCATCAAAAGCTGAAAATATGTAGTCGAGAAGTCATTTCGAGAAATTGACGTTTTAAGT 
TTCGGTTTCCAAATTCAACCGGATGTATCTTCGCCAATAATTGTCAGCAGTTAGAATTTC 
TTTCAACATTATGAAGCCCTTTTTATATATTTTGATTCTGCATCAAAAGCTGAAAATATG 
TAGTCTTGAAGTCATTTCGAGAAATCGACGTTTTAAGTTTCTGTTTCCAAATTCAAACGG 
ATGTATCTTCGCCAATAATTGTCAGAAGTTAGAATTTCTTTCAACATTATGAAGCCCTTT 
TTATATATTTTGATTCTGCATCAAAAGCTGAAAATGTGTAGTCTCGAAGTCATTTCGAGA 
AATTGACGTTTTAAGTTTCTGTTTCCAAATTCAAACGGATGTATCTTCGCCAATAATTGT 
CAGAAGTTAGAATTTCTTTCAACATTATGAAGCCCTTTTTACATATTTTGACCCTGCATC 
AAAAGCTGAAAATATGTAGTCTCGAAGTCATTTTGAGAAGTTAGAATTTCTTTCAACATT 
ATGAAGCCCTTTTTATATATTTTGATTCTGCATCAAAAGCTGAAAATATGTAGTCTCGAA 
GTCWTTTCRAGAAATTGACGTTTTAAGTTTCTGTTTCCAAATTCAAACGGATGTATCTTC 
GCCAATAATTGTCAGAAGTTAGAATTTCTTTCAACATTATGAAGCCCTTTTTATATATTT 
TGACTCTGCATCAAAAGCTGAAAATATGTAGTCTCGAAGTCATTTCGAGAAATTGACGTT 

Я хотел бы извлечь последовательности из позиции 200-300 из последовательности Contig[0001]. Выходной сигнал будет:

>Contig[0001]_200-300 
AGAAATCGACGTTTTAAGTTTCTGTTTCCAAATTCAAACGGATGTATCTTCGCCAATAATT 
GTCAGAAGTTAGAATTTCTTTCAACATTATGAAGCCCTTT 

У меня есть почти 500 последовательностей в моем файле Fasta и у меня есть необходимые postions в разделителями табуляцией файл, содержащий идентификатор начать конец.

Будет здорово, если кто-нибудь сможет мне помочь в этом.

Большое вам спасибо за вашу помощь. Я не уверен, что могу предоставить файл, содержащий информацию о позициях.

новичок

Answer 1

Я предполагаю, что вы знаете, как получить это в вашу программу. Посмотрите на функцию 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; 

Answer 2

В одну сторону. Содержание 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 

Answer 3

Один из способов использования 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 

Answer 4

Я приветствую тех, кто использовал 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 PED от @ARGV, а затем разделитель записи установлен в положение ">". По мере чтения файла - записи fasta за раз - информация заголовка захватывается, оставив последовательность в ${^POSTMATCH}. Все пробелы удаляются из последовательности. Наконец, переформатированный заголовок равен print ed, равно как и диапазон символов в последовательности.

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

Answer 5

Полный сценарий рабочего наилегчайшем:

#!/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; 
}