5

我是 perl 新手，想做一些我认为是对存储在 rtf 文件中的 DNA 序列的基本字符串操作。

本质上，我的文件读取（文件为 FASTA 格式）：

>LM1
AAGTCTGACGGAGCAACGCCGCGTGTATGAAGAAGGTTTTCGGATCGTAA
AGTACTGTCCGTTAGAGAAGAACAAGGATAAGAGTAACTGCTTGTCCCTT
GACGGTATCTAACCAGAAAGCCACGGCTAACTACGTGCCAGCAGCCGCGG
TAATACGTAGGTGGCAAGCGTTGTCCGGATTTATTGGGCGTAAAGCGCGC
GCAGGCGGTCTTTTAAGTCTGATGTGAAAGCCCCCGGCTTAACCGGGGAG
GGTCATTGGAAACTGGAAGACTGGAGTGCAGAAGAGGAGAGTGGAATTCC
ACGTGTAGCGGTGAAATGCGTAGATATGTGGAGGAACACCAGTGGCGAAG
GCGACTCTCTGGTCTGTAACTGACGCTGAGGCGCGAAAGCGTGGGGAGCA
AACAGGATTAGATACCCTGGTAGTCCACGCCGT

我想做的是读入我的文件并打印标题（标题>LM1）然后匹配以下DNA序列GTGCCAGCAGCCGC，然后打印前面的DNA序列。
所以我的输出看起来像这样：

>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 序列分配给我的 $dna 变量，最终也不知道如何避免逐行读取 DNA 序列。我也收到此警告：在stacked.pl 第14 行第1113 行的模式匹配（m//）中使用未初始化的值$dna。

如果有人可以帮助我编写更好的代码或指出正确的方向，我将不胜感激。

4

4 回答 4

3

使用pos 函数：

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";
}
于 2013-03-04T21:08:06.407 回答
2

您的 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. else您可以使用语句更新您的 while 循环。因此，如果它是标题行，则打印它，否则，我们将其添加到$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变量将所有的序列都作为一个。

编辑：添加简单的方法来处理多个序列

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";
}
于 2013-03-04T21:10:52.633 回答
2

我想出了一个使用BioSeqIO的解决方案（以及来自BioPerl发行版的BioSeqtrunc的方法。我还使用索引来查找子序列，而不是使用正则表达式。

如果未找到子序列或子序列从第一个位置开始（因此没有前面的字符），此解决方案不会打印出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
于 2013-03-05T01:57:18.683 回答
0

将整个文件读入内存然后查找正则表达式

while(<FASTA>) {
    chomp($_);
    if ($_ =~  m/^>/ ) {
        my $header = $_;
        print "$header\n";
    } else {
    $dna .= $_;
    }
}
if ($dna =~ /(.*?)GTGCCAGCAGCCGC/) {
    print $1;
}
于 2013-03-04T21:05:23.303 回答