podstawowe manipulowanie ciągiem i ciągiem dla analizy DNA przy użyciu perl

Question

Jestem nowym użytkownikiem języka Perl i chciałbym zrobić coś, co moim zdaniem jest podstawową manipulacją ciągami do sekwencji DNA przechowywanych w pliku rtf.

Zasadniczo mój plik czyta (plik jest w formacie FASTA):

>LM1 
AAGTCTGACGGAGCAACGCCGCGTGTATGAAGAAGGTTTTCGGATCGTAA 
AGTACTGTCCGTTAGAGAAGAACAAGGATAAGAGTAACTGCTTGTCCCTT 
GACGGTATCTAACCAGAAAGCCACGGCTAACTACGTGCCAGCAGCCGCGG 
TAATACGTAGGTGGCAAGCGTTGTCCGGATTTATTGGGCGTAAAGCGCGC 
GCAGGCGGTCTTTTAAGTCTGATGTGAAAGCCCCCGGCTTAACCGGGGAG 
GGTCATTGGAAACTGGAAGACTGGAGTGCAGAAGAGGAGAGTGGAATTCC 
ACGTGTAGCGGTGAAATGCGTAGATATGTGGAGGAACACCAGTGGCGAAG 
GCGACTCTCTGGTCTGTAACTGACGCTGAGGCGCGAAAGCGTGGGGAGCA 
AACAGGATTAGATACCCTGGTAGTCCACGCCGT 

Co chciałbym zrobić, to przeczytać w moim pliku i drukować nagłówek (header jest> LM1), a następnie dopasować następujący DNA sekwencja GTGCCAGCAGCCGC, a następnie wydrukuj poprzednią sekwencję DNA.
Więc moje wyjście będzie wyglądać następująco:

>LM1 
AAGTCTGACGGAGCAACGCCGCGTGTATGAAGAAGGTTTTCGGATCGTAA 
AGTACTGTCCGTTAGAGAAGAACAAGGATAAGAGTAACTGCTTGTCCCTT 
GACGGTATCTAACCAGAAAGCCACGGCTAACTAC 

Napisałem następujący program:

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

Problem polega na tym, że mój program odczytuje pliku linia po linii, a wyjście otrzymuję jest co następuje:

>LM1 
GACGGTATCTAACCAGAAAGCCACGGCTAACTAC 

w zasadzie nie wiem jak przypisać całą sekwencję DNA do mojego zmiennej $ DNA i ostatecznie nie wiem jak aby uniknąć czytania sekwencji sekwencji DNA po linii. Otrzymuję również ostrzeżenie: Użycie niezainicjowanej wartości $ dna w dopasowaniu wzorca (m //) w linii stacked.pl 14, wiersz 1113.

Jeśli ktoś mógłby mi pomóc w napisaniu lepszego kodu lub wskazaniu mi w odpowiednim kierunku byłoby to bardzo cenne.

Answer 1

Używanie 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 

można przetworzyć plik z wieloma wpisami tak:

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

przeczytać cały plik do pamięci, a następnie szukać regexp

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

Answer 3

Instrukcja Your while czyta do końca pliku. Oznacza to, że przy każdej iteracji pętli, $ _ jest następną linią w <FASTA>. Tak więc $dna = <FASTA> nie robi tego, co myślisz, że jest. Czyta więcej niż prawdopodobnie tego chcesz.

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 
} 

Teraz należy przeczytać sekwencję w swoim numerze $dna. Możesz zaktualizować pętlę while za pomocą instrukcji else. Więc jeśli jest to linia nagłówkowa, wydrukuj ją, w przeciwnym razie dodamy ją do $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 .= $_; 
    } 
} 

Po pętli możesz zrobić swoje wyrażenie regularne.

Uwaga: to rozwiązanie zakłada, że ​​w pliku fasta znajduje się tylko 1 sekwencja. Jeśli masz więcej niż jeden, twoja zmienna $dna będzie miała wszystkie sekwencje jako jedna.

Edycja: Dodawanie prosty sposób obsługiwać wiele sekwencji

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

wymyśliłem rozwiązanie używając BioSeqIO (i metody trunc z BioSeq z dystrybucją BioPerl użyłem również index znaleźć podciąg raczej. niż przy użyciu wyrażenia regularnego:

To rozwiązanie nie drukuje id, (linia zaczyna się od>), jeśli podciąg nie został znaleziony lub podciąg zaczyna się od pierwszej pozycji (a więc bez poprzedzających znaków).

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