Biopython SeqIO：如何编写修改后的 SeqRecord header

Question

我想我会尝试使用 Biopython 来挽救合作者提供的一些轻微损坏的 fastq 文件。我只需要修改包含某个子字符串的 header 行（从 @ 开始）。但是，以下代码创建的新 fastq 文件没有变化。毫无疑问，我遗漏了一些明显的东西。

修改后的fastq SeqRecord的正确写法是什么？

import os, sys
from Bio import SeqIO

path_to_reads = sys.argv[1]
if not os.path.exists(path_to_reads + '/fixed'):
    os.mkdir(path_to_reads + '/fixed')

fwd_fastqs = [fn for fn in os.listdir(path_to_reads) if fn.endswith('_F.fastq')]
rev_fastqs = [fn for fn in os.listdir(path_to_reads) if fn.endswith('_R.fastq')]
fastq_pairs = zip(fwd_fastqs, rev_fastqs)

for fastq_pair in fastq_pairs:
    with open(path_to_reads + '/' +  fastq_pair[0], 'rU') as fwd_fastq:
        with open(path_to_reads + '/fixed/' +  fastq_pair[0], 'w') as fixed_fwd_fastq:
            fixed_fwd_records = []
            for fwd_record in SeqIO.parse(fwd_fastq, 'fastq'):
                fwd_record.name = fwd_record.name.replace('/2','/1')
                fixed_fwd_records.append(fwd_record)
            SeqIO.write(fixed_fwd_records, fixed_fwd_fastq, 'fastq')
    # ...

输入数据（两条记录，header行以@开头）：

@MISEQ01:115:000000000-A8FBM:1:1112:18038:15085/1
GATCGGAAGAGCACACGTCTGAACTCCAGTCACTGCCAATCATCTCGTATGCCGTCTTCTGCTTG
+
AAAAAAAAAF4CGGGGGAGGGFGHBHGHC5AGEGFFHGA3F355FGG223FFAEE0GCGA55BAB
@MISEQ01:115:000000000-A8FBM:1:1101:20590:9966/2
GATCACTCCCCTGTGAGGAACTACTGTCTTCACGCAGAAAGCGTCTAGCCATGGCGTTAGTATGA
+
1>A111DFBA1CFA1FFG1BFGB1D1DGFGH3GECA0ABFFG?E///DDGFBB0FEAEEFBDAB2

Answer 1

我找到了一个解决方案，但我认为它不是很明显。读取 header 行可以通过 SeqRecord.id、SeqRecord.name 和 SeqRecord.description.

中的任何一个访问

毫无疑问，它们之间存在细微差别，但我浏览了 SeqIO 文档，并没有明确提及它们。如果我添加 fwd_record.description = ''，我的脚本会按预期将 /1 的出现次数替换为 /2。

因此，工作代码：

for fastq_pair in fastq_pairs:
    with open(path_to_reads + '/' +  fastq_pair[0], 'rU') as fwd_fastq:
        with open(path_to_reads + '/fixed/' +  fastq_pair[0], 'w') as fixed_fwd_fastq:
            fixed_fwd_records = []
            for fwd_record in SeqIO.parse(fwd_fastq, 'fastq'):
                fwd_record.name = fwd_record.name.replace('/2','/1')
                fwd_record.description = ''
                fixed_fwd_records.append(fwd_record)
            SeqIO.write(fixed_fwd_records, fixed_fwd_fastq, 'fastq')

Answer 2

我不是 python 人，但我学的是生物信息学，所以我了解文件格式。我可以解释发生了什么以及为什么：

看看 BioPython Bio.SeqIO.QualityIO fastq writer code BioPython 的 SeqRecord 对象的工作方式是它有 2 个字段来存储部分定义。一个 name 和一个 description。通常人们会认为它会像 FASTA 文件一样工作，并在白色 space 上拆分 defline，左侧拆分为名称，右侧拆分为描述作为可选注释。然而，BioPython 解析器将 defline 的副本作为描述。我的猜测是这是一个 hack（连同我在下面解释的编写器代码）绕过 CASAVA 1.8 读取其中有 spaces。

当编写器写出记录时，它会检查名称和描述是否匹配，如果不匹配，则它会写出 description 行，该行被假定为 CASAVA 1.8 read I猜...

由于您只更改了 name 部分，因此匹配测试失败，因此改用未更改的描述。当您清空 description 时，作者正确地使用了 name 字段。

Biopython SeqIO：如何编写修改后的 SeqRecord header

Biopython SeqIO: how to write modified SeqRecord header

python

bioinformatics

biopython

fastq