从包含唯一分子标识符的 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

我需要一个输出文件，其中删除了给定序列（及其相应信息）的所有精确重复。也就是说，我需要删除文件中已经出现第 2 行的 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')

这会读取每条记录并计算要存储在列表中的序列的校验和。如果其序列没有遇到校验和，则每个后续记录仅添加到输出列表 ('unique')。

Answer 2

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

awk 'FNR==NR {
          if (FNR%4==2) {
              a[]++
              if (a[]>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[]++; if (a[]>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##################################

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

Removing PCR Duplicates From Fastq File Containing Unique Molecular Identifiers

python

bash

awk

sed

测试