Основные регулярные выражения и манипуляции с строкой для анализа ДНК с использованием perl

Question

Я новичок в perl и хотел бы сделать то, что, по моему мнению, является одним из основных манипуляций с строками в последовательности ДНК, хранящиеся в файле rtf.

По сути, мой файл читает (файл в формате FASTA):

>LM1 
AAGTCTGACGGAGCAACGCCGCGTGTATGAAGAAGGTTTTCGGATCGTAA 
AGTACTGTCCGTTAGAGAAGAACAAGGATAAGAGTAACTGCTTGTCCCTT 
GACGGTATCTAACCAGAAAGCCACGGCTAACTACGTGCCAGCAGCCGCGG 
TAATACGTAGGTGGCAAGCGTTGTCCGGATTTATTGGGCGTAAAGCGCGC 
GCAGGCGGTCTTTTAAGTCTGATGTGAAAGCCCCCGGCTTAACCGGGGAG 
GGTCATTGGAAACTGGAAGACTGGAGTGCAGAAGAGGAGAGTGGAATTCC 
ACGTGTAGCGGTGAAATGCGTAGATATGTGGAGGAACACCAGTGGCGAAG 
GCGACTCTCTGGTCTGTAACTGACGCTGAGGCGCGAAAGCGTGGGGAGCA 
AACAGGATTAGATACCCTGGTAGTCCACGCCGT 

То, что я хотел бы сделать, это прочитать в мой файл и напечатать заголовок (заголовок> LM1), то соответствует следующей ДНК последовательности GTGCCAGCAGCCGC, а затем распечатать предыдущую последовательность ДНК.
Так что мой результат будет выглядеть следующим образом:

>LM1 
AAGTCTGACGGAGCAACGCCGCGTGTATGAAGAAGGTTTTCGGATCGTAA 
AGTACTGTCCGTTAGAGAAGAACAAGGATAAGAGTAACTGCTTGTCCCTT 
GACGGTATCTAACCAGAAAGCCACGGCTAACTAC 

Я написал следующую программу:

#!/usr/bin/perl 

use strict; use warnings; 

open(FASTA, "<seq_V3_V6_130227.rtf") or die "The file could not be found.\n"; 

while(<FASTA>) { 
    chomp($_); 
    if ($_ =~ m/^>/) { 
     my $header = $_; 
     print "$header\n"; 
    } 

    my $dna = <FASTA>; 
    if ($dna =~ /(.*?)GTGCCAGCAGCCGC/) { 
     print "$dna"; 
    } 

} 
close(FASTA); 

Проблема заключается в том, что моя программа читает файл построчно, а выход я получаю это следующий:

>LM1 
GACGGTATCTAACCAGAAAGCCACGGCTAACTAC 

