从 fastq 文件中的一对读取中选择一个随机读取

Question

我有一个关于从采样的对端 fastq 文件中随机选择读取的问题。我阅读了一些有关这种方式的主题，但没有一个可以解决我的问题，即： 我有两个 fastq 文件 R1.fastq 和 R2.fastq。我想要实现的是对这些文件进行随机采样，并从每对采样的读取中随机选择一个读取。

到目前为止，我所做的是......

我使用 seqtk 对我的文件进行了采样：

seqtk sample -s100 R1.fastq 10000 > R1_sample.fastq
seqtk sample -s100 R2.fastq 10000 > R2_sample.fastq

然后我按序列 ID 对每个文件进行排序，如下所示：

paste - - - - < R1_sample.fastq | sort -k1 -t " " | tr "\t" "\n" > R1_sample_sorted.fastq

我对 R2_sample.fastq 做了同样的事情。然后我合并了两个排序的文件，这样 R1 将在一个列中，R2 在第二列：

pr -mts R1_sample_sorted.fastq R2_sample_sorted.fastq > merged.fastq

文件如下所示：

@D3YGT8Q1:297:C7T4RACXX:3:1101:1000 @D3YGT8Q1:297:C7T4RACXX:3:1101:1000
TGATGTTTGGATGTAAAGTGAAATATTAGTTGGCG AGCTTTCCTCACTATCTGCTTCATCCGCCAACTAA
+                   +
BBBFFFFFFFFFFFIFFIFFIIIIFIIIFIIFIII B0<FFFFFFFFFFIIIIIIIIIIIIIIIIIIIIII
@D3YGT8Q1:297:C7T4RACXX:3:1101:1000 @D3YGT8Q1:297:C7T4RACXX:3:1101:1000
CCTCCTAGGCGACCCAGACAATTATACCCTAGCCA TGTTTAAGGGGTTGGCTAGGGTATAATTGTCTGGG
+                   +
BBBFFFFFFFFFFIIIIIIIIIIIIIIIIIIIIII BBBFFFFFFFFFFIIIIIIIIBFFIIIIIIIIIII
@D3YGT8Q1:297:C7T4RACXX:3:1101:1000 @D3YGT8Q1:297:C7T4RACXX:3:1101:1000
TTCTATTTATTACCTCAGAAGTTTTTTTCTTCGCA GTAAAAGGCTCAGAAAAATCCTGCGAAGAAAAAAA
+                   +
BBBFFFFFFFFFFIIIIIIIIFIIFIIIFIIIIII BBBFFFFFFFFFFIIIIIIIIIIIIIIIIIIIIII

现在我想从每对中随机选择一个读取。我最初的想法是使用 shuf 从 1-2 范围内获取随机数：

shuf -i1-2 -n1

然后以某种方式选择与我从 shuf 获得的数字相对应的读取。例如，在第一次迭代中我得到 1，所以我从第 1 列中选择读取，在第二次迭代中我得到 2，所以从下一对读取中我选择第二列中的读取等。

我被困在这里了。所以我的问题是，有没有一种巧妙的方法可以做到这一点？也许用 awk 或其他方法？任何帮助将不胜感激。

对 Ashafixs 回答的评论：
感谢您的回复，并对造成的巨大延迟表示歉意...！ 我已经测试了您的解决方案，它们似乎都有缺陷。
对于第一个脚本，我构建了测试 fastq 文件 R1 和 R2，每个文件包含 6 个读取。运行脚本后，我希望它以正确的顺序（ID、seq、desc、qual）输出 6 个读取（24 行），但作为从 R1 或 R2 文件中随机选择的一组读取。我从脚本中得到的是：

@D3YGT8Q1:297:C7T4RACXX:3:1101:10002:27381 2:N:0:ATGCTCGTTCTCTCGT
AGCTTTCCTCACTATCTGCTTCATCCGCCAACTAATATTTCACTTTACATCCAAACATCAAGATC
+
B0<FFFFFFFFFFIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIFIFIFIIIIIIIIII
@D3YGT8Q1:297:C7T4RACXX:3:1101:10004:50631 2:N:0:ATGCTCGTTCTCTCGT
@D3YGT8Q1:297:C7T4RACXX:3:1101:10007:32152 1:N:0:ATGCTCGTTCTCTCGT
GTAAGGTTAGGAGGGTGTTAATTATTAAAATTAAGGCGAAGTTTATTACTCTTTTTTGAATGTTG
+
BBBFFFFFFFFFFIIBFFIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIFFFFFFFF

可以看到输出不正确。第二次阅读缺少三行，总共应该有六个阅读而不是三个。此外，每次我运行脚本时，它都会输出不同数量的读取。

对于第二个脚本，我输入了一个合并的 fastq 文件，如上所述。输出看起来类似于第一个脚本输出：

@D3YGT8Q1:297:C7T4RACXX:3:1101:10002:27381 2:N:0:ATGCTCGTTCTCTCGT
AGCTTTCCTCACTATCTGCTTCATCCGCCAACTAATATTTCACTTTACATCCAAACATCAAGATC
+
B0<FFFFFFFFFFIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIFIFIFIIIIIIIIII
@D3YGT8Q1:297:C7T4RACXX:3:1101:10004:50631 2:N:0:ATGCTCGTTCTCTCGT
@D3YGT8Q1:297:C7T4RACXX:3:1101:10004:50631 2:N:0:ATGCTCGTTCTCTCGT
TGTTTAAGGGGTTGGCTAGGGTATAATTGTCTGGGTCGCCTAGGAGGAGATCGGAAGAGCGTCGT
+
BBBFFFFFFFFFFIIIIIIIIBFFIIIIIIIIIIIFFFIIIIIIFIIIIIFIIIFFFFFFFFFFF
@D3YGT8Q1:297:C7T4RACXX:3:1101:10004:88140 1:N:0:ATGCTCGTTCTCTCGT
ACTGTAACTTAAAAATGATCAAATTATGTTTCCCATGCATCAGGTGCAATGAGAAGCTCTTCATC
+
BBBFFFFFFFFFFIIIIIIIIIIFIIIIIIFIIIIIIIIIIIIIFIIIIIIIIIIIIIIIIIIII
@D3YGT8Q1:297:C7T4RACXX:3:1101:10007:32152 2:N:0:ATGCTCGTTCTCTCGT
CTAGTTTTGACAACATTCAAAAAAGAGTAATAAACTTCGCCTTAATTTTAATAATTAACACCCTC
+
BBBFFFFFFFFFFIIIIIIIIIIIIIIFFFIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIFIII  

但这一次我总是得到五次阅读。还缺一个。并且第二个和第三个读头是一样的。这不应该发生。

Answer 1

您可以尝试以下脚本（它也可以作为一个衬里使用）。首先它从你的第一个 fastq 文件中获取所有标题，然后它随机选择一个 fastq 文件并从中返回 4 行。 请注意：这仅适用于两个文件在相同位置具有相同标题的情况。

#!/bin/bash
headers=$(grep @ R1_sample.fastq)
var=1
for line in $headers ; do 
  r=$(shuf -i1-2 -n1)
  tail -n +$var "R$r"_sample.fastq |  grep -m 1 -A 4 $line
  var=$((var+4))
done

---

或者，您可以扩展合并并选择列方法。 cut 用于从合并输出中删除随机列。

#!/bin/bash
headers=$(grep @ merged.fastq)
var=1

for line in $headers ; do 
  r=$(shuf -i1-2 -n1)
  tail -n +$var merged.fastq | grep -m 1 -A 4 $line | cut -d$'\t' -f$r
  var=$((var+4))
done