根据配对长度处理 FASTQ 文件

Question

下面的文件是双端fastq文件的两个mate，我想根据长度把每个fastq分开。

mate1.fq:

@SRR127.1
TGGTTATGATGTTTGTGTAGGAATAGAAATTTTGATTAAGATATTAGTGAAATTTGAATGTAGTTTATTTGGAAGTTATGGAGAGTTTATATTGTATTTATGTTTATTGTTGTAGATTTATATTTATGTGTATATATTAGTTTTTTTGTGT
+
ABAAAF4FFFFFGGGGGGFFGGFGHGFGHHHHHGGCFFGHHHHH5FDBED55DGGFEGFHHHGBHDDHHHFF3AB3FFG5CBGBEF5BD5DGFEGHFAGAFEDGHGFHHGHGEFFGFGGHFEGHHFHGBEBGHHHHGHBHHFHHGGFGHH2
@SRR127.2
TATGGTAAGAAAATTGAAAATTATAAAAAATGAAAAATGTTTATTTGATGATTTGAAAAATGATGAAATTATTGAAAAATGTGAAAAATGAGAAATGTATATTGTAGGATTTGGAATATGGTGAGATAAATGAAAATTATAGTAAATG
+
AABAA5@D4@5CFFCA55FFGGHDGFHFFCC45DGFA2FA5DD55AAAA55DDBDEDDBGGFF5BA5DDABF5D5B5FF1ADFB5EDGHFG5@BFBD55D5FFB@@5@GBGEFBGHHGB@DBBFHFBDG3B43FFH@FGFHH?FHHHH

mate2.fq:

@SRR127.1
ACCTATAAAAAAACCATATCAATAACTATAAAATCTTTATAAAATCCCACCCAATTAAAAAAAAATAAATTAATACATATAAAACCTTAAACACATAAAACATAATCACATACTATATAAACAATTACTATCACTACTAAACACCTAATA
+
>AA?AF13B@D@1EFCGGGFFG3EBGHHHBB2FGHHGHGFDGHHDFEGFHGGGHG1FFF1GGCGGGBGHHHHHFHHHHFHEGGFHF0BD1FGHHAGEGHFHHHFGGFHHGHHHFHHGGFHBGHFED1FBGFGFHDGHGHFGG1GB0GFHH
@SRR127.2
CTATTTCTCATTTTTTTATAATTTTCAATTCTCTTACCATATTCCACATCCTACACTAAACATTTCTAAATTTTCCACCTTTTTCTATTTTTCTCACCATATTTCATATCCTAAAAAACATATTCCTCATTTACTATAATTTTCAATTATC
+
11>>AFFDFF3@FFF?EFFGFBGHFDFA33D2FF2GGHFE12DD221AF1F1E1BG1GGBFBGGEGHDAABGAGDFABGG1BBDF12A2@2BG@2@DEFFF2B2@2222BB2211FGEE/11@22B2>1B22F2>GBGBD22BGD2>2B22

我编写了以下代码来执行此操作，但我只收到第二个文件 (mate2.fq) 的奇怪错误，而它们都具有 151 bp 读取。

#!/usr/bin/perl

use strict;
use warnings;

my @fh;

my $file_name = $ARGV[0];
my $infile    = $ARGV[1];

#convert every 4-line fastq to 1-line
open(FH, "cat '$infile' | awk '{printf \"%s%s\",$0,(NR%4?FS:RS)}' | ");

while (<FH>) {
  chomp;

  my @line = split(/\s+/, $_);
  my $len  = length($line[1]);

  if ($len >= 100) {

    #print $len,"\n",$_,"\n";
    push @fh, $len;

    if (not defined $fh[$len]) {
      open $fh[$len], '>', "$file_name\_$len";
    }
    print { $fh[$len] } (join("\n", @line), "\n");
  }

}

错误:

Can't use string ("151") as a symbol ref while "strict refs" in use at

如何处理这些文件？

Answer 1

此错误具体意味着您正在做一些需要参考的事情，但没有得到参考。

行：

print {$fh[$len]} (join("\n",@line),"\n");

正在显式打印到一个文件句柄 - 来自一个名为 @fh 的文件句柄列表。

这一行：

push @fh, $len;

将向该列表中插入一个数值。（大概 $line[1] 是 151 个字符长）。所以你实际上是在尝试：

 print {151} (join("\n",@line),"\n");

希望这很明显 - 只是行不通。您看起来像是在尝试打开一个文件句柄，并将其插入到一个数组中：

open $fh[$len], '>', "$file_name\_$len";

我可以建议您最好为此使用哈希吗？否则你会得到一个充满空元素的数组，其中一个被填充。

哪里可以代替：

#further up:
my %fh; 


#and then
open ( $fh{$len}, ">", "$file_name\_$len" ) or warn $!;

别忘了在最后关闭文件句柄：

foreach my $key ( keys %fh ) {
   close ( $fh{$key} );
}

我也建议而不是：

open( FH, "cat '$infile' | awk '{printf \"%s%s\",$0,(NR%4?FS:RS)}' | " );

你可能最好在 perl 中处理它，因为你所做的只是使用外部二进制文件解析文件。（并使用词法文件句柄：`open ( $input, "-|, "cat '$infile' | awk '{printf \"%s%s\",\$0,(NR%4?FS:RS)}'" ）或 warn $!; ）

Answer 2

如您所见，您的问题是由于伪造的 push 将整数值添加到 @fh 数组的末尾。我假设您的目标是将数组扩展到足够长以添加新的文件句柄。你可以通过分配给 $#fh 来做到这一点，所以你会写 $#fh = $len if $#fh < $len;然而这是不必要的，因为当你简单地分配给数组末尾的元素时，Perl 会自动为你扩展数组

我对你的程序有一些评论，希望你觉得有用

shell out 到 awk 命令是不必要和浪费的。 Perl 完全有能力完成所有 awk 可以做的事情
如果您发现自己在写 split /\s+/, $_，那么您几乎肯定是指 split：默认行为是 split ' ', $_。如果您使用 /\s+/ 作为模式并且恰好在您要拆分的字符串上有前导白色 space，那么 split 将 return 一个空字符串作为第一个项目字段列表。如果您改用 ' '（字面上的单个 space，而不是模式 / /），则不会发生这种情况。实际上，split ' ' 等同于 /\S+/g
在字符串中插入变量值时，如果后面有可能是标识符一部分的字符，则将标识符放在大括号内通常会更整洁。所以 "${file_name}_$len" 而不是 "$file_name\_$len"

这就是我编写代码的方式。它将输入记录累积到 $line 中，直到添加了四个记录，然后像以前一样处理该行。

#!/usr/bin/perl

use strict;
use warnings;

my ($file_name, $infile) = @ARGV;

open my $in_fh, '<', $infile or die $!;
my $line;

my @fh;
while ( <$in_fh> ) {
  chomp;
  $line .= $_;

  if ( $. % 4 == 0 or eof ) {

    my @line = split ' ', $line;
    my $len  = length $line[1];
    next if $len < 100;

    open $fh[$len], '>', "${file_name}_$len" unless $fh[$len];
    print { $fh[$len] } "$_\n" for @line;

    $line = undef;
  }
}

根据配对长度处理 FASTQ 文件

Processing FASTQ files based on mate pair length

perl

bioinformatics

bioperl

fastq

sequencing