生存時間分析シリーズが立て続いておりましたが、今回は久しぶりにバイオインフォマティクス関連の記事になります。
以前の記事でも紹介したBiopythonを用いて、簡単なデータの可視化を行ってみます。
Biopythonのインストール
Biopythonのインストールはpipもしくはconda経由で行うことができます。
pip経由
$ pip install biopython
conda経由
$ conda install -c conda-forge biopython
FASTQファイルをとってくる
今回は1000人ゲノムプロジェクトのデータを用います。
1000人ゲノムプロジェクトのデータはFTPで公開されており、
Data portalから検索することができます。
FASTQファイルのフォーマットって何？って方は下記の記事を参考にして下さい。
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スコアに関してもこちらの記事で説明しているので参考にしてみて下さい。
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については下記で説明しています。
Phread Scoreの分布
Ｎ数 / 塩基数
