从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，所以从下一对读取中，我选择第二列中的读取。

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

对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

---

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

#!/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