PCA

PCA实现

简单地说，PCA的过程就是求协方差矩阵特征向量的过程。

下面是教程中很普遍的一个例子：

import numpy as np
from sklearn.datasets import load_iris
from sklearn.decomposition import PCA
import matplotlib.pyplot as plt

data = load_iris()     #加载鸢尾花数据
X = data.data    #X为花的各个特征
Y = data.target    #Y为花的分类标签

newData=X-np.mean(X, axis=0)    #数据去中心化
covMat=np.cov(newData,rowvar=0)     #求协方差矩阵
featValue,featVec=np.linalg.eig(covMat)    #求特征值和特征向量
index=np.argsort(featValue) #特征值排序
n_index=index[-1:-3:-1]    #选择特征值最大的2个数的标签
n_featVec=featVec[:, n_index]    #选择2个特征值对应特征向量
lowData=np.dot(newData,n_featVec)    #计算最后结果
#highData=np.dot(lowData,n_featVec.T)+meanval



plt.figure(figsize=(8,4))
plt.subplot(121)
plt.title("my_PCA")
plt.scatter(lowData[:, 0], lowData[:, 1], c = Y)


群体结构的PCA

方法一：

plink --bfile myfile --pca 10 --out myfile_pca

方法二：

gcta64 --bfile ./myfile  --make-grm --autosome --autosome-num 24  --out ./myfile_grm
gcta64 --grm ./myfile_grm_grm --pca 10 --out ./all.genotype_pca

方法三：

flashpca_x86-64 --bfile myfile

方法四：

library(vcfR)
library(adegenet)
library(adegraphics)
library(pegas)
library(lattice)
library(ape)

vcf <- read.vcfR("test.vcf.gz")
aa.genlight <- vcfR2genlight(vcf, n.cores=1)
pca_result <- glPca(aa.genlight, nf=3, n.cores=1) 
pca_result$eig[1]/sum(pca_result$eig)    # pc1 解释的比例
pca_result$eig[2]/sum(pca_result$eig)    # pc2 解释的比例
pca_result$eig[3]/sum(pca_result$eig)    # pc3 解释的比例
pca_score <- as.data.frame(pca_result$scores) #作图数据

确定显著PCA数量

twstats -t twtable -i pca.eigenval -o eigenvaltw.out

在结果文件查看 p-value列，如果值小于0.05则表示选择这个主成分。