提问人：Dan 提问时间：10/27/2023 最后编辑：tripleeeDan 更新时间：11/6/2023 访问量：81
从一个 fastq 文件中获取在另一个 fastq 文件中找不到的读取
Getting reads from one fastq file that are not found in another fastq file
问：
我想得到一个 fastq 文件，其中包含 fastq 读取，这些读取在 中但不在 .F1.fastqF2.fastq
此处提供了以下代码作为解决方案。
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()
答：
2赞
Timur Shtatland
10/27/2023
#1
简而言之，使用 和 UNIX 工具代替 ，如下所示：BioPythonseqtk
首先，用于从 fastq 文件中提取读取 ID。用于仅提取要保留的已读名称：grepcomm
comm -23 <( cat test2.1.fastq | paste - - - - | grep -Po '^@\K\S+' | sort ) <( cat test1.1.fastq | paste - - - - | grep -Po '^@\K\S+' | sort ) > name.lst
要按 ID 从 fastq 文件中提取读取数据，请使用 seqtk subseq：
seqtk subseq test2.1.fastq name.lst > out.fastq
它也适用于文件。fasta
或者，使用 BBMap filterbyname.sh，请参阅文档：
默认情况下，“filterbyname”会丢弃名称列表中带有名称的读取，并保留其余部分。要包含它们并丢弃其他的，请执行以下操作：
filterbyname.sh in=003.fastq out=filter003.fq names=names003.txt include=t
另请参阅：
要安装这些工具，请使用 ，特别是 ，例如：condaminiconda
conda create --channel bioconda --name seqtk seqtk
conda activate seqtk
# ... use seqtk here ...
conda deactivate
在这里，GNU 使用以下选项：
： 使用 Perl 正则表达式。
：仅打印匹配项（每行 1 个匹配项），而不是整行。grep-P-o
\K：使正则表达式引擎“保留”它之前匹配的所有内容，并且不将其包含在匹配中。具体而言，在打印匹配项时，请忽略正则表达式的前面部分。\K
另请参阅：
上述示例中使用的文件：
cat test1.1.fastq
@M01873:M01873:000000000-B4TLH:1:1101:10001:20001/1
ACGTAACCGGTTAAACCCGGGTTTAAAACCCCGGGGTTTTAAAAACCCCCGGGGGTTTTTAAAAAACCCCCCGGG
+
EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE
@M01873:M01873:000000000-B4TLH:1:1101:10002:20002/1
GGGTTTTTTAAAAAAACCCCCCCGGGGGGGTTTTTTTAAAAAAAACCCCCCCCGGGGGGGGTTTTTTTTAAAAAA
+
EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE
cat test2.1.fastq
@M01873:M01873:000000000-B4TLH:1:1101:10001:20001/1
ACGTAACCGGTTAAACCCGGGTTTAAAACCCCGGGGTTTTAAAAACCCCCGGGGGTTTTTAAAAAACCCCCCGGG
+
EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE
@M01873:M01873:000000000-B4TLH:1:1101:10002:20002/1
GGGTTTTTTAAAAAAACCCCCCCGGGGGGGTTTTTTTAAAAAAAACCCCCCCCGGGGGGGGTTTTTTTTAAAAAA
+
EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE
@M01873:M01873:000000000-B4TLH:1:1101:10003:20003/1
ACGTAACCGGTTAAACCCGGGTTTAAAACCCCGGGGTTTTAAAAACCCCCGGGGGTTTTTAAAAAACCCCCCGGG
+
EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE
@M01873:M01873:000000000-B4TLH:1:1101:10004:20004/1
GGGTTTTTTAAAAAAACCCCCCCGGGGGGGTTTTTTTAAAAAAAACCCCCCCCGGGGGGGGTTTTTTTTAAAAAA
+
EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE
cat name.lst
M01873:M01873:000000000-B4TLH:1:1101:10003:20003/1
M01873:M01873:000000000-B4TLH:1:1101:10004:20004/1
cat out.fastq
@M01873:M01873:000000000-B4TLH:1:1101:10003:20003/1
ACGTAACCGGTTAAACCCGGGTTTAAAACCCCGGGGTTTTAAAAACCCCCGGGGGTTTTTAAAAAACCCCCCGGG
+
EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE
@M01873:M01873:000000000-B4TLH:1:1101:10004:20004/1
GGGTTTTTTAAAAAAACCCCCCCGGGGGGGTTTTTTTAAAAAAAACCCCCCCCGGGGGGGGTTTTTTTTAAAAAA
+
EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE
1赞
Pierre
10/27/2023
#2
假设两个 fastq 读取具有相同的 NAME，则它们是相同的：
在 NAME 上使用 、 线性化
恢复 fastq 格式pastejointr
类似的东西（未测试）
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"
0赞
tripleee
11/2/2023
#3
就像 BioPython 文档告诉你的那样，它们改变了对象的比较方式。（这是循环遍历它时返回的内容。这实际上是一个你到底想比较什么的问题。SeqRecordSeqIO
这是一个重构，它检查实际序列是否相同，如果是，则跳过较小文件中已有的序列。
from Bio import SeqIO
chosen_seqs = [y.seq for y in SeqIO.parse("F2.fastq", "fastq")]
needed_reads = []
for x in SeqIO.parse("F1.fastq", "fastq"):
if x.seq not in chosen_seqs:
needed_reads.append(x)
with open("DIFF.fastq", "w") as output_handle:
SeqIO.write(needed_reads,output_handle, "fastq")
这被略微重构，以避免将所有信息保留在内存中;您只需要在内存中设置较小的集合，然后一次一个项目地循环访问较大的集合。
要更改选择逻辑，请从较小的文件中更改为要检查的任何特征，并相应地更改为从每个文件中提取的任何其他特征。y.seqx.seqSeqRecord
来自 Biopython 1.67（2016 年 6 月 8 日）的条目，来自 https://raw.githubusercontent.com/biopython/biopython/master/NEWS.rst：
到目前为止，SeqRecord 对象的比较使用了默认的 Python 对象 比较（它们在内存中是同一个实例吗？这可能令人惊讶，但是 比较所有属性太复杂了。截至此版本 尝试比较 SeqRecord 对象应引发异常。如果 你想要旧的行为，改用。
id(record1) == id(record2)
评论