歡迎關(guān)注公眾號：oddxix

主成分分析簡介

主成分分析 (PCA, principal component analysis)是一種數(shù)學(xué)降維方法, 利用正交變換(orthogonal transformation)把一系列可能線性相關(guān)的變量轉(zhuǎn)換為一組線性不相關(guān)的新變量，也稱為主成分，從而利用新變量在更小的維度下展示數(shù)據(jù)的特征。

主成分是原有變量的線性組合，其數(shù)目不多于原始變量。組合之后，相當(dāng)于我們獲得了一批新的觀測數(shù)據(jù)，這些數(shù)據(jù)的含義不同于原有數(shù)據(jù)，但包含了之前數(shù)據(jù)的大部分特征，并且有著較低的維度，便于進一步的分析。

在空間上，PCA可以理解為把原始數(shù)據(jù)投射到一個新的坐標(biāo)系統(tǒng)，第一主成分為第一坐標(biāo)軸，它的含義代表了原始數(shù)據(jù)中多個變量經(jīng)過某種變換得到的新變量的變化區(qū)間；第二成分為第二坐標(biāo)軸，代表了原始數(shù)據(jù)中多個變量經(jīng)過某種變換得到的第二個新變量的變化區(qū)間。這樣我們把利用原始數(shù)據(jù)解釋樣品的差異轉(zhuǎn)變?yōu)槔眯伦兞拷忉寴悠返牟町悺?/p>

這種投射方式會有很多，為了最大限度保留對原始數(shù)據(jù)的解釋，一般會用最大方差理論或最小損失理論，使得第一主成分有著最大的方差或變異數(shù) (就是說其能盡量多的解釋原始數(shù)據(jù)的差異)；隨后的每一個主成分都與前面的主成分正交，且有著僅次于前一主成分的最大方差 (正交簡單的理解就是兩個主成分空間夾角為90°，兩者之間無線性關(guān)聯(lián)，從而完成去冗余操作)。

主成分分析的意義

1.簡化運算

在問題研究中，為了全面系統(tǒng)地分析問題，我們通常會收集眾多的影響因素也就是眾多的變量。這樣會使得研究更豐富，通常也會帶來較多的冗余數(shù)據(jù)和復(fù)雜的計算量。

比如我們我們測序了100種樣品的基因表達譜借以通過分子表達水平的差異對這100種樣品進行分類。在這個問題中，研究的變量就是不同的基因。每個基因的表達都可以在一定程度上反應(yīng)樣品之間的差異，但某些基因之間卻有著調(diào)控、協(xié)同或拮抗的關(guān)系，表現(xiàn)為它們的表達值存在一些相關(guān)性，這就造成了統(tǒng)計數(shù)據(jù)所反映的信息存在一定程度的冗余。另外假如某些基因如持家基因在所有樣本中表達都一樣，它們對于解釋樣本的差異也沒有意義。這么多的變量在后續(xù)統(tǒng)計分析中會增大運算量和計算復(fù)雜度，應(yīng)用PCA就可以在盡量多的保持變量所包含的信息又能維持盡量少的變量數(shù)目，幫助簡化運算和結(jié)果解釋。

2.去除數(shù)據(jù)噪音

