我正在尝试编辑包含基因组数据和每个序列两侧的唯一分子标识符的 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 行 block ,其中第二行已经出现在文件中。
所以在上面的例子中,因为序列在第 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
最佳答案
这似乎是我们遍历文件两次的完美案例:首先计算重复项,然后打印适当的行:
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##################################
关于python - 从包含唯一分子标识符的 Fastq 文件中删除 PCR 重复项,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/29477233/