从包含唯一分子标识符的 Fastq 文件中删除 PCR 重复

Question

我正在尝试编辑包含基因组数据和每个序列两侧的唯一分子标识符的 Fastq 文件。

前两个读取的示例如下所示：

1 @HISEQ:230:C6G45ANXX:3:1101:1395:2141 1:N:0:ACAGTGGTTGAACCTT
2 TGACGGCACTTTCTCTTCCCAACCACGTGGCTGCAGACTTCTTGCTCTCAAGTTGTCCTGACATGCTCTGAGAGCACACACAACATACATACAACACCTGGATCTGTGAATTAATTACTGCCTAGG
3 +
4 BB//<<BFBFFF<FFFFBBB<<<F/FBBB<FF/B<FFFFFFFFFFFFFFBFFFBFB/FBFFB//F//B<FFF</</BF<BBBFFFFF//B<FBFF/77F/B/BF7/FF/<BF/7FFFFBBF//B7B
5 @HISEQ:230:C6G45ANXX:3:1101:1498:2162 1:N:0:ACAGTGGTTGAACCTT
6 TGACGGCACTTTCTCTTCCCAACCACGTGGCTGCAGACTTCTTGCTCTCAAGTTGTCCTGACATGCTCTGAGAGCACACACAACATACATACAACACCTGGATCTGTGAATTAATTACTGCCTAGG
7 +
8 BBB<B<F<FFFFFFFBFFFFFFBFFFFBFF/F<FFFFBBFFFFFFFFFFBFB/BFFFFFFFFFFFBFFB/<<<FFFFFFFFFFFFFFBFFFF##################################

这些行解释如下：

1 Information
2 Sequence
3 +
4 Quality Scoring
5 Information
6 Sequence
7 +
8 Quality Scoring

我需要一个输出文件，其中已删除给定序列的所有精确重复（及其相应信息）。也就是说，我需要删除文件中已经出现第二行的 4 行块。

因此，在上面的示例中，由于第 2 行和第 6 行中的序列匹配，输出文件应该包含第 1、2、3 和 4 行，而不是 5、6、7 和 8 行。

生成的输出文件：

1 @HISEQ:230:C6G45ANXX:3:1101:1395:2141 1:N:0:ACAGTGGTTGAACCTT
2 TGACGGCACTTTCTCTTCCCAACCACGTGGCTGCAGACTTCTTGCTCTCAAGTTGTCCTGACATGCTCTGAGAGCACACACAACATACATACAACACCTGGATCTGTGAATTAATTACTGCCTAGG
3 +
4 BB//<<BFBFFF<FFFFBBB<<<F/FBBB<FF/B<FFFFFFFFFFFFFFBFFFBFB/FBFFB//F//B<FFF</</BF<BBBFFFFF//B<FBFF/77F/B/BF7/FF/<BF/7FFFFBBF//B7B

Answer 1

我认为您在 FASTQ 订单上偏离了一行。在您的示例中：

1 Junk
2 Junk
3 Junk
4 Information        +->
5 Sequence           |     these four lines constitute a single record
6 +                  |
7 Quality Scoring    +->
9 Information
10 Sequence
11 +

所以第 1-3 行（垃圾）实际上是 上一个 记录的前 3 行，而第 9-11 行是 下一个 记录的前 3 行。

无论如何，我建议您使用 BioPython 的 SeqIO 来解析您的 FASTQ 文件并进行重复数据删除。

http://biopython.org/wiki/SeqIO

一种基本的方法是：

from Bio import SeqIO
from Bio.SeqUtils.CheckSum import seguid
detected = []
unique = []
for rec in SeqIO.parse(open('inputfile.fastq', 'rU'), 'fastq'):
   cksum = seguid(rec.seq)
   if cksum not in detected:
       unique.append(rec)
       detected.append(cksum)
SeqIO.write(unique, open('deduplicated.fastq','w'), 'fastq')

这会读取每条记录并计算序列的校验和以存储在列表中。每个后续记录仅在其序列没有遇到校验和时才添加到输出列表（“唯一”）。

Answer 2

这似乎是我们循环遍历文件两次的完美案例：首先计算重复项，然后打印适当的行：

awk 'FNR==NR {
          if (FNR%4==2) {
              a[$2]++
              if (a[$2]>1) b[int(FNR/4)]=1
             }
          next}
      b[int(FNR/4)]==0' file file

这里的关键是播放文件中的 4K+2 行，并跟踪到目前为止出现了哪些行。如果是这样，我们存储K（来自4K+2），以便在文件的下一个循环中避免这些行出现在4K+0/1/2/3 形式上。

为了清楚起见，我假设第一列中的行不存在（我不知道它们是为了澄清还是真的存在）。删除它们应该是微不足道的。

测试

$ awk 'FNR==NR {if (FNR%4==2) {a[$2]++; if (a[$2]>1) b[int(FNR/4)]=1} next} b[int(FNR/4)]==0' a a
@HISEQ:230:C6G45ANXX:3:1101:1395:2141 1:N:0:ACAGTGGTTGAACCTT
TGACGGCACTTTCTCTTCCCAACCACGTGGCTGCAGACTTCTTGCTCTCAAGTTGTCCTGACATGCTCTGAGAGCACACACAACATACATACAACACCTGGATCTGTGAATTAATTACTGCCTAGG
+
BBB<B<F<FFFFFFFBFFFFFFBFFFFBFF/F<FFFFBBFFFFFFFFFFBFB/BFFFFFFFFFFFBFFB/<<<FFFFFFFFFFFFFFBFFFF##################################