Pythonバイオ? Pythonバイオ/ツール?
22   2019-07-10 (水) 16:05:41

TPM以降の解析

Cufflink・Cuffdiff・あたりでTPM (FPRM) を出してからの処理。

全体にデータはそこそこまともか?
発現パターンの似た遺伝子は?

RとBioconductorでNGS解析://catway.jp/lecture/tohoku-u/NGS-R-Bioconductor-2nd.html (二階堂愛)を追う形で

「Cuffdiffの結果をcummeRbundで読込む」 ⇒ FRPMデータ

こちらは、既にここまでできている。

全体の発現量データのチェック

発現量の分布の表示

# 東北大 二階堂愛 http://catway.jp/lecture/tohoku-u/NGS-R-Bioconductor-2nd.html

%matplotlib inline
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
sns.set()
sns.set_style('whitegrid')

# count_tpm.tsvの読込み

df = pd.read_csv('count_tpm.tsv', sep='\t', index_col=0)
all_zero_index = df.index[df.sum(axis=1) == 0]
df = df.drop(all_zero_index)
dfl = np.log10(df)
dfl = dfl.replace([np.inf, -np.inf], -1)  # 値がinf/-infのデータは-1に置換
print(dfl.head())
#dfl.plot.kde(style=['-', '--', ':', '-', '-.', ':'])  # 線のパターンは5種類しかない
sns.kdeplot(dfl['batch_1'])  # 別々に書いて重ねることにした
sns.kdeplot(dfl['batch_2'])
sns.kdeplot(dfl['batch_3'])
sns.kdeplot(dfl['chemostat_1'])
sns.kdeplot(dfl['chemostat_2'])
sns.kdeplot(dfl['chemostat_3'])
plt.xlim(-1, 5)
plt.xlabel('$log_{10}(TPM)$')
plt.ylabel('density')
plt.show()

結果は density-distribution.png

同様のデータを、箱ひげ図で描く。

%matplotlib inline
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set()
sns.set(style='whitegrid', palette='Set3', font=['IPAexGothic'])
# count_tpm.tsvの読込み
df = pd.read_csv('count_tpm.tsv', sep='\t', index_col=0)
all_zero_index = df.index[df.sum(axis=1) == 0]
df = df.drop(all_zero_index)
dfl = np.log10(df)
dfl = dfl.replace([np.inf, -np.inf], -1)
print(dfl.head())
#dfl.boxplot()
#dfl.plot.box()
sns.boxplot(data=dfl)
#plt.xlim(-0.3, 5)
plt.ylabel('$log_{10}(TPM)$')
plt.show()

結果は box-chart.png

サンプル間の散布図を描く

seabornの得意技pairplotを使うと何も考えなくてよい。

%matplotlib inline
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
sns.set()
#sns.set(style='whitegrid', palette='Set3', font=['IPAexGothic'], markers='.')
sns.set(style='whitegrid', font=['IPAexGothic'])
# count_tpm.tsvの読込み
df = pd.read_csv('count_tpm.tsv', sep='\t', index_col=0)
all_zero_index = df.index[df.sum(axis=1) == 0]
df = df.drop(all_zero_index)
dfl = np.log10(df)
dfl = dfl.replace([np.inf, -np.inf], -1)
print(dfl.head())
#dfl.boxplot()
#dfl.plot.box()
sns.pairplot(dfl, kind='reg', height=3)  # size=0.1とすると図間が詰まる
#plt.xlim(-0.3, 5)
#plt.xlabel('$log_{10}(TPM)$')
plt.show()
dfl['batch'] = dfl[['batch_1', 'batch_2', 'batch_3']].mean(axis=1) # 横方向の平均
dfl['chemostat'] = dfl[['chemostat_1', 'chemostat_2', 'chemostat_3']].mean(axis=1)
print(dfl.head())
sns.relplot(kind='scatter', x='batch', y='chemostat', data=dfl) # size=0.1とすると点が小さくなる
plt.show()

結果は、 scatter-plot-large.png

scatter-plot-average.png

ヒートマップ

項目をクラスタ化した結果を使って寄せる。dendrogramの戻り値を使う。

