3

我有两个以下 Fasta 文件：

文件1.fasta

>0
GAATAGATGTTTCAAATGTACCAATTTCTTTCGATT
>1
GTTAAGTTATATCAAACTAAATATACATACTATAAA
>2
GGGGCTGTGGATAAAGATAATTCCGGGTTCGAATAC

文件 2.qual

>0
40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40
40 40 40 40 40 40 40 40 15 40 40
>1
40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 20 40 40 40
40 40 40 40 40 40 40 40 40 40 40
>2
40 40 40 40 7 40 40 5 40 40 40 40 40 40 40 40 37 13 31 20 15 40 10 11 4
40 8 3 29 10 19 18 40 19 15 5

请注意每个 fasta 标头的“qual”文件中的换行符 - 用“>”标记。两个文件的文件头 ('>') 的数量相同。数字质量的数量=序列长度。

我想要做的是附加这两个文件产生：

GAATAGATGTTTCAAATGTACCAATTTCTTTCGATT  40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 15 40 40
GTTAAGTTATATCAAACTAAATATACATACTATAAA  40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 20 40 40 40 40 40 40 40 40 40 40 40 40 40 40
GGGGCTGTGGATAAAGATAATTCCGGGTTCGAATAC  40 40 40 40 7 40 40 5 40 40 40 40 40 40 40 40 37 13 31 20 15 40 10 11 4 40 8 3 29 10 19 18 40 19 15 5

但不知何故，我下面的代码无法正确执行？特别是 'qual' 文件中每个条目的第二行不会被打印出来。

use strict;
use Data::Dumper;        
use Carp;
use File::Basename;      

my $fastafile = $ARGV[0] || "reads/2039F.2.fasta"; 
my $base      = basename( $fastafile, ".fasta" );
my $qualfile  = "reads/" . $base . ".qual";
print "$qualfile\n";

open SEQ, '<', $fastafile or die $!; #Seq
open PRB, '<', $qualfile or die $!; #quality


while (my $seq = <SEQ>) {
     my $qual = <PRB>;
     chomp($seq);
     chomp($qual);

     if ($seq =~ /^>/ || $qual =~ /^>/) {
         next;
     }
     else {
         print "$seq\t$qual\n";      
     }

}

正确的方法是什么？

4

3 回答 3

8

问题是您正在并行推进文件，因此当一个文件中的行是“>”时，下一个文件中可能不是“>”。

您读取数据的方式是成对的，如下所示：

1：>0
2：>0
1：GAATAGATGTTTCAAATGTACCAATTTCTTTCGATT
2：40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40
1：>1
2：40 40 40 40 40 40 40 40 15 40 40
1：GTTAAGTTATATCAAACTAAATATACATACTATAAA
2：>1
1：>2
2：40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 20 40 40 40
1：GGGGCTGTGGATAAAGATAATTCCGGGTTCGAATAC
2：40 40 40 40 40 40 40 40 40 40 40
1：EOF
2：>2
1：EOF
2：40 40 40 40 7 40 40 5 40 40 40 40 40 40 40 40 37 13 31 20 15 40 10 11 4
1：EOF
2：40 8 3 29 10 19 18 40 19 15 5

应用您的循环规则的同一组数据将执行此操作：

1：GAATAGATGTTTCAAATGTACCAATTTCTTTCGATT
2：40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40
1：GGGGCTGTGGATAAAGATAATTCCGGGTTCGAATAC
2：40 40 40 40 40 40 40 40 40 40 40

因此，您需要将循环逻辑分开或找到使文件匹配的方法。

这是分离搜索的尝试，但我还没有测试过。

fileIO: {
  while( 1 ){ 
   my $seq; 
   my $qual  = q{};
   while( 1 ){ 
     $seq = <SEQ>; 
     last fileIO if not $seq;  # stop at end of file
     last if $seq !~ /^>/; 
  }
  while( 1 ){ 
     my $qual_in = <PRB>;
     last fileIO if not $qual_in; # stop at end of file 
     last if $qual_in =~ /^>/ and $qual ne q{}; 
     next if $qual_in =~ /^>/ and $qual eq q{}; 
     $qual .= $qual_in;
  }
  print "$seq \n $qual \n";

 }
}

