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

3 回答

• 8

  您错过了质量得分的第二行（以及随后的每一行），并且还会错过其他序列行 . 为了这个和代码重用目的，处理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标头 .

  回复于 2024-05-03T06:36:56+08:00
• 4

  问题是你正在并行浏览文件，所以当一行中的行为“>”时，下一行可能不是“>” .

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

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

  回复于 2024-05-03T06:36:56+08:00
• 3

  这是一个不使用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>
  

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

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

  回复于 2024-05-03T06:36:56+08:00