获取F1.fastq中不在F2.fastq中的fastq读取

Question

我想获得一个由 F1.fastq 中但不在 F2.fastq 中的 fastq 读取组成的 fastq 文件。

以下代码作为解决方案提供here。

needed_reads = []

reads_array = []

chosen_array = []

for x in SeqIO.parse("F1.fastq","fastq"):

        reads_array.append(x)

for y in SeqIO.parse("F2.fastq","fastq"):

        chosen_array.append(y)

needed_reads = list(set(reads_array) - set(chosen_array))

output_handle = open("DIFF.fastq","w")

SeqIO.write(needed_reads,output_handle,"fastq")

output_handle.close()

但是，当 set() 函数在 SeqIO.parse() 的输出上运行时，它不再起作用并产生“不可散列类型：SeqRecord”错误。 SeqIO.parse() 仍然根据 type() 函数生成一个列表类型对象。然而，该列表似乎是不可散列的。如何解决这个问题？

Answer 1

简单来说，使用

BioPython

和 UNIX 工具代替

seqtk

，如下所示：

如何从fastq文件中删除读取列表？

首先，使用

grep

从 fastq 文件中提取读取 ID。使用

comm

仅提取要保留的已读名称。

要按 ID 从 fastq 文件中提取读数，请使用

seqtk subseq

。

提取文件name.lst中名称的序列，每行一个序列名称：

seqtk subseq in.fq name.lst > out.fq

它也适用于

fasta

文件。

或者，使用

BBMap filterbyname.sh

，请参阅 docs:

默认情况下，“filterbyname”会丢弃名称列表中名称的读取，并保留其余内容。要包含它们并丢弃其他的，请执行以下操作：

filterbyname.sh in=003.fastq out=filter003.fq names=names003.txt include=t

另请参阅：

• 使用 Seqtk 从 FASTQ 中提取特定序列
• 如何提取给定 Fastq Id 和 Fastq 文件的 Fastq 序列
• 如何使用 ID 列表提取 Fastq 中的读取子集？

要安装这些工具，请使用

conda

，特别是

miniconda

，例如：

conda create --channel bioconda --name seqtk seqtk
conda activate seqtk
# ... use seqtk here ...
conda deactivate

参考资料：

• seqtk

• conda

• conda create