更新

我将上面的代码重新分解为一个函数，它将根据需要从任意文件句柄中读取一个块，它似乎可以根据需要工作。当然请注意，我在这里尝试了一些我一直打算用于实际操作的技巧。

use strict;
use warnings;

# 
#  readUntilNext( $fileHandle, \$scalar_ref ); 
#
#  returns 0 when nothing could be read from the fileHandle. 
#  otherwise returns 1; 
#

sub readUntilNext {
    my ($fh)            = shift;
    my ($output)        = shift;
    my ($output_buffer) = '';
    while (1) {
        my $line = <$fh>;
        if ( !$line ) { # No more data
            # No data to flush to user, return false.
            return 0 if $output_buffer eq q{};
            last;  # data to  flush to user, loop exit. 
        }
        if ( $line =~ /^>/ ) {
            # Didn't get anything, keep looking. 
            next if $output_buffer eq q{};
            # Got something, flush data to user. 
            last;
        }
        chomp($line);
        $output_buffer .= $line;
    }
    # Data to flush to user 
    # Write to the scalar-reference 
    $$output .= $output_buffer;
    return 1;
}

open my $m, '<', 'a.txt';
open my $n , '<', 'b.txt';
# Creates 2 scalar references every loop, and only loops as long 
# as both files have data. 
while ( readUntilNext( $m, \my $seq ) && readUntilNext( $n, \my $qual ) ) {
    print "$seq\t$qual\n";
}

上面的代码，经过测试，完全符合您的要求。

请注意\我的东西

while( readUntilNext( $m, \my $seq ) ) { 
}

基本相同

my $seq; 
while( readUntilNext( $m, \$seq ) ) { 
}

除了前者每次都创建一个新的标量，保证相同的值不会被连续循环看到；

所以它变得更像：

while( 1 ){ 
 my $seq; 
 last if not readUntilNext($m, \$seq);
 do { 
    # loop body here
 }
}
于 2009-04-10T06:46:15.257 回答
4

这是一个不使用 perl 而是使用普通 shell 命令的解决方案：

prompt>grep -v '^>[0-9]' file1.fasta > tmp1
prompt>(tr '\012' ' ' < file2.qual; echo) | sed 's/>[0-9]* /\n/g' | sed 1d > tmp2
prompt>paste tmp1 tmp2
GAATAGATGTTTCAAATGTACCAATTTCTTTCGATT    40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 15 40 40
GTTAAGTTATATCAAACTAAATATACATACTATAAA    40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 20 40 40 40 40 40 40 40 40 40 40 40 40 40 40
GGGGCTGTGGATAAAGATAATTCCGGGTTCGAATAC    40 40 40 40 7 40 40 5 40 40 40 40 40 40 40 40 37 13 31 20 15 40 10 11 4 40 8 3 29 10 19 18 40 19 15 5
prompt>

我搜索了很多年的粘贴命令（知道“这是一个超级基本的操作，一定有人已经实现了一些东西来解决这个问题”）。

第二个命令行首先将所有换行符转换为空格，然后添加 echo 命令以将最终换行符添加到输入中（因为 sed 将忽略缺少 EOL 的行），从而将所有输入行合并为一行，然后 sed 命令再次分裂（可移植性说明：并非所有 sed 程序都可以使用任意行长度，但 GNU sed 可以）。

于 2009-04-10T15:54:07.843 回答
3

您错过了质量分数的第二行（以及随后的每一行），并且还会错过其他序列行。出于此和代码重用目的，处理 FASTA 序列的方法是作为整个条目/记录：

local $/ = "\n>";
while (my $seq = <SEQ>) {
     my $qual = <PRB>;
     chomp($seq);  $seq =~ s/^>*.+\n//;  $seq =~ s/\n//g;
     chomp($qual);  $qual =~ s/^>*.+\n//;  $qual =~ s/\n/ /g;

     print "$seq\t$qual\n";      

}

您还可以在第一次替换时轻松捕获 FASTA 标头。

于 2009-04-11T21:26:07.880 回答