バイオインフォマティクス Pythonによる実践レシピ (読書メモ)

バイオインフォマティクスに入門したいと思ったので、手にとって読んでみた本。

生物学の基本的な知識がかなり必要になるが、

調べながら手を動かしてみるとある種の現場の苦労がわかって楽しかった。

冒頭部分、手を動かしながら以下メモをとった。

書籍の Githubリンク

書籍にはGithubのリンクが書かれており、ここで Jupyter notebook形式でコードが提供されている。その殆どは軽く改変すると Google Colaboratory で動作する。

https://github.com/PacktPublishing/Bioinformatics-with-Python-Cookbook-Second-Edition

次世代シーケンス・データ解析

GenBankへのアクセス方法

National Center for Biotechnology information (NCBI)が運営するゲノムデータベースのうちの一つ、GenBankへのアクセス方法

NCBIから提供されるAPIをラップしたBiopythonを利用する。ピーク時間(米国東部標準時9-17時)を避ける、一度に100件以上のクエリを送らない、同時に4個以上のクエリを投入しない、実行時は必ずメールアドレスを設定する、などの注意が必要。

(書籍中では、ローカルにJupyter notebookを立ち上げて説明しているが、ここではColaboratory を使う)

Colaboratoryに Biopythonをインストールし、利用するモジュールをimportする

!pip install biopython

from Bio import Entrez, Medline, SeqIO
Entrez.email = "メールアドレス" 

handle = Entrez.einfo()
rec = Entrez.read(handle)
print(rec)

{'DbList': ['pubmed', 'protein', 'nuccore', 'ipg', 'nucleotide', 'structure', 'genome', 'annotinfo', 'assembly', 'bioproject', 'biosample', 'blastdbinfo', 'books', 'cdd', 'clinvar', 'gap', 'gapplus', 'grasp', 'dbvar', 'gene', 'gds', 'geoprofiles', 'homologene', 'medgen', 'mesh', 'ncbisearch', 'nlmcatalog', 'omim', 'orgtrack', 'pmc', 'popset', 'proteinclusters', 'pcassay', 'protfam', 'pccompound', 'pcsubstance', 'seqannot', 'snp', 'sra', 'taxonomy', 'biocollections', 'gtr']}

取得可能なデータベースの中から今回はnucleotideを使う。

