从一个 fastq 文件中获取在另一个 fastq 文件中找不到的读取内容

Question

我想获得一个由 fastq 读取组成的 fastq 文件，这些读取位于

F1.fastq

但不在

F2.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()

但是，当

"unhashable type: SeqRecord"

函数在

set()

的输出上运行时，它不再起作用并产生

SeqIO.parse()

错误。

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

Answer 2

如果两个 fastq 读取共享相同的 NAME，则它们是相同的：

用

paste

进行线性化，

join

NAME， 使用

tr

恢复 fastq 格式

类似（未测试）

 join -t $'\t' -v1 -1 1 -2 1 \
      <(cat F1.fastq | paste - - - - | sort -t $'\t' -k1,1) \
      <(cat F2.fastq | paste - - - - | sort -t $'\t' -k1,1) |\
 tr "\t" "\n"