今回の発表に用いたjupyter notebookは練習用データと一緒にhttps://github.com/yyoshiaki/drbono.chap4 に置いてあります。¶
今回のノートブックを開くには、以下の手順。
- anaconda(python3)をインストール。
$ git clone https://github.com/yyoshiaki/drbono.chap4.gitとターミナルでうち、今回のリポジトリをcloneしてくる。gitがない場合は適宜インストール。$ jupyter notebookを立ち上げて、20180414_drbonobon_chap4.ipynbというファイルを実行する。
自己紹介¶
- 安水 良明
- 大阪大学医学部医学科6年
- 免疫学フロンティアセンター 実験免疫学（坂口研）
- Python使ってます。
- 最近Bioinformatics contest 2018(ITMO大学)で入賞して喜んでいます。
4.1配列データ解析¶
- 配列アラインメントとツール
- ペアワイズアラインメント
- ドットプロット
- 大域的アラインメント（Global Alignment）
- 局所的アラインメント（Local Alignment）
- 多重配列アラインメントと系統樹
ペアワイズアラインメント, ドットプロット¶
- ２本の配列を比較すること。
- 同じ生物で比べることもあるし、違う生物間で比べることもある。
- 配列解析の出発点。
- それを可視化するのがドットプロット。
だいたいのことはEMBOSSで実行できます。
! head seq/L07770.1.fasta
! head seq/U23808.2.fasta
大域的アラインメント（Global Alignment）と 局所的アラインメント（Local Alignment）¶
大域的アラインメント（Global Alignment）¶
大域的アラインメント とは配列中の全塩基やアミノ酸を並べるようにする方法。 現在Needleman-Wunch法と呼ばれる方法が広く使われている。
局所的アラインメント（Local Alignment）¶
大域的アライメントの方法を改良して、部分的な類似性が見つけられるようにしたのが局所的アライメント。局所的アライメントの方法として最もよく用いられているのが Smith-Waterman法。query配列をDBから検索してくるという作業はこの局所的アラインメントのおかげで可能になった。
EMBOSSでの実行例は大石さんのスライドが参考になります。挙動をpythonライブラリのscikit-bioで見てみましょう。
from skbio.sequence import DNA
from skbio.alignment import *
s1 = DNA("AGGTGCGTGCCTAAGGTATGCAAGAGATGG")
s2 = DNA("ACGTGCCTAGGTACGCAAG")
# 教育用の関数なので、動作は遅い。下記URLを参照。
# http://scikit-bio.org/docs/0.2.0/alignment.html?highlight=alignment#module-skbio.alignment
r = global_pairwise_align_nucleotide(s1, s2)
# r = local_pairwise_align_nucleotide(s1, s2)
print(r)
局所アラインメントのバリエーション¶
Smith-Waterman法 先程の例で使用。FASTAパッケージに含まれているssearchが並列計算を用いたDB配列検索によく用いられていた。
FASTA法 かつてよく使われた。現在はBLASTがよく使われるが、FANTOMプロジェクトでcDNA配列解析をしていた際は、塩基配列をアミノ酸に翻訳して配列比較を行うプログラムのfastxとfastyが有用だったそう。
BLAST 現在おなじみの配列類似性検索ツール。
BLAT 非常に高速に、しかもエクソンとイントロンの境界も考慮して、ゲノムに対するアラインメントがなされる。ほぼ一致する領域を探すことに特化している。BLASTはもともと遠く離れた類縁関係を検出することを目的としていた。営利目的の使用は有償。
BLASTとBLATを試してみます。
データベース作成¶
! head seq/chr22.hg19.fa
! makeblastdb -in seq/chr22.hg19.fa -dbtype nucl
!ls seq
検索¶
今回はDiGeorge症候群に関連するDGCR6のmRNAから配列を取ってきました。
! head seq/dgcr6.fa
! blastn -query seq/dgcr6.fa -db seq/chr22.hg19.fa
扱いやすい形で保存することも出来る¶
! blastn -query seq/dgcr6.fa -db seq/chr22.hg19.fa -outfmt 7 -out blast.out.txt
! head blast.out.txt
BLAST使い方まとめ¶
- makeblastdbでDBを作る
- blast?で検索
| query\DB | 塩基配列 | アミノ酸配列 |
|---|---|---|
| 塩基配列 | blastn, tblastx | blastx |
| アミノ酸配列 | tblastn | blastp |
NCBIのwebサイトやbiopythonを使ってpythonプログラムに組み込むことも可能。
!blat seq/chr22.hg19.fa seq/dgcr6.fa output.psl
!head output.psl
ラウス肉腫ウイルスが宿主のニワトリから持ち出したと考えられる配列をニワトリゲノムに対して検索する、といったことも可能である。
多重配列アラインメントと系統樹¶
多重配列アラインメント(multiple sequence alignment)は3本以上の配列を並べて、出来る限りギャップをいれないようにして、似たアミノ酸（もしくは塩基）を並べる手法である。アラインメントした結果から分子系統樹を推定し、その結果からアラインメントを更に改良するという作業が行われる。clustal-omegaやmafftが使われる。どちらともwebインターフェースとしても公開されている。
多重配列アラインメントの結果はJalview等のソフトウェアで可視化出来る。詳しくは書籍もしくはtogo tvを参照。
4.1配列データ解析¶
- マッピング（Suffix Array） 重要!!!
- BWAとBowtie
- GGRNとGGGenome
BWAとBowtie¶
BWAとBowtieはNGSから出てきたデータをreference genome配列にマッピングすることに特化したソフトウェア。検索する対象DB（reference genome配列）が予めBurrows Wheeler Transform(BWT)で処理されており、そのため高速に、短い配列に対しても処理を行うことが出来る。
マッピングの流れは基本的に以下。
- インデックスの作成（referenceに対して一回で良い。）
- マッピング
インストール¶
condaがおすすめ。個人的にmacでマッピングはしたくない。condaだとsudo権限のないlinuxでも簡単にインストール出来る。
$ conda install -c bioconda bowtie
$ conda install -c bioconda bwa
alignerについて諸々¶
bwaは伝統的にゲノム解析で使われる印象。bowtieは少し古く、bowtie2、更にはhisat、hisat2へと開発が進んでいる。さらに、STARやLASTなど様々な選択肢がある。トランスクリプトーム解析だと最近ではkallistoやsalmonなど、pseudo-alignmentと呼ばれる高速でトランスクリプトームの定量まで行えるソフトウェアが人気。どのワークフローが良いかは日夜議論の対象となっている。
出典 Mohammad, S. et al. Gaining comprehensive biological insight into the transcriptome by performing a broad-spectrum RNA-seq analysis. doi:10.1038/s41467-017-00050-4
GGRNAとGGGenome¶
上述のBWAやBowtieと同様の検索をGoogle検索のように行えるDBCLSのサービス。CRISPR/CasのガイドRNA設計ツールにも応用されている。
! curl http://gggenome.dbcls.jp/hg38/1/AGGTVANNTGACCT.gff
4.1配列データ解析¶
- アセンブル
最近の次世代シーケンサーが長いreadが読めるようになったとはいえ、一本で完全長のDNAなどを読むのはまだ難しい。そこで、次世代シーケンサーによって得られたreadをコンピューター上でつなぎ合わせて再構築する。それがAssemble。
transcriptome assembleならTrinityがよく使われる¶
Trinityはguid referenceなしのアセンブルに用いられる。あまり解析の進んでいない生物にくわえ、ヒトにおけるfusion geneやvirusゲノムの組み込みなどを解析する際にも使える。癌ゲノムに焦点を当てたTrinity Cancer Transcriptome Analysis Toolkit( https://github.com/NCIP/Trinity_CTAT/wiki )なども公開されている。guid referenceが信頼できる場合はStringtieを使うこともできる。
Pythonとは¶
文法を極力単純化してコードの可読性を高め、読みやすく、また書きやすくしてプログラマの作業性とコードの信頼性を高めることを重視してデザインされた、汎用の高水準言語である。 wikipedia 「Python」より
フリーかつオープンソースの言語。バイオインフォマティクスだけではなく幅広いコミュニティをもち、web構築やGUI構築、科学演算、機械学習、データベース、ひいてはドローンの操作まで出来る万能言語。
Rとの比較¶
様々なところでRとpythonの比較はされている。
統計やバイオインフォマティクス特化ならRで十分。DESeq2など高機能ライブラリもRで用意されている。最新の機械学習やドローンなど、他のこともしたい欲があるならpythonがおすすめ。juliaもじわじわ人気。
4.2 数値データ解析¶
- 階層クラスタリング
- PCA
まとめてpythonでやってみましょう。 データはImmuno-NavigatorからTregのマイクロアレイデータをダウンロードした。
Vandenbon, A. et al. Immuno-Navigator, a batch-corrected coexpression database, reveals cell type-specific gene networks in the immune system. Proc. Natl. Acad. Sci. U. S. A. 113, E2393-402 (2016).
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
% matplotlib inline
# immuno-NavigatorからダウンロードしたTregの発現ファイルを読み込み。
df_exp = pd.read_csv('data/expression_3_human.txt', sep='\t')
# Affymetrixのprobe IDの詳細データをダウンロード。これを使って発現テーブルのProbe IDを置き換える。
df_annot = pd.read_csv('data/HG-U133_Plus_2-na36-annot-csv/HG-U133_Plus_2.na36.annot.csv', comment='#')
df_annot = df_annot[['Probe Set ID', 'Gene Symbol']]
df_annot.columns = ['ProbeID', 'GeneSymbol']
df_exp['ProbeID'] = df_exp.index
df_exp = pd.merge(df_exp, df_annot, on='ProbeID', how='inner')
df_exp = df_exp.drop(['ProbeID'], axis=1)
# 一つの遺伝子に複数のProbeが割り当てられている事があるので、平均を取った。
df_exp = df_exp.groupby('GeneSymbol').mean()
df_exp = df_exp[df_exp.index != '---']
df_exp
# 変動係数を出してみる
(df_exp.std(axis=1) / df_exp.mean(axis=1)).sort_values(ascending=False).head(10)
sns.clustermap(df_exp.loc[['FOXP3', 'CTLA4', 'IL2RA', 'STAT5A', 'STAT5B', 'CD4', 'RORA', 'IL3', 'IL4']],
z_score=0, cmap='bwr', method='ward')
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
# from sklearn.manifold import TSNE
# PCAの前に正規化したほうが良いです。
# http://scikit-learn.org/stable/auto_examples/preprocessing/plot_scaling_importance.html
# 正規化
X = StandardScaler().fit_transform(df_exp.T)
# PCA
pca = PCA(n_components=50)
pca.fit(X)
result = pca.transform(X)
# plot
plt.scatter(result[:,0], result[:,1])
# plt.savefig('img/pca.png')
# 固有ベクトル(Loadings)を分析
df_comp = pd.DataFrame(pca.components_)
df_comp.columns = df_exp.index
df_comp.loc[0].abs().sort_values(ascending=False).head(20)
df_comp.loc[1].abs().sort_values(ascending=False).head(20)
# 寄与率
plt.bar(range(50),pca.explained_variance_ratio_)
plt.title('explained variance ratio')
plt.xlabel('PC')
# plotlyで更にわかりやすく
import plotly.offline as offline
import plotly.graph_objs as go
offline.init_notebook_mode()
# Create a trace
trace = go.Scatter(
x = result[:,0],
y = result[:,1],
mode = 'markers',
text = df_exp.columns,
marker = dict(
color = df_exp.loc['FOXP3'],
size='16',
colorscale='Viridis',
showscale=True,
colorbar = dict(title = 'FOXP3')
)
)
data = [trace]
layout=dict(title='PCA')
fig=dict(data=[trace], layout=layout)
# はみ出るので次のスライドいきます。
# offline.iplot(fig, filename='PCA')
offline.iplot(fig, filename='PCA')
使ったライブラリまとめ¶
- scikit-learn : Google社製機械学習ライブラリ。PCAのために利用。死角なし。
- pandas : 同じくデータサイエンス3本柱の一つ。データフレームオブジェクトと豊富な機能が利用可能。
- matplotlib : 鉄板可視化ライブラリ。これを知らずにpythonを知っているとは言えない。
- seaborn : クラスターマップに使用。matplotlibをリッチにしてくれる。
plotly : インタラクティブな可視化が出来るライブラリ。これも熱い。
scikit-bio : alignmentの例で使用。
Pythonの勉強を始めたいと思った方へ¶
- udemy 【世界で5万人が受講】実践 Python データサイエンス : 17.5hのビデオ講座。私はこれで始めました。有料。
- rosalind : bioinformaticsの練習サイト。少し慣れてきたら。
- その他書籍、オンライン学習等色々あります。
情報収集手段（Pythonに限らないが。）¶
- qiita : プログラマのための技術情報共有サービス。
- kaggle : データサイエンスコンペティションのプラットフォーム。kernelなど、可視化やデータ処理の参考に。
- twitter,blog : バイオインフォマティクス界隈の人や情報界隈の人は割とやっている人が多い。「bounohu blog」「六本木で働くデータサイエンティストのブログ」「驚異のアニヲタ社会復帰への道」「macでインフォマティクス」などなど他分野にわたる。twitterで有名論文公式や有名研究者などをフォローしたり、RSSでまとめてブログをよんだり。
- google : 言わずもがな。なんでも調べる。エラーが出たらエラーをコピペ。