в принципе я не знаю, как назначить всю последовательность ДНК моего переменный $ днк и в конечном счете, не знаю, как чтобы избежать прослеживания последовательности ДНК последовательно. Также я получаю это предупреждение: Использование неинициализированного значения $ dna в соответствии с шаблоном (m //) в stacked.pl строке 14, строка 1113.

Если кто-нибудь может мне помочь с написанием лучшего кода или указать мне в правильном направлении это было бы очень признательно.

Answer 1

Используя pos function:

use strict; 
use warnings; 

my $dna = ""; 
my $seq = "GTGCCAGCAGCCGC"; 
while (<DATA>) { 
    if (/^>/) { 
    print; 
    } else { 
    if (/^[AGCT]/) { 
     $dna .= $_; 
    } 
    } 

} 

if ($dna =~ /$seq/g) { 
    print substr($dna, 0, pos($dna) - length($seq)), "\n"; 
} 

__DATA__ 
>LM1 

AAGTCTGACGGAGCAACGCCGCGTGTATGAAGAAGGTTTTCGGATCGTAA 
AGTACTGTCCGTTAGAGAAGAACAAGGATAAGAGTAACTGCTTGTCCCTT 
GACGGTATCTAACCAGAAAGCCACGGCTAACTACGTGCCAGCAGCCGCGG 
TAATACGTAGGTGGCAAGCGTTGTCCGGATTTATTGGGCGTAAAGCGCGC 
GCAGGCGGTCTTTTAAGTCTGATGTGAAAGCCCCCGGCTTAACCGGGGAG 
GGTCATTGGAAACTGGAAGACTGGAGTGCAGAAGAGGAGAGTGGAATTCC 
ACGTGTAGCGGTGAAATGCGTAGATATGTGGAGGAACACCAGTGGCGAAG 
GCGACTCTCTGGTCTGTAACTGACGCTGAGGCGCGAAAGCGTGGGGAGCA 
AACAGGATTAGATACCCTGGTAGTCCACGCCGT 

Вы можете обработать файл с несколькими записями вроде так:

while (<DATA>) { 
    if (/^>/) { 
    if ($dna =~ /$seq/g) { 
     print substr($dna, 0, pos($dna) - length($seq)), "\n"; 
     $dna = ""; 
    } 
    print; 
    } elsif (/^[AGCT]/) { 
    $dna .= $_; 
    } 
} 

if ($dna && $dna =~ /$seq/g) { 
    print substr($dna, 0, pos($dna) - length($seq)), "\n"; 
} 

Answer 2

читать весь файл в память, то обратите внимание на регулярное выражение

while(<FASTA>) { 
    chomp($_); 
    if ($_ =~ m/^>/) { 
     my $header = $_; 
     print "$header\n"; 
    } else { 
    $dna .= $_; 
    } 
} 
if ($dna =~ /(.*?)GTGCCAGCAGCCGC/) { 
    print $1; 
} 

Answer 3

Ваш оператор while читает до конца файла. Это означает, что на каждой итерации цикла, $ _ - следующая строка в <FASTA>. Так что $dna = <FASTA> не делает то, что вы думаете. Он читает больше, чем вы, вероятно, хотите.

while(<FASTA>) { #Reads a line here 
    chomp($_); 
    if ($_ =~ m/^>/) { 
    my $header = $_; 
    print "$header\n"; 
    } 
    $dna = <FASTA> # reads another line here - Causes skips over every other line 
} 

Теперь, вы должны прочитать последовательность в ваш $dna. Вы можете обновить цикл while с помощью инструкции else. Так что, если его линия заголовка, напечатайте его, иначе мы добавим его в $dna.

while(<FASTA>) { 
    chomp($_); 
    if ($_ =~ m/^>/) { 
    # It is a header line, so print it 
    my $header = $_; 
    print "$header\n"; 
    } else { 
    # if it is not a header line, add to your dna sequence. 
    $dna .= $_; 
    } 
} 

После цикла вы можете выполнить регулярное выражение.

Примечание. Это решение предполагает, что в файле fasta имеется только 1 последовательность. Если у вас более одного, ваша переменная $dna будет иметь все последовательности как одну.

Edit: Добавление простой способ обработки нескольких последовательностей

my $dna = ""; 
while(<FASTA>) { 
    chomp($_); 
    if ($_ =~ m/^>/) { 

    # Does $dna match the regex? 
    if ($dna =~ /(.*?)GTGCCAGCAGCCGC/) { 
     print "$1\n"; 
    } 

    # Reset the sequence 
    $dna = ""; 

    # It is a header line, so print it 
    my $header = $_; 
    print "$header\n"; 

    } else { 
    # if it is not a header line, add to your dna sequence. 
    $dna .= $_; 
    } 
} 

# Check the last sequence 
if ($dna =~ /(.*?)GTGCCAGCAGCCGC/) { 
    print "$1\n"; 
} 

Answer 4

Я придумал решение с использованием BioSeqIO (и метод trunc из BioSeq от BioPerl распределения Я также использовал index найти подпоследовательность а. чем использование регулярного выражения.

Это решение не распечатывает id, (строка начинается с>), если подпоследовательность не была найдена или если подпоследовательность начинается с первого положения (и, следовательно, не имеет предшествующих символов).

#!/usr/bin/perl 
use strict; 
use warnings; 
use Bio::SeqIO; 

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

my $out = Bio::SeqIO->new(-file => '>test.dat', 
          -format => 'fasta'); 

my $lookup = 'GTGCCAGCAGCCGC'; 

while (my $seq = $in->next_seq()) { 
    my $pos = index $seq->seq, $lookup; 

    # if $pos != -1, ($lookup not found), 
    # or $pos != 0, (found $lookup at first position, thus 
    # no preceding characters). 
    if ($pos > 0) { 
     my $trunc = $seq->trunc(1,$pos); 
     $out->write_seq($trunc); 
    } 
} 

__END__ 
*** fasta_junk.fasta 
>LM1 
AAGTCTGACGGAGCAACGCCGCGTGTATGAAGAAGGTTTTCGGATCGTAA 
AGTACTGTCCGTTAGAGAAGAACAAGGATAAGAGTAACTGCTTGTCCCTT 
GACGGTATCTAACCAGAAAGCCACGGCTAACTACGTGCCAGCAGCCGCGG 
TAATACGTAGGTGGCAAGCGTTGTCCGGATTTATTGGGCGTAAAGCGCGC 
GCAGGCGGTCTTTTAAGTCTGATGTGAAAGCCCCCGGCTTAACCGGGGAG 
GGTCATTGGAAACTGGAAGACTGGAGTGCAGAAGAGGAGAGTGGAATTCC 
ACGTGTAGCGGTGAAATGCGTAGATATGTGGAGGAACACCAGTGGCGAAG 
GCGACTCTCTGGTCTGTAACTGACGCTGAGGCGCGAAAGCGTGGGGAGCA 
AACAGGATTAGATACCCTGGTAGTCCACGCCGT 

*** contents of test.dat 
>LM1 
AAGTCTGACGGAGCAACGCCGCGTGTATGAAGAAGGTTTTCGGATCGTAAAGTACTGTCC 
GTTAGAGAAGAACAAGGATAAGAGTAACTGCTTGTCCCTTGACGGTATCTAACCAGAAAG 
CCACGGCTAACTAC