[[Pythonバイオ]] [[Pythonバイオ/ツール]]~
&counter();   &lastmod();~

*(遺伝研セミナー 3-1)発現解析結果の次元圧縮(3-1)と統計的検定(3-2) [#oeb2204f]

やること:
-前節で得た変異の値の表を、次元圧縮する
-検定


2-2で得た前提ファイル: diff_ex.tsv
 Geneid	batch	chemostat	log2fold	product
 gene_2989	0.5113545477349243	1609.2291342813207	11.619755389936394	Rgi2p
 gene_4740	3.714171096311252	5568.713898594607	10.550087794327904	Sip18p
 gene_4667	5.970597813523303	5108.597498100943	9.740835924216602	Spg4p
 gene_5965	8.609101140765873	5809.572900018754	9.398353606047285	Gre1p
 gene_4237	1.158010357867918	775.5514280425233	9.387429238134654	hypothetical protein
 (以下略)

***PCAによる次元圧縮 [#ge06e7b5]
count_tpm.tsvで得られたTMPデータで、batch_1, 2, 3, chemostat_1, 2, 3の6次元を
PCAにより圧縮する。処理の基本は
 df = pd.read_table('count_tpm.tsv', index_col=0)
 pca = sklearn.decomposition.PCA()
 coords = pca.fit_transform(df.transpose().values)
だが、得られたcoordsをプロットすること、入力データのうちオール0は除外すること
などの処理が必要。
 %matplotlib inline
 import pandas as pd
 import numpy as np
 import matplotlib.pyplot as plt
 import seaborn as sns
 import sklearn.decomposition
 df = pd.read_table('count_tpm.tsv', index_col=0)
 all_zero_index = df.index[df.sum(axis=1) == 0]
 df = df.drop(all_zero_index)
 
 sample_colors = {'batch_1':'b', # b : blue
     'batch_2':'b',
     'batch_3':'b',
     'chemostat_1':'g', # g: green
     'chemostat_2':'g',
     'chemostat_3':'g'}
 #colors = (df.columns).map(sample_colors)
 colors = [sample_colors[u] for u in df.columns]
 print(colors)
 # PCA 実行
 pca = sklearn.decomposition.PCA()
 coords = pca.fit_transform(df.transpose().values)
 
 def scatter_plot(coords, sample_labels, colors, xlabel=None, ylabel=None, title=''):
     fig = plt.figure()
     ax = fig.add_subplot(111)
     ax.scatter(coords[:, 0], coords[:, 1], color=colors)
     for i, sample_label in enumerate(sample_labels):
         ax.annotate(sample_label, xy=(coords[i, :2]), xytext=(10,10),
             textcoords='offset points', color=colors[i],
             arrowprops={'arrowstyle':'-', 'edgecolor':colors[i]})
     ax.set_xlabel(xlabel)
     ax.set_ylabel(ylabel)
     ax.set_title(title)
     plt.show()
     
 scatter_plot(coords, df.columns, colors, xlabel='PC1', ylabel='PC2')
出力は
#ref(SRR453566-71-PCA.png):

Zスコアによる標準化(データ群の当該数値から平均値を引いて標準偏差で割ったもの)
をしてPCA。

%matplotlib inline
 # Z-score 標準化
 import sklearn.preprocessing
 values = df.transpose().values
 scaler = sklearn.preprocessing.StandardScaler(with_mean=True, with_std=True)
 std_values = scaler.fit_transform(values)
 std_df = pd.DataFrame(std_values.T, index=df.index, columns=df.columns)
 print('\nRaw values')
 print('TPM Average :')
 print(df.mean(axis=1)[:10])
 print('\nTPM Standard deviation :')
 print(df.std(ddof=0, axis=1)[:10])
 print('\n\nStandardized values')
 print('TPM Average :')
 print(std_df.mean(axis=1)[:10])
 print('\nTPM Standard deviation :')
 print(std_df.std(ddof=0, axis=1)[:10])
 pca = sklearn.decomposition.PCA()
 coords = pca.fit_transform(std_df.transpose().values)
 scatter_plot(coords, std_df.columns, colors, xlabel='PC1', ylabel='PC2')

出力は
#ref(SRR453566-71-PCA-Zscore.png);

***統計的検定の活用〜相関・独立性検定 [#vba49c01]

6つのサンプルbatch_1、batcn_2、batch_3、chemstat_1, chemstat_2, chemstat_3の間の
相関を見る。まずは散布図を描く。
 %matplotlib inline
 import pandas as pd
 import seaborn as sns
 sns.set() #seaborn
 # 全組み合わせの散布図を表⽰
 tpm = pd.read_table("count_tpm.tsv") 
 sns.pairplot(tpm) 

結果は
#ref(SRR453566-71-scatterplot.png);

batch間、chemostat間は高い相関、batch-chemostat館はばらつきがあることがわかる。

相関係数を求めてみる。

 # Spearman相関係数を計算してみる
 from scipy import stats
 corr, p = stats.spearmanr(tpm.batch_1, tpm.batch_2) 
 print("Spearman correlation between batch_1 and batch_2:  " + str(corr)) 
 print("p value: " + str(p)) 
  
 corr, p = stats.spearmanr(tpm.batch_1, tpm.chemostat_1) 
 print("Spearman correlation between batch_1 and chemostat_1:  " + str(corr)) 
 print("p value: " + str(p)) 
結果は
 Spearman correlation between batch_1 and batch_2:  0.9903475599652781
   p value: 0.0
 Spearman correlation between batch_1 and chemostat_1:  0.8483821121349234
   p value: 0.0
相関係数で見ると、異種間の相関係数が明らかに小さいが、それでもそこそこ相関がある。

違いがあるかどうかを、Mann-Whitney U検定を用いて調べる。
2つの母集団が同じであるとするという仮定に対するp値。p<0.05(0.01)だと棄却→同じでない

 # 違いがあるかどうかを検定してみる
 # 母集団が正規分布しているとはいいにくいので、Mann-Whitney U検定を用いる
 b1c1 = stats.mannwhitneyu(tpm.batch_1, tpm.chemostat_1, alternative="two-sided")
 #                alternative="two-sided"  を指定しないと&#12122;側検定になってしまう
 b1c2 = stats.mannwhitneyu(tpm.batch_1, tpm.chemostat_2, alternative="two-sided") 
 b1c3 = stats.mannwhitneyu(tpm.batch_1, tpm.chemostat_3, alternative="two-sided") 
 b2c1 = stats.mannwhitneyu(tpm.batch_2, tpm.chemostat_1, alternative="two-sided") 
 b2c2 = stats.mannwhitneyu(tpm.batch_2, tpm.chemostat_2, alternative="two-sided") 
 b2c3 = stats.mannwhitneyu(tpm.batch_2, tpm.chemostat_3, alternative="two-sided") 
 b3c1 = stats.mannwhitneyu(tpm.batch_3, tpm.chemostat_1, alternative="two-sided") 
 b3c2 = stats.mannwhitneyu(tpm.batch_3, tpm.chemostat_2, alternative="two-sided") 
 b3c3 = stats.mannwhitneyu(tpm.batch_3, tpm.chemostat_3, alternative="two-sided") 
 print(b1c1.pvalue, b1c2.pvalue, b1c3.pvalue, "\n", b2c1.pvalue, b2c2.pvalue, 
 b2c3.pvalue, "\n", b3c1.pvalue, b3c2.pvalue, b3c3.pvalue) 
 
 #  ついでに、同種類batch_1とbatch_2の間ではどうか?
 print()
 b1b2 = stats.mannwhitneyu(tpm.batch_1, tpm.batch_2, alternative="two-sided")
 b1b3 = stats.mannwhitneyu(tpm.batch_1, tpm.batch_3, alternative="two-sided")
 b2b3 = stats.mannwhitneyu(tpm.batch_2, tpm.batch_3, alternative="two-sided")
 c1c2 = stats.mannwhitneyu(tpm.chemostat_1, tpm.chemostat_2, alternative="two-sided")
 c1c3 = stats.mannwhitneyu(tpm.chemostat_1, tpm.chemostat_3, alternative="two-sided")
 c2c3 = stats.mannwhitneyu(tpm.chemostat_2, tpm.chemostat_3, alternative="two-sided")
 print(b1b2.pvalue, b1b3.pvalue, b2b3.pvalue, '\n', \
    c1c2.pvalue, c1c3.pvalue, c2c3.pvalue)

結果は、異種間の場合
 7.443032122149763e-45 1.872248071308331e-35 5.5364679561011494e-30 
 6.146536587950258e-55 1.2997112467897929e-44 1.6801015887922584e-38 
 1.5468893669458462e-52 1.9165599784456836e-42 1.8061342891881474e-36
で、いずれもp<<0.01で棄却される(同じではない)

他方、同種間(b1-b2間)の場合
 0.08065503501024336 0.15188401219837053 0.7603592816804289 
 0.07125032415042876 0.006632138001689073 0.3395999834216805
で、いずれも棄却しないと言いたいが、chemostat_1と3の間はpがかなり小さく、棄却される。


カイ2乗独立性検定により、遺伝子Xと他の遺伝子との独立性を検定し、独立なものを探す。

 # カイ2乗独立性の検定により、遺伝子ごとの発現量の有意差を見る
 # 次の分割表に対するカイ2乗による独立性検定
 # |        |batch    | chemostat
 # |遺伝子1 |50    | 100
 # |それ以外 |12,000,000  |13,000,000
 # 帰無仮説:独立である を検定→p<<0.05なら棄却(独立でない)
 
 #カイ&#12038;乗独&#12148;性の検定でとりあえず最初の10遺伝&#12070;を検定してみる
 
 import numpy as np
 #カイ&#12038;乗独&#12148;性の検定はカウントデータを使うので、rawtagを使う
 rawtag = pd.read_table("count_raw.tsv")
 gene_ids = rawtag['Geneid']
 #batch 3 replicateのタグカウントの合計を計算
 sx = sum(rawtag.batch_1) + sum(rawtag.batch_2) + sum(rawtag.batch_3) 
 #chemostat 3 replicateのタグカウントの合計を計算
 sy = sum(rawtag.chemostat_1) + sum(rawtag.chemostat_2) + sum(rawtag.chemostat_3)
  
 ar = [] 
 for i in range(10): 
     #i番&#12140;の遺伝&#12070;のbatch 3 replicateでのタグカウントの合計を計算
     #i番目の遺伝子のbatch 3 replicateでのタグカウントの合計を計算
     x = rawtag.batch_1[i] + rawtag.batch_2[i] + rawtag.batch_3[i] 
     #i番&#12140;の遺伝&#12070;のchemostat 3 replicateでのタグカウントの合計を計算
     #i番目の遺伝子のchemostat 3 replicateでのタグカウントの合計を計算
     y = rawtag.chemostat_1[i] + rawtag.chemostat_2[i] + rawtag.chemostat_3[i] 
     
     #i番&#12140;の遺伝&#12070;以外の全遺伝&#12070;の、batch 3 replicateでのタグカウントの合計を計算
     #i番目の遺伝子以外の全遺伝子の、batch 3 replicateでのタグカウントの合計を計算
     sxr = sx - x
     #i番&#12140;の遺伝&#12070;以外の全遺伝&#12070;の、chemostat 3 replicateでのタグカウントの合計を計算
     #i番目の遺伝子以外の全遺伝子の、chemostat 3 replicateでのタグカウントの合計を計算
     syr = sy - y
 
     if ( rawtag.batch_1[i] == 0 or rawtag.batch_2[i] == 0 or rawtag.batch_3[i] == 0 or \
         rawtag.chemostat_1[i] == 0 or rawtag.chemostat_2[i] == 0 or \
         rawtag.chemostat_3[i] == 0): 
         ar.append(1) 
         #カイ&#12038;乗独&#12148;性の検定では、各値のどれかが0だと検定の精度が悪くなるため、
         #カイ2乗独立性の検定では、各値のどれかが0だと検定の精度が悪くなるため、
         #  そのようなデータの場合には検定をせずにp値を1とする
     else: 
         forchi = np.array([[x,y],[sxr,syr]]) 
         chi = stats.chi2_contingency(forchi) 
         ar.append(chi[1]) 
 print(ar) 
 
 print()
 # 上だと検定を繰り返している。それに対するp値の補正としてBonferroni補正を使う
 # Bonferroni補正は、0.05等の個々の検定のp値を、検定回数で割って得られた値を有意性の判定に使う
 import statsmodels.stats.multitest as multi
 sx = sum(rawtag.batch_1) + sum(rawtag.batch_2) + sum(rawtag.batch_3) 
 sy = sum(rawtag.chemostat_1) + sum(rawtag.chemostat_2) + sum(rawtag.chemostat_3) 
 ar = [] 
 for i in range(10): 
     x = rawtag.batch_1[i] + rawtag.batch_2[i] + rawtag.batch_3[i] 
     y = rawtag.chemostat_1[i] + rawtag.chemostat_2[i] + rawtag.chemostat_3[i] 
     sxr = sx - x 
     syr = sy - y 
     if ( rawtag.batch_1[i] == 0 or rawtag.batch_2[i] == 0 or rawtag.batch_3[i] == 0 or \
         rawtag.chemostat_1[i] == 0 or rawtag.chemostat_2[i] == 0 or rawtag.chemostat_3[i] == 0): 
         ar.append(1) 
     else: 
         forchi = np.array([[x,y],[sxr,syr]]) 
         chi = stats.chi2_contingency(forchi) 
         ar.append(chi[1]) 
 print(ar) 
 #α=0.05として、Bonferroni補正を&#12175;う
 #α=0.05として、Bonferroni補正を行う
 pval_corr = multi.multipletests(ar, alpha=0.05, method='Bonferroni') 
 print('Bonferroni:\n', pval_corr[1]) 
 
 # Bonferroni補正はかなり厳しい補正&#12101;法である。 
 # Bonferroni補正はかなり厳しい補正手法である。 
 # そこで、familywise error rateの代わりに、有意差ありと判定された結果の中に本当は差が無いものを含む確率
 # (False Discovery Rate)を&#12032;定の&#12116;準以下にする、FDRをp値の代わりに&#12132;いる&#12095;法での補正&#12101;法を試してみる。 
 # (False Discovery Rate)を&#12032;定の&#12116;準以下にする、FDRをp値の代わりに&#12132;いる&#12095;法での補正手法を試してみる。 
 # FDRの計算&#12101;法にはいくつか種類があるが、&#12032;番基本的なBenjamini Hochberg法では、各検定 で得られたp値を
 # &#12073;さい順に並べて、N/順位をp値にかけることで個別の検定におけるFDRを求める。
 # 多くの場合、FDRはp値と区別するためにq値と呼ばれる。 
 pval_corr = multi.multipletests(ar, alpha=0.05, method='fdr_bh') #Benjamini Hochberg法 
 print('FDR-BH', pval_corr[1]) 

結果は
 [1, 1, 1, 1, 0.019436546571414606, 1, 1, 1, 2.0214198154233473e-13, 0.0]
 
 [1, 1, 1, 1, 0.019436546571414606, 1, 1, 1, 2.0214198154233473e-13, 0.0]
 Bonferroni:
  [1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00
  1.94365466e-01 1.00000000e+00 1.00000000e+00 1.00000000e+00
  2.02141982e-12 0.00000000e+00]
 FDR-BH [1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00
  6.47884886e-02 1.00000000e+00 1.00000000e+00 1.00000000e+00
  1.01070991e-12 0.00000000e+00]


これを踏まえて、全遺伝子でトライ

 # 全遺伝で検定をし、FDRを求めてファイルに出力する
 # 独立であるという帰無仮定に対するp値。p<0.05(0.01)だと棄却→独立でない
 sx = sum(rawtag.batch_1) + sum(rawtag.batch_2) + sum(rawtag.batch_3) 
 sy = sum(rawtag.chemostat_1) + sum(rawtag.chemostat_2) + sum(rawtag.chemostat_3) 
 ai = len(rawtag.batch_1) 
 ar = [] 
 for i in range(ai): 
     x = rawtag.batch_1[i] + rawtag.batch_2[i] + rawtag.batch_3[i] 
     y = rawtag.chemostat_1[i] + rawtag.chemostat_2[i] + rawtag.chemostat_3[i] 
     sxr = sx - x 
     syr = sy - y 
     if ( rawtag.batch_1[i] <= 0 or rawtag.batch_2[i] <= 0 or rawtag.batch_3[i] <= 0 \
         or rawtag.chemostat_1[i] <= 0 or rawtag.chemostat_2[i] <= 0 or\
         rawtag.chemostat_3[i] <= 0): 
         ar.append(1) 
     else: 
         forchi = np.array([[x,y],[sxr,syr]]) 
         chi = stats.chi2_contingency(forchi) 
         ar.append(chi[1]) 
 pval_corr = multi.multipletests(ar, alpha=0.05, method='fdr_bh') 
 np.savetxt('BHtest.txt', pval_corr[1]) 
 gene_ids = list(rawtag['Geneid'])
 output = pd.DataFrame(list(zip(gene_ids, list(pval_corr[1]))))
 output.columns = ['Geneid', 'pvalue']
 output = output[output['pvalue']<0.05]
 print(output[:29])
 print('complete')

結果は
        Geneid         pvalue
 4   gene_0005   2.643073e-02
 8   gene_0009   4.575206e-13
 9   gene_0010   0.000000e+00
 10  gene_0011   0.000000e+00
 12  gene_0013   6.566874e-13
 13  gene_0014   1.237817e-08
 14  gene_0015   6.762139e-64
 15  gene_0016   6.777083e-18
 16  gene_0017   0.000000e+00
 17  gene_0018   1.301251e-34
 18  gene_0019   6.370921e-08
 19  gene_0020   8.823167e-71
 20  gene_0021   4.864046e-25
 (以下略)

トップ   編集 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS