Pythonでバイオインフォマティクス - FASTQファイルを扱う

f:id:kimoppy126:20201018231834p:plain

生存時間分析シリーズが立て続いておりましたが、今回は久しぶりにバイオインフォマティクス関連の記事になります。
以前の記事でも紹介したBiopythonを用いて、簡単なデータの可視化を行ってみます。

www.kimoton.com

Biopythonのインストール

Biopythonのインストールはpipもしくはconda経由で行うことができます。

pip経由

$ pip install biopython

conda経由

$ conda install -c conda-forge biopython

FASTQファイルをとってくる

今回は1000人ゲノムプロジェクトのデータを用います。
1000人ゲノムプロジェクトのデータはFTPで公開されており、 Data portalから検索することができます。

FASTQファイルのフォーマットって何?って方は下記の記事を参考にして下さい。

www.kimoton.com

wgetコマンドでダウンロードできます。

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

配列情報を取得

Bio.SeqIOパッケージを使うと、ファイルの読み書き含め配列情報を直感的に扱うことができます。

from Bio import SeqIO

# SeqIO.parseはFASTAやFASTQファイルのパースに使用される
recs = SeqIO.parse(gzip.open('SRR003265.filt.fastq.gz', 'rt', encoding='utf-8'), 'fastq')
# SeqIOの返り値はrecordのジェネレータ
rec = next(recs)

print(type(rec))
# .id内に1行目@以降のIDが格納されている
print(f"ID: {rec.id}")
# .description内に1行目全体が格納されている
print(f"Description: {rec.description}")
# .seqに2行目の配列が格納されている
print(f"Seq: {rec.seq}")
# .letter_annotationsに4行目のPhread qualityが格納されている
print(f"Phread Quality: {rec.letter_annotations['phred_quality']}")
<class 'Bio.SeqRecord.SeqRecord'>
ID: SRR003265.9754138
Description: SRR003265.9754138 3042NAAXX:3:100:1896:1223 length=51
Seq: ATTTATAAAAAAAATAAAAAACAAAAAACCGACACAAAATAAAAACCATNA
Phread Quality: [14, 8, 3, 2, 8, 3, 4, 2, 9, 5, 10, 2, 7, 2, 3, 5, 3, 4, 3, 3, 4, 10, 5, 4, 6, 5, 3, 5, 5, 6, 2, 9, 3, 11, 4, 3, 3, 5, 3, 2, 3, 2, 3, 3, 3, 3, 4, 9, 2, 0, 2]

ファイル中の各塩基割合を算出

配列の統計情報として基礎中の基礎が各塩基がどれだけ含まれているかのカウントです。
以下では、A, T, G, Cそれぞれの塩基がファイル中の全リード中にどれだけ含まれているのかをカウントします。

from collections import Counter, defaultdict
from Bio import SeqIO

# SeqIO.parseはFASTAやFASTQファイルのパースに使用される
recs = SeqIO.parse(gzip.open('SRR003265.filt.fastq.gz', 'rt', encoding='utf-8'), 'fastq')
# defaultdict型を使ってデフォルト値を整数型に指定
cnt = defaultdict(int)
# 各配列中の出現文字数をカウント
for rec in recs:
    for letter in rec.seq:
        cnt[letter] += 1
# 各配列の合計出現文字数を算出
tot = sum(cnt.values())
for letter, cnt in cnt.items():
    print('%s: %.2f %d' % (letter, 100. * cnt / tot, cnt))
A: 68.63 35
T: 13.73 7
C: 13.73 7
G: 1.96 1
N: 1.96 1

「N」数をカウント

FASTQファイル中の「N」は、シーケンサーが読み取りに失敗した塩基を示しています。
ここで用いているFASTQファイルは、フィルタリング後のファイル(*.filt.fastq.gz)を使用しているため、「N」の頻度はかなり低くなっています。

また、1000人ゲノムプロジェクトでは、25番目より前に「N」が出現しないようにフィルタリングされているため、25番目までの「N」数は0となっています。

from collections import defaultdict
from Bio import SeqIO
import matplotlib.pyplot as plt

# SeqIO.parseはFASTAやFASTQファイルのパースに使用される(3度目の登場)
recs = SeqIO.parse(gzip.open('SRR003265.filt.fastq.gz', 'rt', encoding='UTF-8'), 'fastq')
# defaultdict型を使ってデフォルト値を整数型に指定
n_cnt = defaultdict(int)
# 各配列に渡って、「N」数をカウント
for rec in recs:
    for i, letter in enumerate(rec.seq):
        pos = i + 1
        if letter == 'N':
            n_cnt[pos] += 1
# プロットの最大x軸に使用するため、最大の配列開始位置からの距離を取得
seq_len = max(n_cnt.keys())
# プロットのx軸を作成
positions = range(1, seq_len + 1)

# x軸を配列開始位置からの距離、y軸を「N」の出現頻度とした折れ線グラフを作成
plt.figure()
fig, ax = plt.subplots(figsize=(16,9))
ax.set_title("'N' count distribution")
ax.plot(positions, [n_cnt[x] for x in positions])
ax.set_xlim(1, seq_len)
plt.show()

Phread Scoreの分布を可視化

Phredスコアに関してもこちらの記事で説明しているので参考にしてみて下さい。

www.kimoton.com

Seabornを使うと、かっちょいいFigureを作成することができます。

from collections import defaultdict
from Bio import SeqIO
import seaborn as sns
import matplotlib.pyplot as plt

# SeqIO.parseはFASTAやFASTQファイルのパースに使用される(4度目の登場)
recs = SeqIO.parse(gzip.open('SRR003265.filt.fastq.gz', 'rt', encoding='utf-8'), 'fastq')
# defaultdict型を使ってデフォルト値をリスト型に指定
qual_pos = defaultdict(list)
# 各配列に渡って、ポジションごとのPhred Qualityをqual_posに格納
for rec in recs:
    for i, qual in enumerate(rec.letter_annotations['phred_quality']):
        # スコアが40未満のみ可視化の対象とする
        if qual == 40:
            continue
        pos = i + 1
        qual_pos[pos].append(qual)
vps = []
poses = list(qual_pos.keys())
# ポジションを昇順にソート
poses.sort()
# Phred Scoreをポジション順にプロット用のリスト(vps)に格納
for pos in poses:
    vps.append(qual_pos[pos])  
# x軸を配列開始位置からの距離、y軸をPhred Scoreとしたボックスプロットを作成
plt.figure()
fig, ax = plt.subplots(figsize=(16,9))
sns.boxplot(data=vps, ax=ax)
ax.set_xticklabels([str(x + 1) for x in range(max(qual_pos.keys()))])
plt.show()

FastQCを使って同様の可視化

FastQCを使うと、上記のようにわざわざスクリプトを作成しなくても配列のクオリティを確認するために必要なプロットを出力してくれます。上記の情報以外にも配列ごとの平均クオリティスコアやGC含量、配列長の分布などが含まれます。

FastQCについては下記で説明しています。

www.kimoton.com

Phread Scoreの分布

N数 / 塩基数

参考