Pythonバイオ/ツール/RNAseq④-2 Python版発現解析（続きの続き）の変更点

追加された行はこの色です。
削除された行はこの色です。
Pythonバイオ/ツール/RNAseq④-2 Python版発現解析（続きの続き）へ行く。
Pythonバイオ/ツール/RNAseq④-2 Python版発現解析（続きの続き）の差分を削除
[[Pythonバイオ]]　[[Pythonバイオ/ツール]]~
&counter();　　　&lastmod();~

*TPM以降の解析 [#h2dccce6]
Cufflink・Cuffdiff・あたりでTPM (FPRM) を出してからの処理。

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

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

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

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

**全体の発現量データのチェック [#v06dfb17]
***発現量の分布の表示 [#lf15dc44]
# 東北大　二階堂愛　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=['-', '--', ':', '-', '-.', ':'])  # 線のパターンは５種類しかない
 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()

結果は
&ref(density-distribution.png);

同様のデータを、箱ひげ図で描く。
***同様のデータを、箱ひげ図で描く。 [#r15dda03]
 %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()

結果は
&ref(box-chart.png);

***サンプル間の散布図を描く [#zc9f847b]
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()

結果は、
&ref(scatter-plot-large.png);

&ref(scatter-plot-average.png);

***ヒートマップ [#u69ad453]

項目をクラスタ化した結果を使って寄せる。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()

結果は
&ref(heatmap-dendrogram.png);

&ref(heatmap.png);

***JSD (Jense-Shannon divergence) 距離によるgene間のふるまいの違いの算定 [#j0b52a21]
たとえば、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()

結果は

&ref(JSD-similar.png);
Pythonバイオ/ツール/RNAseq④-2 Python版発現解析（続きの続き） の変更点

Pythonバイオ/ツール/RNAseq④-2 Python版発現解析（続きの続き）の変更点
