cebach561 cebach561 - 5 months ago 19
Perl Question

basic regex and string manipulation for DNA analysis using perl

I am new to perl and would like to do what I think is some basic string manipulation to DNA sequences stored in an rtf file.

Essentially, my file reads (file is in FASTA format):

>LM1
AAGTCTGACGGAGCAACGCCGCGTGTATGAAGAAGGTTTTCGGATCGTAA
AGTACTGTCCGTTAGAGAAGAACAAGGATAAGAGTAACTGCTTGTCCCTT
GACGGTATCTAACCAGAAAGCCACGGCTAACTACGTGCCAGCAGCCGCGG
TAATACGTAGGTGGCAAGCGTTGTCCGGATTTATTGGGCGTAAAGCGCGC
GCAGGCGGTCTTTTAAGTCTGATGTGAAAGCCCCCGGCTTAACCGGGGAG
GGTCATTGGAAACTGGAAGACTGGAGTGCAGAAGAGGAGAGTGGAATTCC
ACGTGTAGCGGTGAAATGCGTAGATATGTGGAGGAACACCAGTGGCGAAG
GCGACTCTCTGGTCTGTAACTGACGCTGAGGCGCGAAAGCGTGGGGAGCA
AACAGGATTAGATACCCTGGTAGTCCACGCCGT


What I would like to do is read into my file and print the header (header is >LM1) then match the following DNA sequence
GTGCCAGCAGCCGC
and then print the preceding DNA sequence.

So my output would look like this:

>LM1
AAGTCTGACGGAGCAACGCCGCGTGTATGAAGAAGGTTTTCGGATCGTAA
AGTACTGTCCGTTAGAGAAGAACAAGGATAAGAGTAACTGCTTGTCCCTT
GACGGTATCTAACCAGAAAGCCACGGCTAACTAC


I have written the following 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);


The problem is that my program reads the file line by line and the output I am receiving is the following:

>LM1
GACGGTATCTAACCAGAAAGCCACGGCTAACTAC


Basically I don't know how to assign the entire DNA sequence to my $dna variable and ultimately don't know how to avoid reading the DNA sequence line by line. Also I am getting this warning:
Use of uninitialized value $dna in pattern match (m//) at stacked.pl line 14, line 1113.

If anyone could give me some help with writing better code or point me in the correct direction it would be much appreciated.

Answer

Using the 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

You can process a file with multiple entries like so:

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";
}
Comments