%matplotlib inline
# heatmap クラスター化によって寄せる
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.cluster.hierarchy import linkage, dendrogram, leaves_list
from scipy.spatial.distance import pdist
sns.set()
#sns.set(style='whitegrid', palette='Set3', font=['IPAexGothic'])
sns.set(style='whitegrid', font=['IPAexGothic'], palette='Spectral')
# count_tpm.tsvの読込み
df = pd.read_csv('count_tpm.tsv', sep='\t', index_col=0)
#all_zero_index = df.index[df.sum(axis=1) == 0]
#df = df.drop(all_zero_index)
#df = df.replace(0, 0.00001)
print(df.head())
print(df[-5:])
# 絶対値が大きい100genes(+に大きい50と-に大きい50)を取る
#dfr = pd.concat([df[:200], df[-200:]])
 dfr = pd.concat([df[:50], df[-50:]])
 print(dfr.head())
 print(dfr[-5:])
#df = df.applymap(lambda x: x+1)
dfr = dfr + 1  # logを取る前に1を足す
print('gene_0008\n', dfr.loc['gene_0008'])
# log10を取る
dfl = np.log10(dfr)
#dfl = dfl.replace([np.inf, -np.inf], -5)

dfl_t = dfl.T
print('dfl\n', dfl.head())
#linkage_result_gene = linkage(dfl, method='average', metric='correlation')
linkage_result_gene = linkage(dfl)
plt.figure(figsize=(16,20))
res = dendrogram(linkage_result_gene, labels=dfl.index, leaf_rotation=90, 
leaf_font_size=10, get_leaves=True)
leavesy = leaves_list(linkage_result_gene)
print('leavesy', leavesy)
lvs = pd.DataFrame([res['ivl'], res['leaves']], columns=res['ivl']).T
lvs.columns = ['gene_id2', 'leaves']
print(lvs)

dfl2 = pd.merge(lvs, dfl, left_on='gene_id2', right_index=True)
print('dfl2\n', dfl2.head())
print('dfl2\n', dfl2[-5:])

#lvsの並び順のままの方がいい
print(dfl2.head())
plt.figure(figsize=(10, 50))
sns.set_palette('hot')
#sns.heatmap(dfl2[['batch_1','batch_2','batch_3', 'chemostat_1','chemostat_2','chemostat_3']], vmax=1, vmin=-1, center=0)
sns.heatmap(dfl2[['batch_1','batch_2','batch_3', 'chemostat_1','chemostat_2','chemostat_3']], vmax=3, vmin=0, center=1.0, cmap='YlOrBr')
plt.show()

結果は heatmap-dendrogram.png

heatmap.png

JSD (Jense-Shannon divergence) 距離によるgene間のふるまいの違いの算定

たとえば、gene_0013に似たふるまいのgeneを探す。最も近い10個を表示。

%matplotlib inline
# JSD Jensen-Shannon divergenceによる、違いの算定
# ある遺伝子と似た発言を示す遺伝子を探す
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.cluster.hierarchy import linkage, dendrogram, leaves_list
from scipy.spatial import distance
import itertools
sns.set()
sns.set(style='whitegrid', font=['IPAexGothic'], palette='Spectral')
# count_tpm.tsvの読込み
df = pd.read_csv('count_tpm.tsv', sep='\t', index_col=0)

df = df + 1
# log10を取る
dfl = np.log10(df)
dfl_t = dfl.T
# すべての行の組合せで、距離を計算する
# オール0の行は除かなければならない(確率なので全て足して1に正規化する)
dfl_copy = dfl.replace(0.0, np.nan).dropna(how='all', axis=0)
dfl = dfl.loc[dfl_copy.index]
# print('dfl\n', dfl.head()); print()

# dfl = dfl[:10]
# targetと他のgeneとの対の距離を計算する
target = 'gene_0013'
dfl['d'] = dfl.apply(lambda x: distance.jensenshannon(x, dfl.loc[target]), axis=1)
#print(dfl)
dfx = dfl.sort_values('d', ascending=True)[:10]
print('dfx\n', dfx)

dfg = dfx.drop('d', axis=1).T
dfg.plot()
plt.ylabel('$log_{10}(TPM-1)$')
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.show()

結果は

JSD-similar.png


添付ファイル: filescatter-plot-large.png 4件 [詳細] filescatter-plot-average.png 2件 [詳細] fileJSD-similar.png 3件 [詳細] fileheatmap-dendrogram.png 2件 [詳細] fileheatmap.png 2件 [詳細] filedensity-distribution.png 2件 [詳細] filebox-chart.png 2件 [詳細]

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2019-07-10 (水) 16:05:41 (12d)