欢迎关注公众号：oddxix

主成分分析简介

主成分分析 (PCA, principal component analysis)是一种数学降维方法, 利用正交变换(orthogonal transformation)把一系列可能线性相关的变量转换为一组线性不相关的新变量，也称为主成分，从而利用新变量在更小的维度下展示数据的特征。

主成分是原有变量的线性组合，其数目不多于原始变量。组合之后，相当于我们获得了一批新的观测数据，这些数据的含义不同于原有数据，但包含了之前数据的大部分特征，并且有着较低的维度，便于进一步的分析。

在空间上，PCA可以理解为把原始数据投射到一个新的坐标系统，第一主成分为第一坐标轴，它的含义代表了原始数据中多个变量经过某种变换得到的新变量的变化区间；第二成分为第二坐标轴，代表了原始数据中多个变量经过某种变换得到的第二个新变量的变化区间。这样我们把利用原始数据解释样品的差异转变为利用新变量解释样品的差异。

这种投射方式会有很多，为了最大限度保留对原始数据的解释，一般会用最大方差理论或最小损失理论，使得第一主成分有着最大的方差或变异数 (就是说其能尽量多的解释原始数据的差异)；随后的每一个主成分都与前面的主成分正交，且有着仅次于前一主成分的最大方差 (正交简单的理解就是两个主成分空间夹角为90°，两者之间无线性关联，从而完成去冗余操作)。

主成分分析的意义

1.简化运算

在问题研究中，为了全面系统地分析问题，我们通常会收集众多的影响因素也就是众多的变量。这样会使得研究更丰富，通常也会带来较多的冗余数据和复杂的计算量。

比如我们我们测序了100种样品的基因表达谱借以通过分子表达水平的差异对这100种样品进行分类。在这个问题中，研究的变量就是不同的基因。每个基因的表达都可以在一定程度上反应样品之间的差异，但某些基因之间却有着调控、协同或拮抗的关系，表现为它们的表达值存在一些相关性，这就造成了统计数据所反映的信息存在一定程度的冗余。另外假如某些基因如持家基因在所有样本中表达都一样，它们对于解释样本的差异也没有意义。这么多的变量在后续统计分析中会增大运算量和计算复杂度，应用PCA就可以在尽量多的保持变量所包含的信息又能维持尽量少的变量数目，帮助简化运算和结果解释。

2.去除数据噪音