NCBIの nucleotide データベースに登録された Plasmodium falciparum (熱帯熱マラリア原虫) の配列データから chloroquine resistance transporter (クロロキン耐性トランスポーター遺伝子参考に読んだ https://www.okayama-u.ac.jp/tp/release/release_id272.html)の遺伝子を取得する。

handle = Entrez.esearch(db="nucleotide", term='CRT[Gene Name] AND "Plasmodium falciparum"[Organism]')
rec_list = Entrez.read(handle)
if rec_list['RetMax'] < rec_list['Count']:
    handle = Entrez.esearch(db="nucleotide", term='CRT[Gene Name] AND "Plasmodium falciparum"[Organism]',
                            retmax=rec_list['Count'])
    rec_list = Entrez.read(handle)

1,2行目でデータベース名、生物名、遺伝子名を指定して検索している、3行目以降も同じことをしているが、デフォルトではretmaxが20なのでデータを全て取得できないため、rec_list['Count']に書き直して再実行している。次に取得された ID を引数に渡してレコード本体を取得。

id_list = rec_list['IdList']
hdl = Entrez.efetch(db='nucleotide', id=id_list, rettype='gb', retmax=rec_list['Count'])

書籍では id_list は 481 件となっているが、現在は 2022件であった。次に、結果を読み込み構文解析を実施していく、まずイテレータ形式で得られるデータをリスト形式に変換する(データが大きすぎるときはメモリを圧迫するのでこの方法は使わない)。少し時間がかかる。

recs = list(SeqIO.parse(hdl, 'gb'))

次に、リスト形式の中から１レコードを指定して結果を出力してみる

or rec in recs:
    if rec.name == 'KM288867':
        break
print(rec.name)
print(rec.description)

KM288867
Plasmodium falciparum clone PF3D7_0709000 chloroquine resistance transporter (CRT) gene, complete cds

レコードのdescriptionには説明が入っている。

例えばレコードの中からexon(RNAスプライシングの後に残る領域のこと https://ja.wikipedia.org/wiki/Pre-mRNA_%E3%82%B9%E3%83%97%E3%83%A9%E3%82%A4%E3%82%B7%E3%83%B3%E3%82%B0 ) の開始位置と終了位置を出力するためには次のように書く。

for feature in rec.features:

    if feature.type == 'exon':

        loc = feature.location

        print('Exon', loc.start, loc.end, loc.strand)

Exon 4513 4646 1
Exon 4799 4871 1
Exon 4994 5070 1
Exon 5166 5249 1

...

タンパク質への翻訳を見るには

for feature in rec.features:
    if feature.type == 'CDS':
        print(feature.qualifiers['translation'])
['MKFASKKNNQKNSSKNDERYRELDNLVQEGNGSRLGGGSCLGKCAHVFKLIFKEIKDNIFIYILSIIYLSVCVMNKIFAKRTLNKIGNYSFVTSETHNFICMIMFFIVYSLFGNKKGNSKERHRSFNLQFFAISMLDACSVILAFIGLTRTTGNIQSFVLQLSIPINMFFCFLILRYRYHLYNYLGAVIIVVTIALVEMKLSFETQEENSIIFNLVLISALIPVCFSNMTREIVFKKYKIDILRLNAMVSFFQLFTSCLILPVYTLPFLKQLHLPYNEIWTNIKNGFACLFLGRNTVVENCGLGMAKLCDDCDGAWKTFALFSFFNICDNLITSYIIDKFSTMTYTIVSCIQGPAIAIAYYFKFLAGDVVREPRLLDFVTLFGYLFGSIIYRVGNIILERKKMRNEENEDSEGELTNVDSIITQ']

基本的な配列解析方法

ヒトのLactase(LCT)遺伝子を解析する。

以下のようなコードで id を指定してレコードを取得する

hdl = Entrez.efetch(db='nucleotide', id=['NM_002299'], rettype='fasta')  # Lactase gene
seq = SeqIO.read(hdl, 'fasta')

ファイルへの保存は以下のようにする。FASTA形式で保存する。

w_seq = seq[11:5795]

w_hdl = open('example.fasta', 'w')
SeqIO.write([w_seq], w_hdl, 'fasta')
w_hdl.close()

ファイルからの読み込みは以下のようにする。

recs = SeqIO.parse('example.fasta', 'fasta')
for rec in recs:
    seq = rec.seq

(書籍にある

seq = Seq.Seq(str(seq), IUPAC.unambiguous_dna)

という記法は利用できなくなったとのこと 参考：https://biopython.org/wiki/Alphabet )

転写RNA配列を生成するには

rna = seq.transcribe()
rna
Seq('GAAAAUGGAGCUGUCUUGGCAUGUAGUCUUUAUUGCCCUGCUAAGUUUUUCAUG...AUU')

タンパク質配列データに変換するには

prot = seq.translate()
prot
Seq('ENGAVLACSLYCPAKFFMLGVRLGV**KFHFHRWSSNQ*LAAQPEWSPGRPEF*...VFI')

RNAからタンパク質配列に変換するときは、正しい遺伝暗号表を使っているか注意する。
例えばヒトのミトコンドリアは異なる遺伝暗号表を用いる。

FASTQフォーマットの利用
FASTQファイルは高速シーケンサーからの出力ファイル。サンプルに含まれる DNA, RNA などの断片の塩基配列
の出力フォーマット。
(次世代シーケンサーからの出力がどんなものかは、以下のチュートリアル参照
https://jp.illumina.com/science/technology/next-generation-sequencing/beginners.html )

ここでは1000人ゲノムプロジェクトのデータを利用する。1000人ゲノムプロジェクトには、
アフリカのヨルバ人(YRI)、北及び西ヨーロッパに祖先を持つアメリカユタ州住民(CEU)
東京の日本人(JPT)、北京の漢人(CHB) などが含まれている。

ファイルのダウンロード
!wget -nd ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/NA18489/sequence_read/SRR003265.filt.fastq.gz

FASTQファイルは配列データとクオリティスコアが合わせて格納されている。中身は

@SRR003265.31 3042NAAXX:3:1:1252:1819 length=51
GGGAAAAGAAAAACAAACAAACAAAAACAAAACACAGAAACAAAAAAACCA
+
IIIIIIIIIIIIIIIIIII?8IAD>I1IIAD@IIH7I95=@-@+7=.;588

1行目は @に続いて配列のID番号と配列の説明文字列。

2行目は FASTAファイルなどと同じ配列データ

3行目は +

4行目は2行目に書かれた配列データのクオリティスコア

読み込むには以下のようにする。以下は１レコードについて内容を出力する。

recs = SeqIO.parse(gzip.open('SRR003265.filt.fastq.gz', 'rt', encoding='utf-8'), 'fastq')
rec = next(recs)
print(rec.id, rec.description, rec.seq)
print(rec.letter_annotations)
SRR003265.31 SRR003265.31 3042NAAXX:3:1:1252:1819 length=51 GGGAAAAGAAAAACAAACAAACAAAAACAAAACACAGAAACAAAAAAACCA
{'phred_quality': [40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 40, 30, 23, 40, 32, 35, 29, 40, 16, 40, 40, 32, 35, 31, 40, 40, 39, 22, 40, 24, 20, 28, 31, 12, 31, 10, 22, 28, 13, 26, 20, 23, 23]}

全レコードについて塩基のカウント

recs = SeqIO.parse(gzip.open('SRR003265.filt.fastq.gz', 'rt', encoding='utf-8'), 'fastq')
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))
G: 20.68 5359334
A: 28.60 7411965
C: 21.00 5444053
T: 29.58 7666885
N: 0.14 37289

リード位置ごとの"N"の分布を調べる。N はシーケンサによる読み取りの失敗を表す。
今回使っているファイルはフィルタリング済みなので N の数は本来よりずっと少なくなっている。
recs = SeqIO.parse(gzip.open('SRR003265.filt.fastq.gz', 'rt', encoding='UTF-8'), '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(figsize=(16,9))
ax.plot(positions, [n_cnt[x] for x in positions])
ax.set_xlim(1, seq_len)
これを見ると 25 番目まで N が出現しないことがわかるが、これは 1000人ゲノムプロジェクトで 
25番目より前に N が出現しないようにフィルタリングしているためである。
その他の特徴として、N の出現頻度は一様ではなく、位置によって N の出現頻度が大きく変化することがわかる。
次に N ではなくクオリティスコアをプロットする。なお フィルタリングされているため 25番目より前は無視する。
recs = SeqIO.parse(gzip.open('SRR003265.filt.fastq.gz', 'rt', encoding='utf-8'), '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:
            continue
        pos = i + 1
        qual_pos[pos].append(qual)
vps = []
poses = list(qual_pos.keys())
poses.sort()
for pos in poses:
    vps.append(qual_pos[pos])
fig, ax = plt.subplots(figsize=(16,9))
sns.boxplot(data=vps, ax=ax)
ax.set_xticklabels([str(x) for x in range(26, max(qual_pos.keys()) + 1)])
pass



x軸はリードの開始位置からの距離。クオリティスコア Phred がだんだん下がっているのがわかる。

アラインメント・データ
シーケンサーからの配列データを取得したあとは通常 Burrows-Wheeler Aligner (bwa) などのアラインメントツール
を用いて配列をリファレンスゲノムにマッピングする。
アラインメントされたデータを扱うためには pysam などを用いる。
bam = pysam.AlignmentFile('NA18489.chrom20.ILLUMINA.bwa.YRI.exome.20121211.bam', 'rb')

BAMファイルの中でアラインメントによってうまくマップされた位置の分布を見てみる
counts = [0] * 76
for n, rec in enumerate(bam.fetch('20', 0, 10000000)):
    for i in range(rec.query_alignment_start, rec.query_alignment_end):
        counts[i] += 1
freqs = [x / (n + 1.) for x in counts]
fig, ax = plt.subplots(figsize=(16,9))
ax.plot(range(1, 77), freqs)

マップされる頻度は両端ではかなり低く、中央付近では落ち込んでいることがわかる。

VCFファイルGATKやSAMToolsの出力として、VCFファイルが出力される。VCFファイルの中には一塩基多型(SNP), 
挿入欠失(InDel), コピー数多型(CNV) などのゲノム変異をまとめたファイルである。

VCFファイルの一部だけを取り扱うためのコマンド tabix のインストール
# install tabix
!apt install tabix

VCFファイルをPython で取り扱うための pyvcf のインストール

!pip install pyvcf

VCFファイル内の変異の個数をカウントする
f = vcf.Reader(filename='genotypes.vcf.gz')
my_type = defaultdict(int)
num_alts = defaultdict(int)
for rec in f:
    my_type[rec.var_type, rec.var_subtype] += 1
    if rec.is_snp:
        num_alts[len(rec.ALT)] += 1
print(my_type)
print(num_alts)
defaultdict(<class 'int'>, {('snp', 'ts'): 10054, ('snp', 'tv'): 5917, ('sv', 'CNV'): 2, ('indel', 'del'): 273, ('snp', 'unknown'): 79, ('indel', 'ins'): 127, ('indel', 'unknown'): 13, ('sv', 'DEL'): 6, ('sv', 'SVA'): 1})
defaultdict(<class 'int'>, {1: 15971, 2: 79})

SNPのうち、転位(transition) が 2/3 転座(transversion) が 1/3 他に3個のアレルを持つ SNP が 79 個
ということがわかる。

ゲノムアクセシビリティとSNPデータのフィルタリング
ハマダラカ1000匹ゲノムプロジェクトのデータを利用
!tabix -fh ftp://ngs.sanger.ac.uk/production/ag1000g/phase1/preview/ag1000g.AC.phase1.AR1.vcf.gz 3L:1-200000 |bgzip -c > centro.vcf.gz
!tabix -fh ftp://ngs.sanger.ac.uk/production/ag1000g/phase1/preview/ag1000g.AC.phase1.AR1.vcf.gz 3L:21000000-21200000 |bgzip -c > standard.vcf.gz       
!tabix -p vcf centro.vcf.gz
!tabix -p vcf standard.vcf.gz
ハマダラカの3L染色体の中央部とセントロメア周辺のデータを取得




変異(両アレルSNP)の発生数をそれぞれのデータの各位置でカウントしたもの。
セントロメア周辺は染色体中央部に比べて変異が少ないことがわかる。

ゲノム解析

高品質レファレンスゲノムを用いた解析

熱帯熱マラリア原虫のゲノムデータのGC含有率(グアニンとシトシンの含有率)を調べる。

通常レファレンスゲノムはFASTAファイル形式で提供される。

!wget http://plasmodb.org/common/downloads/release-9.3/Pfalciparum3D7/fasta/data/PlasmoDB-9.3_Pfalciparum3D7_Genome.fasta
genome_name = 'PlasmoDB-9.3_Pfalciparum3D7_Genome.fasta'
recs = SeqIO.parse(genome_name, 'fasta')
chroms = {}
for rec in recs:
    print(rec.description)
Pf3D7_05_v3 | organism=Plasmodium_falciparum_3D7 | version=2012-02-01 | length=1343557 | SO=chromosome
Pf3D7_10_v3 | organism=Plasmodium_falciparum_3D7 | version=2012-02-01 | length=1687656 | SO=chromosome
Pf3D7_07_v3 | organism=Plasmodium_falciparum_3D7 | version=2012-02-01 | length=1445207 | SO=chromosome
Pf3D7_03_v3 | organism=Plasmodium_falciparum_3D7 | version=2012-02-01 | length=1067971 | SO=chromosome
Pf3D7_13_v3 | organism=Plasmodium_falciparum_3D7 | version=2012-02-01 | length=2925236 | SO=chromosome
Pf3D7_11_v3 | organism=Plasmodium_falciparum_3D7 | version=2012-02-01 | length=2038340 | SO=chromosome
Pf3D7_14_v3 | organism=Plasmodium_falciparum_3D7 | version=2012-02-01 | length=3291936 | SO=chromosome
Pf3D7_09_v3 | organism=Plasmodium_falciparum_3D7 | version=2012-02-01 | length=1541735 | SO=chromosome
Pf3D7_01_v3 | organism=Plasmodium_falciparum_3D7 | version=2012-02-01 | length=640851 | SO=chromosome
Pf3D7_12_v3 | organism=Plasmodium_falciparum_3D7 | version=2012-02-01 | length=2271494 | SO=chromosome
Pf3D7_08_v3 | organism=Plasmodium_falciparum_3D7 | version=2012-02-01 | length=1472805 | SO=chromosome
Pf3D7_06_v3 | organism=Plasmodium_falciparum_3D7 | version=2012-02-01 | length=1418242 | SO=chromosome
Pf3D7_04_v3 | organism=Plasmodium_falciparum_3D7 | version=2012-02-01 | length=1200490 | SO=chromosome
Pf3D7_02_v3 | organism=Plasmodium_falciparum_3D7 | version=2012-02-01 | length=947102 | SO=chromosome
M76611 | organism=Plasmodium_falciparum_3D7 | version=2012-02-01 | length=5967 | SO=mitochondrial_chromosome
PFC10_API_IRAB | organism=Plasmodium_falciparum_3D7 | version=2012-02-01 | length=34242 | SO=apicoplast_chromosome

レコードは染色体に対応しており、length は各染色体のサイズとなっている。

なお、SO=mitochondrial_chromosome はミトコンドリア染色体 SO=apicoplast_chromosome は アピコプラスト

各染色体の各部分における GC 含有率 を青の濃さで表す(最大値：22%, 最小値: 17.5% )

上側外れ値は赤、下側外れ値は黄。ウインドウサイズ 20kbp

aaatxt

このブログを検索