Pythonバイオ? Pythonバイオ/ツール?
88   2019-06-09 (日) 09:07:15

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

やること:

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による次元圧縮

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')

出力は

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')

出力は

SRR453566-71-PCA-Zscore.png

統計的検定の活用〜相関・独立性検定

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) 

結果は

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番目の遺伝子のbatch 3 replicateでのタグカウントの合計を計算
    x = rawtag.batch_1[i] + rawtag.batch_2[i] + rawtag.batch_3[i] 
    #i番目の遺伝子のchemostat 3 replicateでのタグカウントの合計を計算
    y = rawtag.chemostat_1[i] + rawtag.chemostat_2[i] + rawtag.chemostat_3[i] 
    
    #i番目の遺伝子以外の全遺伝子の、batch 3 replicateでのタグカウントの合計を計算
    sxr = sx - x
    #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) 
        #カイ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補正を行う
pval_corr = multi.multipletests(ar, alpha=0.05, method='Bonferroni') 
print('Bonferroni:\n', pval_corr[1]) 

# Bonferroni補正はかなり厳しい補正手法である。 
# そこで、familywise error rateの代わりに、有意差ありと判定された結果の中に本当は差が無いものを含む確率
# (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
(以下略)

添付ファイル: fileSRR453566-71-scatterplot.png 23件 [詳細] fileSRR453566-71-PCA-Zscore.png 25件 [詳細] fileSRR453566-71-PCA.png 24件 [詳細]

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