比如说我们在样品的制备过程中，由于不完全一致的操作，导致样品的状态有细微的改变，从而造成一些持家基因也发生了相应的变化，但变化幅度远小于核心基因 (一般认为噪音的方差小于信息的方差）。而PCA在降维的过程中滤去了这些变化幅度较小的噪音变化，增大了数据的信噪比。

3.利用散点图实现多维数据可视化

在上面的表达谱分析中，假如我们有1个基因，可以在线性层面对样本进行分类；如果我们有2个基因，可以在一个平面对样本进行分类；如果我们有3个基因，可以在一个立体空间对样本进行分类；如果有更多的基因，比如说n个，那么每个样品就是n维空间的一个点，则很难在图形上展示样品的分类关系。利用PCA分析，我们可以选取贡献最大的2个或3个主成分作为数据代表用以可视化。这比直接选取三个表达变化最大的基因更能反映样品之间的差异。（利用Pearson相关系数对样品进行聚类在样品数目比较少时是一个解决办法）

4.发现隐性相关变量

我们在合并冗余原始变量得到主成分过程中，会发现某些原始变量对同一主成分有着相似的贡献，也就是说这些变量之间存在着某种相关性，为相关变量。同时也可以获得这些变量对主成分的贡献程度。对基因表达数据可以理解为发现了存在协同或拮抗关系的基因。

问题描述

下表1是某些学生的语文、数学、物理、化学成绩统计：

image

首先，假设这些科目成绩不相关，也就是说某一科目考多少分与其他科目没有关系。那么一眼就能看出来，数学、物理、化学这三门课的成绩构成了这组数据的主成分（很显然，数学作为第一主成分，因为数学成绩拉的最开）。为什么一眼能看出来？因为坐标轴选对了！下面再看一组学生的数学、物理、化学、语文、历史、英语成绩统计，见表2，还能不能一眼看出来：

image.gif

数据太多了，以至于看起来有些凌乱！也就是说，无法直接看出这组数据的主成分，因为在坐标系下这组数据分布的很散乱。究其原因，是因为无法拨开遮住肉眼的迷雾~如果把这些数据在相应的空间中表示出来，也许你就能换一个观察角度找出主成分。如下图1所示：

image

但是，对于更高维的数据，能想象其分布吗？就算能描述分布，如何精确地找到这些主成分的轴？如何衡量你提取的主成分到底占了整个数据的多少信息？所以，我们就要用到主成分分析的处理方法。

数据降维

为了说明什么是数据的主成分，先从数据降维说起。假设三维空间中有一系列点，这些点分布在一个过原点的斜面上，如果你用自然坐标系x,y,z这三个轴来表示这组数据的话，需要使用三个维度，而事实上，这些点的分布仅仅是在一个二维的平面上，那么，问题出在哪里？如果你再仔细想想，能不能把x,y,z坐标系旋转一下，使数据所在平面与x,y平面重合？这就对了！如果把旋转后的坐标系记为x’,y’,z’，那么这组数据的表示只用x’和y’两个维度表示即可！当然了，如果想恢复原来的表示方式，那就得把这两个坐标之间的变换矩阵存下来。这样就能把数据维度降下来了！但是，我们要看到这个过程的本质，如果把这些数据按行或者按列排成一个矩阵，那么这个矩阵的秩就是2！这些数据之间是有相关性的，这些数据构成的过原点的向量的最大线性无关组包含2个向量，这就是为什么一开始就假设平面过原点的原因！那么如果平面不过原点呢？这就是数据中心化的缘故！将坐标原点平移到数据中心，这样原本不相关的数据在这个新坐标系中就有相关性了！有趣的是，三点一定共面，也就是说三维空间中任意三点中心化后都是线性相关的，一般来讲n维空间中的n个点一定能在一个n-1维子空间中分析！

上一段文字中，认为把数据降维后并没有丢弃任何东西，因为这些数据在平面以外的第三个维度的分量都为0。现在，假设这些数据在z’轴有一个很小的抖动，那么我们仍然用上述的二维表示这些数据，理由是我们可以认为这两个轴的信息是数据的主成分，而这些信息对于我们的分析已经足够了，z’轴上的抖动很有可能是噪声，也就是说本来这组数据是有相关性的，噪声的引入，导致了数据不完全相关，但是，这些数据在z’轴上的分布与原点构成的夹角非常小，也就是说在z’轴上有很大的相关性，综合这些考虑，就可以认为数据在x’,y’ 轴上的投影构成了数据的主成分！

PCA的思想是将n维特征映射到k维上（k<n），这k维是全新的正交特征。这k维特征称为主成分，是重新构造出来的k维特征，而不是简单地从n维特征中去除其余n-k维特征。

利用GCAT做主成分分析（PCA）实践

1.对样品文件进行处理

原始文件：vcf格式

#对raw.vcf文件处理
cp ../haplotype_GATK/*raw.vcf* ./   

#生成GZ格式
for i in *vcf ; do bgzip $i;done

#生成tbi文件格式
for i in *vcf.gz ; do tabix -p vcf $i;done
#合并
ls *vcf.gz | xargs vcf-merge > merged.vcf

image

2.merge.vcf转二进制

将merge.vcf转二进制文件，测序结果call 的snp一般都是vcf格式，所以我们用到一个vcftools的软件，该软件只能在linux下使用。

##vcf转格式，生成了tmp.map和tmp.ped两个文件。
vcftools --vcf merged.vcf --plink --out tmp

##用plink转成二进制文件，生成了分别以bed bim fam为后缀的tmp文件。
plink --noweb --file tmp --make-bed --out tmp_bfile

ps：plink里面—file 选项后跟文件名，不用加后缀~

image

3.生成matrix

gcta跟eigenstrat软件包做pca的效果是一样的，而且gcta比eigenstrat容易使用的多了，所以单纯做pca的话用gcta就好了，做gcta分两步。

gcta64 --bfile tmp_bfile --make-grm --autosome --out tmp_grm

说明：

1）tmp_bfile是你上一步plink生成的二进制文件（不包括后缀名）

2）tmp_grm是你指定输出的文件名

3）—autosome 意思是只选出常染色体来运行（对应1-22号染色体）

gcta64 --grm tmp_grm --pca 3 --out tmp_pca

说明：
1）tmp_grm是你上一步生成的文件名，不包含.grm.gz这个后缀

2）tmp_pca是输出文件

3）这样得到两个文件一个是tmp_pca.eigenval另一个是tmp_pca.eigenvec。在后者行首加入一行：1 2 pc1 pc2 pc3（分隔符为空格），并保存。
3、我们这就已经准备好了R作图用的matrix了，接下来用R作图。下载好R，然后两步走：

tmp_pca.eigenvec

4.R作图

先把matrix读入到R中

plot(data$pc1,data$pc2,pch=c(rep(c(1,2,3),8),rep(c(4,5),6)),col=c(rep("#FF0000FF",3),rep("#FF6D00FF",3),rep("#FFDB00FF",3),rep("#B6FF00FF",3),rep("#49FF00FF",3),rep("#00FF24FF",3),rep("#00FF92FF",3),rep("#00FFFFFF",3),rep("#0092FFFF",2),rep("#0024FFFF",2),rep("#4900FFFF",2),rep("#B600FFFF",2),rep("#FF00DBFF",2),rep("#FF006DFF",2)),lwd=2.4,cex=2.5,main="PCA",xlab="pc1",ylab="pc2")

##增加图示信息
legend("bottomright",c("chb","jpt"),data$X1,pch=c(rep(c(1,2,3),8),rep(c(4,5),6)),col=c(rep("#FF0000FF",3),rep("#FF6D00FF",3),rep("#FFDB00FF",3),rep("#B6FF00FF",3),rep("#49FF00FF",3),rep("#00FF24FF",3),rep("#00FF92FF",3),rep("#00FFFFFF",3),rep("#0092FFFF",2),rep("#0024FFFF",2),rep("#4900FFFF",2),rep("#B600FFFF",2),rep("#FF00DBFF",2),rep("#FF006DFF",2)),lwd=0.2,cex=0.3)

说明：
我的样品数据是8个样本需要3维分析，6个样品进行2维分析

1）pch=c(rep(c(1,2,3),8),rep(c(4,5),6))意思是把8个样本用pch=1，pch=2，pch=3的图案来表示，6个样本用pch=4，pch=5表示