Grep 不工作 - 用变量搜索文件

Question 1

上面的脚本将 grep 您的MKD_nsi_lib1_R1_001.fq文件 48 次。如果该文件的大小不小，您的脚本将非常慢。

它还针对 Barcodes.txt 运行cat48sed次，虽然速度不快，但不会像读取 .fq 文件 48 次那样“昂贵”（就时间和磁盘 I/O 而言）。

与其grep在同一个输入文件上运行多次，不如编写一个 awk 或 perl 脚本来一次完成您需要的操作（并且文件越大barcodes.txt，MKD_nsi_lib1_R1_001.fq您的情况就越好）。

像这样的东西：

#!/usr/bin/perl

use strict;
# %patterns is a hash where the keys are fixed-text
# strings, and the values are file-handles to 
# files opened for append.
my %patterns;

# First open the barcodes.txt file and read it into
# the %patterns hash
my $barcodes;
open($barcodes,'<','barcodes.txt') || 
  die "Couldn't open 'barcodes.txt' for read: $!\n";

while(<$barcodes>) {
  chomp; # strip the newline at the end of each line
  my $outfile = "GrepBarcode_$_.txt";
  open($patterns{$_}, ">>", $outfile) ||
    die "Couldn't open '$outfile' for append: $!\n";};
close($barcodes);

# Now process the .fq file(s) listed on the command line.
# also works with stdin.
while(<>) {
  # this assumes that the keyword is at the start
  # of the line and is followed by whitespace. This
  # is only a guess on my part, since you didn't describe
  # or provide a sample of your file.  If there's a different
  # delimiter in the input file, adjust the regex in the split
  # function.
  my ($p,undef) = split /\s+/, $_, 2;

  if (defined($patterns{$p})) {
    print { $patterns{$p} } $_;
  };
};

要运行它，您可以将其保存到一个文件（例如split-fq.pl），使用使其可执行chmod +x split-fq.pl，然后使用要处理的文件名运行它（barcodes.txt 文件被硬编码到脚本中），例如

./split-fq.pl MKD_nsi_lib1_R1_001.fq

这是为了使用固定字符串而编写的，因为它比对MKD_nsi_lib1_R1_001.fq.它只是从每个输入行中提取第一个“单词”，并检查它是否是散列中的键%patterns- 如果是，则将当前行写入关联的文件句柄。

但是，可以使用正则表达式（但速度较慢），例如

#!/usr/bin/perl

use strict;

# %patterns is a hash where the keys are pre-compiled
# regular expressions anchored to the start of line ^,
# and the values are handles to files opened for append.
my %patterns;

my $barcodes;
open($barcodes,'<','barcodes.txt') ||
    die "Couldn't open 'barcodes.txt' for read: $!\n";

while(<$barcodes>) {
  chomp;
  my $outfile = "GrepBarcode_$_.txt";
  open($patterns{qr/^$_/}, ">>", $outfile) ||
    die "Couldn't open '$outfile' for append: $!\n";
};

close($barcodes);

while(<>) {
  MATCH: foreach my $re (keys %patterns) {
    if (m/$re/) {
      print { $patterns{$re} } $_;
      last MATCH; # no need to test any more patterns against current line
    };
  };
};

这将比上面的固定文本版本慢，但仍然比grep在 shellfor循环中运行 48 次快得多 - 它只需读取 .fq 文件一次，而不是 48 次。

注意：这些只是您如何做类似事情的示例。我不知道它们是否能正确处理您的数据，因为我不知道您的文件中有什么 - 您没有提供barcodes.txt 或 .fq 文件的示例。您几乎肯定需要修改脚本以适合您的实际数据。

另请注意，完全有可能已经存在更好的用于分割 fastq 文件的工具。事实上，有一个巨大的用 Perl 编写的生物信息学脚本和工具库，位于https://bioperl.org/

如果您更喜欢 python，请参阅https://biopython.org/

当然，还有一个专门针对生物信息学问题的堆栈交换站点：https://bioinformatics.stackexchange.com/

以下版本应该适用于您提供的示例数据。

它的工作方式与第一个固定字符串版本类似（并且速度应该差不多），但它使用冒号 ( ) 作为字段分隔符将.fq 文件的每个输入行拆分为两个字段（变量$num和）。$data:

然后它使用 perl 的substr()函数将的前 5 个字母提取$data到另一个名为的变量中$start。

$start如果数组中存在值为的键%patterns，则会将当前行 ( $_) 写入关联的输出文件（中的文件句柄$patterns{$start}）。

#!/usr/bin/perl

use strict;
my %patterns;
my $barcodes;

open($barcodes,'<','barcodes.txt') || 
  die "couldn't open 'barcodes.txt' for read: $!\n";
while(<$barcodes>) {
  chomp;
  my $outfile = "GrepBarcode_$_.txt";
  open($patterns{$_},">>","$outfile") ||
    die "couldn't open '$outfile' for append: $!\n";
};
close($barcodes);

while(<>) {
  my ($num,$data) = split /:/, $_, 2;
  my $start = substr($data,0,5);

  if (defined($patterns{$start})) {
    print { $patterns{$start} } $_;
  };
};

当我运行它来测试它时，它只生成空的GrepBarcode_?????.txt输出文件 - 这是因为你的代码中没有任何行与examplefile.txt任何 5 个字母的代码匹配。我添加到 Barcodes.txt 中，它生成了包含以下内容的AAAGA文件：GrepBarcode_AAAGA.txt

$ cat GrepBarcode_AAAGA.txt 
6:AAAGAGAAATGTAATTTATACATACAGTACATATATATATGGCAGCTGTCTCCCCAAATCCTGCTCTACTGCGTCATTGTTGTGGGAATTATTCCTGGGAGGGATGCGTGAAAAATGCAAGGATATGTGCCAAGAGTACTGCAGCACTA
10:AAAGACACTGCAGATAAACCCTGTGTAATAAATACATAAAATATGTTCCAACCATTTTTATAAATTTTCTGAGTAATCTGTGTTGGATTTTCAGAGTAAGCAAATGAGAAATTAGAGTATTTGATTCCCTGTTGCTTATCCAGGACTTT

Answer