Pythonバイオ? Pythonバイオ/ツール?
2584 2019-06-09 (日) 09:07:15
やること：
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 （以下略）
count_tpm.tsvで得られたTMPデータで、batch_1, 2, 3, chemostat_1, 2, 3の６次元を PCAにより圧縮する。処理の基本は
![[Toho University]](image/CLogo_Kihon_process.jpg "[Toho University]")
![[YYLAB]](image/yamanouc2.png "[YYLAB]")
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')
出力は
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')
出力は
６つのサンプル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)
結果は
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" を指定しないと⽚側検定になってしまう 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がかなり小さく、棄却される。
カイ２乗独立性検定により、遺伝子Xと他の遺伝子との独立性を検定し、独立なものを探す。
# カイ２乗独立性の検定により、遺伝子ごとの発現量の有意差を見る
# 次の分割表に対するカイ２乗による独立性検定
# ｜ ｜batch ｜ chemostat
# ｜遺伝子１ ｜５０ ｜ １００
# ｜それ以外 ｜12,000,000 ｜13,000,000
# 帰無仮説：独立である を検定→p<<0.05なら棄却（独立でない）
#カイ⼆乗独⽴性の検定でとりあえず最初の10遺伝⼦を検定してみる
import numpy as np
#カイ⼆乗独⽴性の検定はカウントデータを使うので、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)
#カイ２乗独立性の検定では、各値のどれかが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）を⼀定の⽔準以下にする、FDRをp値の代わりに⽤いる⼿法での補正手法を試してみる。
# FDRの計算⽅法にはいくつか種類があるが、⼀番基本的なBenjamini Hochberg法では、各検定 で得られたp値を
# ⼩さい順に並べて、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 （以下略）
SRR453566-71-scatterplot.png 727件
[詳細]
SRR453566-71-PCA-Zscore.png 741件
[詳細]
SRR453566-71-PCA.png 829件
[詳細]
