Python 解析 FastQ 文件 - 序列和质量得分修剪

Question

我正在尝试根据修剪后的质量得分修剪序列。由于我对 python 比较陌生，所以我一直在寻找一些可以解决问题的简单方法。

我有以下基于 phred 分数转换为数值的质量分数：

1816 28 32 32 26 164 30 32 1816 1618 216 22 1616 216 218 20 28 1816 1816 24 1816 1620 20 24 28 28 28 216 26 28 1622 216 28 24 18 24 28 1622 30 222 24 18 218 216 26 26 218 28 1624 24 16 26 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2

下面是对应的序列：

ACCGAGCCGAAGGAGACCGCATTCACCCGGATGCCCTTCGAGGCCAGCGCCACCGCCATCGACCGCGTCATCTGCTCCACCGCGGCCCAGCTGATGGAACA

本质上，我希望能够删除/剪切/修剪质量分数末尾的所有值“2”，比较两个字符串并删除/剪切/修剪序列末尾的相应序列核苷酸.

我尝试将质量分数分成两部分，并使用它来打印序列的第一部分：

cutquality = actualquality.split(" 2 ",1)
newquality = cutquality[0]
lenofseq = len(cutquality[0].strip("  "))
newseq = actualseq[:lenofseq]

但它似乎没有打印剪切序列。我能够将质量分数降低到我想要的位置。您可能已经注意到的另一件事是质量得分值没有适当间隔。我不确定为什么会发生这种情况。我使用" ".join() 作为质量分值，一旦它们被转换，这就是我得到的输出。

非常感谢您的建议！

Answer 1

dropwhile() 来自 itertools，当应用于反转数据时，会做你想做的事——它会过滤直到过滤器不再为真，然后完全退出，低质量分数将不再触发它：

from itertools import dropwhile

qualities = '1816 28 32 32 26 164 30 32 1816 1618 216 22 1616 216 218 20 28 1816 1816 24 1816 1620 20 24 28 28 28 216 26 28 1622 216 28 24 18 24 28 1622 30 222 24 18 218 216 26 26 218 28 1624 24 16 26 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2'.split()

sequence = 'ACCGAGCCGAAGGAGACCGCATTCACCCGGATGCCCTTCGAGGCCAGCGCCACCGCCATCGACCGCGTCATCTGCTCCACCGCGGCCCAGCTGATGGAACA'

QUALITY_LIMIT = 2

QUALITY_SCORE, BASE = 0, 1

filter = lambda pair: pair[QUALITY_SCORE] <= QUALITY_LIMIT

reversed_quality_sequence = zip(reversed([int(quality) for quality in qualities]), reversed(sequence))

filtered_sequence = "".join(reversed([pair[BASE] for pair in dropwhile(filter, reversed_quality_sequence)]))

print(filtered_sequence)

输出

ACCGAGCCGAAGGAGACCGCATTCACCCGGATGCCCTTCGAGGCCAGCGCCA

同样的技术也可以用于（更容易）在序列开始时清理低质量数据。

Answer 2

要仅从质量分数的末尾进行修剪，您可以在反向的quality_scores 列表上使用itertools.dropwhile() 来删除尾随的“2”项（感谢@cdlane 的这个想法）。然后将bases 切成修剪后的质量分数的长度。无需拉链：

from itertools import dropwhile

quality_scores = '1816 28 32 32 26 164 30 32 1816 1618 216 22 1616 216 218 20 28 1816 1816 24 1816 1620 20 24 28 28 28 216 26 28 1622 216 28 24 18 24 28 1622 30 222 24 18 218 216 26 26 218 28 1624 24 16 26 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2'.split()
bases = 'ACCGAGCCGAAGGAGACCGCATTCACCCGGATGCCCTTCGAGGCCAGCGCCACCGCCATCGACCGCGTCATCTGCTCCACCGCGGCCCAGCTGATGGAACA'

length = len(list(dropwhile(lambda x: x == '2', reversed(quality_scores))))
newseq = bases[:length]
print(newseq)

输出

ACCGAGCCGAAGGAGACCGCATTCACCCGGATGCCCTTCGAGGCCAGCGCCA