分析现代序列文件

FSATA 格式是一个古老的格式了，现在生物信息学从业者更多地使用的是FASTQ 格式。
首先，我们先来获取我们此次要分析的 FASTQ 文件。

wget  ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/NA18489/sequence_read/SRR003265.filt.fastq.gz

现在我们来打开文件。

import gzip
from Bio import SeqIO
recs = SeqIO.parse(gzip.open('SRR003265.filt.fastq.gz','rt'),
'fastq')
rec = next(recs)
print(rec.id, rec.description, rec.seq)
print(rec.letter_annotations)

我们现在看看 reads 的分布。

from collections import defaultdict # defaultdict 是字典的一种，不需要对字典进行初始化
recs = SeqIO.parse(gzip.open('SRR003265.filt.fastq.gz'), 'fastq')
cnt = defaultdict(int) #此时字典的默认值为整数
for rec in recs: # 分析每一个 read
    for letter in rec.seq:
        cnt[letter] += 1 #统计每一个碱基的个数
tot = sum(cnt.values()) #碱基的总数
for letter, cnt_value in cnt.items():
    print('%s: %.2f %d' % (letter, 100. * cnt_value / tot, cnt_value)) #每个碱基的百分数和绝对数量

在这里面我们关心测序结果为 N 的碱基，N 表示测序没有测准。我们以 reads 的位置为自变量，N 的个数为因变量作图。

import seaborn as sns
import matplotlib.pyplot as plt
recs = SeqIO.parse(gzip.open('SRR003265.filt.fastq.gz'), 'fastq')
n_cnt = defaultdict(int)
for rec in recs:
    for i, letter in enumerate(rec.seq):
        pos = i + 1
        if letter == 'N':
            n_cnt[pos] += 1
seq_len = max(n_cnt.keys()) #用于设置横坐标长度
positions = range(1, seq_len + 1)
fig, ax = plt.subplots()
ax.plot(positions, [n_cnt[x] for x in positions])
ax.set_xlim(1, seq_len)

N 的数量随 read 位置的变化

我们进一步分析 Phred 分数的分布，Phred 分数衡量了测序的质量。

recs = SeqIO.parse(gzip.open('SRR003265.filt.fastq.gz'), 'fastq')
cnt_qual = defaultdict(int)
for rec in recs:
    for i, qual in enumerate(rec.letter_annotations['phred_quality']):
    if i < 25: #在上一个图中可以发现当 read 在25bp内时都测得很准，因此不需要得到其 Phred 分数
        continue
    cnt_qual[qual] += 1
tot = sum(cnt_qual.values())
for qual, cnt in cnt_qual.items():
    print('%d: %.2f %d' % (qual, 100. * cnt / tot, cnt))

对之前得 Phred 分数分布作图。

recs = SeqIO.parse(gzip.open('SRR003265.filt.fastq.gz'), 'fastq')
qual_pos = defaultdict(list)
for rec in recs:
    for i, qual in enumerate(rec.letter_annotations['phred_quality']):
        if i < 25 or qual == 40: # Phred 分数为40是这里面的最大值了，因此可以移除
            continue
        pos = i + 1
        qual_pos[pos].append(qual)
vps = []
poses = sorted(qual_pos) #通过键进行排序
for pos in poses:
    vps.append(qual_pos[pos]) #由字典转化为列表，便于作图
fig, ax = plt.subplots()
sns.boxplot(data=vps, ax=ax)
ax.set_xticklabels([str(x) for x in range(26, max(qual_pos.keys()) + 1)])

Phred 分数随 read 位置的变化

小结

在这一章中我们学习了：

• 如何读取 FASTQ 格式文件
• 如何分析 FASTQ 格式文件