比如說我們在樣品的制備過程中，由于不完全一致的操作，導(dǎo)致樣品的狀態(tài)有細微的改變，從而造成一些持家基因也發(fā)生了相應(yīng)的變化，但變化幅度遠小于核心基因 (一般認為噪音的方差小于信息的方差）。而PCA在降維的過程中濾去了這些變化幅度較小的噪音變化，增大了數(shù)據(jù)的信噪比。

3.利用散點圖實現(xiàn)多維數(shù)據(jù)可視化

在上面的表達譜分析中，假如我們有1個基因，可以在線性層面對樣本進行分類；如果我們有2個基因，可以在一個平面對樣本進行分類；如果我們有3個基因，可以在一個立體空間對樣本進行分類；如果有更多的基因，比如說n個，那么每個樣品就是n維空間的一個點，則很難在圖形上展示樣品的分類關(guān)系。利用PCA分析，我們可以選取貢獻最大的2個或3個主成分作為數(shù)據(jù)代表用以可視化。這比直接選取三個表達變化最大的基因更能反映樣品之間的差異。（利用Pearson相關(guān)系數(shù)對樣品進行聚類在樣品數(shù)目比較少時是一個解決辦法）

4.發(fā)現(xiàn)隱性相關(guān)變量

我們在合并冗余原始變量得到主成分過程中，會發(fā)現(xiàn)某些原始變量對同一主成分有著相似的貢獻，也就是說這些變量之間存在著某種相關(guān)性，為相關(guān)變量。同時也可以獲得這些變量對主成分的貢獻程度。對基因表達數(shù)據(jù)可以理解為發(fā)現(xiàn)了存在協(xié)同或拮抗關(guān)系的基因。

問題描述

下表1是某些學(xué)生的語文、數(shù)學(xué)、物理、化學(xué)成績統(tǒng)計：

image

首先，假設(shè)這些科目成績不相關(guān)，也就是說某一科目考多少分與其他科目沒有關(guān)系。那么一眼就能看出來，數(shù)學(xué)、物理、化學(xué)這三門課的成績構(gòu)成了這組數(shù)據(jù)的主成分（很顯然，數(shù)學(xué)作為第一主成分，因為數(shù)學(xué)成績拉的最開）。為什么一眼能看出來？因為坐標(biāo)軸選對了！下面再看一組學(xué)生的數(shù)學(xué)、物理、化學(xué)、語文、歷史、英語成績統(tǒng)計，見表2，還能不能一眼看出來：

image.gif

數(shù)據(jù)太多了，以至于看起來有些凌亂！也就是說，無法直接看出這組數(shù)據(jù)的主成分，因為在坐標(biāo)系下這組數(shù)據(jù)分布的很散亂。究其原因，是因為無法撥開遮住肉眼的迷霧~如果把這些數(shù)據(jù)在相應(yīng)的空間中表示出來，也許你就能換一個觀察角度找出主成分。如下圖1所示：

image

但是，對于更高維的數(shù)據(jù)，能想象其分布嗎？就算能描述分布，如何精確地找到這些主成分的軸？如何衡量你提取的主成分到底占了整個數(shù)據(jù)的多少信息？所以，我們就要用到主成分分析的處理方法。

數(shù)據(jù)降維

為了說明什么是數(shù)據(jù)的主成分，先從數(shù)據(jù)降維說起。假設(shè)三維空間中有一系列點，這些點分布在一個過原點的斜面上，如果你用自然坐標(biāo)系x,y,z這三個軸來表示這組數(shù)據(jù)的話，需要使用三個維度，而事實上，這些點的分布僅僅是在一個二維的平面上，那么，問題出在哪里？如果你再仔細想想，能不能把x,y,z坐標(biāo)系旋轉(zhuǎn)一下，使數(shù)據(jù)所在平面與x,y平面重合？這就對了！如果把旋轉(zhuǎn)后的坐標(biāo)系記為x’,y’,z’，那么這組數(shù)據(jù)的表示只用x’和y’兩個維度表示即可！當(dāng)然了，如果想恢復(fù)原來的表示方式，那就得把這兩個坐標(biāo)之間的變換矩陣存下來。這樣就能把數(shù)據(jù)維度降下來了！但是，我們要看到這個過程的本質(zhì)，如果把這些數(shù)據(jù)按行或者按列排成一個矩陣，那么這個矩陣的秩就是2！這些數(shù)據(jù)之間是有相關(guān)性的，這些數(shù)據(jù)構(gòu)成的過原點的向量的最大線性無關(guān)組包含2個向量，這就是為什么一開始就假設(shè)平面過原點的原因！那么如果平面不過原點呢？這就是數(shù)據(jù)中心化的緣故！將坐標(biāo)原點平移到數(shù)據(jù)中心，這樣原本不相關(guān)的數(shù)據(jù)在這個新坐標(biāo)系中就有相關(guān)性了！有趣的是，三點一定共面，也就是說三維空間中任意三點中心化后都是線性相關(guān)的，一般來講n維空間中的n個點一定能在一個n-1維子空間中分析！

上一段文字中，認為把數(shù)據(jù)降維后并沒有丟棄任何東西，因為這些數(shù)據(jù)在平面以外的第三個維度的分量都為0?，F(xiàn)在，假設(shè)這些數(shù)據(jù)在z’軸有一個很小的抖動，那么我們?nèi)匀挥蒙鲜龅亩S表示這些數(shù)據(jù)，理由是我們可以認為這兩個軸的信息是數(shù)據(jù)的主成分，而這些信息對于我們的分析已經(jīng)足夠了，z’軸上的抖動很有可能是噪聲，也就是說本來這組數(shù)據(jù)是有相關(guān)性的，噪聲的引入，導(dǎo)致了數(shù)據(jù)不完全相關(guān)，但是，這些數(shù)據(jù)在z’軸上的分布與原點構(gòu)成的夾角非常小，也就是說在z’軸上有很大的相關(guān)性，綜合這些考慮，就可以認為數(shù)據(jù)在x’,y’ 軸上的投影構(gòu)成了數(shù)據(jù)的主成分！

PCA的思想是將n維特征映射到k維上（k<n），這k維是全新的正交特征。這k維特征稱為主成分，是重新構(gòu)造出來的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轉(zhuǎn)二進制

將merge.vcf轉(zhuǎn)二進制文件，測序結(jié)果call 的snp一般都是vcf格式，所以我們用到一個vcftools的軟件，該軟件只能在linux下使用。

##vcf轉(zhuǎn)格式，生成了tmp.map和tmp.ped兩個文件。
vcftools --vcf merged.vcf --plink --out tmp

##用plink轉(zhuǎn)成二進制文件，生成了分別以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 意思是只選出常染色體來運行（對應(yīng)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、我們這就已經(jīng)準備好了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)

說明：
我的樣品數(shù)據(jù)是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表示

附1：常用顏色col的代號

附2：圖形符號pch的代號

參考資料

• PCA主成分分析原理及分析實踐詳細介紹 | Public Library of Bioinformatics
• 利用GCAT做主成分分析（PCA） | Public Library of Bioinformatics
• 主成分分析（PCA）原理詳解 | Public Library of Bioinformatics
• PCA 文字化描述
• pca1

轉(zhuǎn)載請注明出處

歡迎關(guān)注公眾號：oddxix