从测序 Fastq 文件中删除很少出现的行

Question

我有一个文本文件，其中包含 4 行长组的实验数据。我想删除那些很少见的数据点。所以下面是前 8 行文件格式，这重复了数千行（文件中不存在行号）：

1 Info_1
2 Sequence_1
3 Info_1
4 Info_1
5 Info_2
6 Sequence_2
7 Info_2
8 Info_2
9 Info_3
10 Sequence_3
11 Info_3
12 Info_3

因此，第 1-4 行包含序列 1 的信息，第 5-8 行包含序列 2 的信息，第 9-12 行包含序列 3 的信息...等等。在某些情况下，通常会删除包含完全唯一或发现次数少于 3 次的序列的任何一组四行。

我想做的是将第 2 行与第 6、10、14、18 行进行比较，如果发现大于 3 次，则什么也不做。如果找到 3 次或更少，则删除第 1-4 行和包含匹配序列的每组 4 行。然后运行对文件中的每一行进行相同的比较。

因此，如果在上述文件中，序列 1 和序列 3 匹配，并且因为该序列只重复了 < 3 次，删除每组四行，生成的文件应如下所示：

1 Info_2
2 Sequence_2
3 Info_2
4 Info_2

这是我的开头：

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

这不会删除所有发现次数少于 3 次的行。我将不胜感激任何帮助。谢谢。

根据要求，这是一个真实的可测试示例：输入：

Info1
AAGC
Info1
Info1
Info2
AACT
Info2
Info2
Info3
AAGC
Info3
Info3
Info4
AAGC
Info4
Info4
Info5
AACT
Info5
Info5

因为 AAGC 出现 >= 3 次但 AACT 出现 <3 次输出应该是：

Info1
AAGC
Info1
Info1
Info3
AAGC
Info3
Info3
Info4
AAGC
Info4
Info4

希望这有助于澄清。

Answer 1

这道题最好用高级编程语言和专门的生物信息库，我用的是python和biopython

from Bio import SeqIO

input_file = open("input.fastq")

#to store sequence in dictionary using dna sequences as keys
#it consume RAM (if fastq is very very large, it could be a problem)
dict_seq = {}
for seq in SeqIO.parse(input_file, "fastq"):
  if not str(seq.seq) in dict_seq:
    dict_seq[str(seq.seq)] = []
  dict_seq[str(seq.seq)].append(seq)

#filter sequences
list_filter = [ dna for dna in dict_seq.keys() if len(dict_seq[dna]) >= 3]

#print output in fastq format
output_file = open("filter.fastq", "w")
for dna in list_filter:
  for seq in dict_seq[dna]:
    output_file.write(seq.format("fastq"))

output_file.close()

input.fastq

@EAS54_6_R1_2_1_413_324
CCCTTCTTGTCTTCAGCGTTTCTCC
+
;;3;;;;;;;;;;;;7;;;;;;;88
@EAS54_6_R1_2_1_540_792
TTGGCAGGCCAAGGCCGATGGATCA
+
;;;;;;;;;;;7;;;;;-;;;3;83
@EAS54_6_R1_2_1_443_348
GTTGCTTCTGGCGTGGGTGGGGGGG
+EAS54_6_R1_2_1_443_348
;;;;;;;;;;;9;7;;.7;393333
@EAS54_6_R1_2_1_413_789
CCCTTCTTGTCTTCAGCGTTTCTCC
+
;;3;;;;;;;;;;;;7;;;;;;;88
@EAS54_6_R1_2_1_413_329
CCCTTCTTGTCTTCAGCGTTTCTCC
+
;;3;;;;;;;;;;;;7;;;;;;;88

filter.fastq

@EAS54_6_R1_2_1_413_324
CCCTTCTTGTCTTCAGCGTTTCTCC
+
;;3;;;;;;;;;;;;7;;;;;;;88
@EAS54_6_R1_2_1_413_789
CCCTTCTTGTCTTCAGCGTTTCTCC
+
;;3;;;;;;;;;;;;7;;;;;;;88
@EAS54_6_R1_2_1_413_329
CCCTTCTTGTCTTCAGCGTTTCTCC
+
;;3;;;;;;;;;;;;7;;;;;;;88

Answer 2

$ cat tst.awk
{
    numRecs += (NR%4 == 1 ? 1 : 0)
    rec[numRecs,(NR-1)%4 + 1] = [=10=]
    cnt[[=10=]]++
}
END {
    for (recNr=1; recNr<=numRecs; recNr++) {
        if (cnt[rec[recNr,2]] > 2) {
            for (fldNr=1; fldNr<=4; fldNr++) {
                print rec[recNr,fldNr]
            }
        }
    }
}

$ awk -f tst.awk file
Info1
AAGC
Info1
Info1
Info3
AAGC
Info3
Info3
Info4
AAGC
Info4
Info4

Answer 3

正在扩展您的脚本：

$ awk '
    (NR==FNR){ # First pass, record the count.
        if(FNR%4==2){a[]++}
        next;
    }

    (FNR%4){ # remember the group of 4 lines.
        r[FNR%4]=[=10=];
        next;
    }

    { # On every 4th line, check if the count is >= 3 & print the entire 4 line group.
        if (a[r[2]] >= 3) print r[1] "\n" r[2] "\n" r[3] "\n" [=10=] 
    }' inputFile inputFile

从测序 Fastq 文件中删除很少出现的行

Remove Rarely Occurring Lines From Sequencing Fastq File

bash

awk

fastq