01
PCA in Python
本文介绍如下内容:
1 构建可以用PCA的数据集
2 利用scikit-learn库的PCA函数做PCA工作
3 计算每个主成分的方差
4 利用matplotlib库做PCA图
5 通过loading scores分析变量的影响度
02
构建数据集
导入Python库
代码
import random as rd
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from sklearn.decomposition import PCA
from sklearn import preprocessing
import warnings
warnings.filterwarnings('ignore')
生成数据集
代码
genes = ['gene' + str(i) for i in range(1, 101)]
wt = ['wt' + str(i) for i in range(1, 6)]
ko = ['ko' + str(i) for i in range(1, 6)]
data = pd.DataFrame(columns=[*wt, *ko], index=genes)
for gene in data.index:
data.loc[gene, 'wt1':'wt5'] = np.random.poisson(lam=rd.randrange(10, 1000), size=5)
data.loc[gene, 'ko1':'ko5'] = np.random.poisson(lam=rd.randrange(10, 1000), size=5)
print(data.head())
结果
03
对数据集做PCA
利用sklearn库的PCA函数对数据集做PCA,进行PCA之前,对数据集做scale处理。
代码
scaled_data = preprocessing.scale(data.T)
pca = PCA()
pca.fit(scaled_data)
pca_data = pca.transform(scaled_data)
04
计算每个主成分的方差
计算每个主成分的方差和方差贡献度。
代码
# 计算每个主成分的方差
per_var = np.round(pca.explained_variance_, decimals=2)
print(per_var)
# 计算每个主成分方差的贡献度
per_var1 = np.round(pca.explained_variance_ratio_*100, decimals=1)
print(per_var1)
结果
05
PCA可视化
利用matplot库对PCA结果做可视化分析
代码
labels = ['PC' + str(x) for x in range(1, len(per_var) + 1)]
plt.bar(x=range(1, len(per_var1) + 1), height=per_var1, tick_label=labels)
plt.ylabel(u'方差贡献度', fontproperties='SimHei', fontsize=16)
plt.xlabel(u'主成分', fontproperties='SimHei', fontsize=16)
plt.title(u'Scree Plot', fontproperties='SimHei', fontsize=16)
plt.show()
结果
代码
pca_df = pd.DataFrame(pca_data, index=[*wt, *ko], columns=labels)
plt.scatter(pca_df.PC1, pca_df.PC2)
plt.title(u'PCA 图', fontproperties='SimHei', fontsize=16)
plt.xlabel(u'PC1 -{0}%'.format(per_var[0]))
plt.ylabel(u'PC2 -{0}%'.format(per_var[1]))
for sample in pca_df.index:
plt.annotate(sample, (pca_df.PC1.loc[sample], pca_df.PC2.loc[sample]))
plt.show()
结果
06
通过loading score分析变量影响度
通过第一主成分的loading score分析变量的影响度,即那些基因对第一主成分有着显著影响。
代码
loading_scores = pd.Series(pca.components_[0], index=genes)
sorted_loading_scores = loading_scores.abs().sort_values(ascending=False)
top_10_genes = sorted_loading_scores[0:10].index.values
print(top_10_genes)
print(loading_scores[top_10_genes])
结果
思考题:
1 Python做PCA和R做PCA有什么差异?