如何在 Perl 中合并两个 FASTA 文件（一个带换行符的文件）？

Question

我有以下两个 Fasta 文件：

file1.fasta

>0
GAATAGATGTTTCAAATGTACCAATTTCTTTCGATT
>1
GTTAAGTTATATCAAACTAAATATACATACTATAAA
>2
GGGGCTGTGGATAAAGATAATTCCGGGTTCGAATAC

file2.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

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

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

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

}

正确的做法是什么？

Answer 1

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

你读取数据的方式是成对的，像这样：

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
 }
}

Answer 2

这是一个不使用 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 可以）。

Answer 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 标头。