Pythonバイオ? Pythonバイオ/ツール?
1921 2019-07-10 (水) 16:05:41
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
![[Toho University]](image/CLogo_Kihon_process.jpg "[Toho University]")
![[YYLAB]](image/yamanouc2.png "[YYLAB]")
%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=['-', '--', ':', '-', '-.', ':']) # 線のパターンは５種類しかない
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()
結果は
%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()
結果は
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()
結果は、
項目をクラスタ化した結果を使って寄せる。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を取る前に１を足す
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()
結果は
たとえば、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の行は除かなければならない（確率なので全て足して１に正規化する）
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()
結果は
scatter-plot-large.png 574件
[詳細]
scatter-plot-average.png 542件
[詳細]
JSD-similar.png 534件
[詳細]
heatmap-dendrogram.png 559件
[詳細]
heatmap.png 567件
[詳細]
density-distribution.png 575件
[詳細]
box-chart.png 496件
[詳細]
