近日，做差異分析的時(shí)候，想著看一下樣本本身的特征是以什么分類的，除了計(jì)算樣本之間的距離，還用到的PCA（主成分分析）。在DESeq2包中專門由一個(gè)PCA分析的函數(shù)，即plotPCA，里面的參數(shù)也比較簡(jiǎn)單。
plotPCA參數(shù)
object：對(duì)象
在進(jìn)行PCA之前是要進(jìn)行歸一化的，因?yàn)樵诮稻S映射的過程中， 存在映射誤差， 所有在對(duì)高維特征降維之前， 需要做特征歸一化（feature normalization）, 這個(gè)歸一化操作包括： （1） feature scaling (讓所有的特征擁有相似的尺度， 要不然一個(gè)特征特別小， 一個(gè)特征特別大會(huì)影響降維的效果) （2） zero mean normalization (零均值歸一化)。、
但是在DESeq2包中實(shí)際上已經(jīng)有了歸一化的方法，rlog和vst，具體是怎么做到歸一化的，這個(gè)原理不了解。
在使用的需要根據(jù)樣本量的多少來選擇方法。樣本量少于30的話，選擇rlog，多于30的話，建議選擇vst。親測(cè)！樣本量多于30的時(shí)候，用了rlog，出現(xiàn)以下提示：
rlog() may take a few minutes with 30 or more samples,
vst() is a much faster transformation
所以官方就會(huì)建議樣本量大于30的時(shí)候用vst，而且vst速度是飛起的狀態(tài)，而rlog是老牛拉慢車。
intgroup：分組
官方解釋：interesting groups: a character vector of names in colData(x) to use for grouping。也就是在構(gòu)建dds的時(shí)候的分組，實(shí)際上這個(gè)分組最后影響的是PCA圖中的顏色，但是并不影響PCA圖中各個(gè)樣本的位置。換句話說，PCA降維的是時(shí)候并不會(huì)考慮這個(gè)分組。這個(gè)也是我搞了好多天才弄清楚的，所以原理還是要搞懂的，不能只會(huì)跑流程。
ntop：用于主成分分析的時(shí)候基因數(shù)
官方解釋：number of top genes to use for principal components, selected by highest row variance。也就是說PCA分析的時(shí)候，是根據(jù)基因表達(dá)的counts值來的。如果給出了ntop，那么意思就是只選出表達(dá)量高的這幾個(gè)基因去計(jì)算。沒想通作者為什么設(shè)置這個(gè)，我覺得應(yīng)該是先篩選出差異的基因，然后再設(shè)置這個(gè)數(shù)會(huì)比較好，這樣就是看篩選出的基因是不是符合預(yù)期。反正感覺ntop目前來說用的不是很多。
returnData：返回PC1和PC2的dataframe
官方解釋：should the function only return the data.frame of PC1 and PC2 with intgroup covariates for custom plotting (default is FALSE)。注意只返回PC1和PC2，其他的成分不返回，而且還返回分組情況和樣本名。是這種格式：
PC1 PC2 group condition name
X11_T 23.7381223 -3.9537594 II II X11_T
X77_T -1.7273395 29.0222667 II II X77_T
X14_T -14.2359870 1.8616294 II II X14_T
X76_T -0.2251326 27.9834964 II II X76_T
X72_T -9.9218391 19.7866404 II II X72_T
X20_T 12.1552460 -19.1133494 II II X20_T
X71_T -45.6232065 -1.7467169 II II X71_T
這個(gè)參數(shù)的作用就是當(dāng)你的PCA圖需要加樣本名的時(shí)候，也就是你希望在PCA圖上知道哪個(gè)點(diǎn)是哪個(gè)樣本，你就需要導(dǎo)出這個(gè)。
其實(shí)我的初衷也就是想看在PCA圖上哪個(gè)點(diǎn)是哪個(gè)樣本，因此我當(dāng)然要設(shè)置returnData=T。
#構(gòu)建dds
dds <- DESeqDataSetFromMatrix(countData=indata, colData = state, design= ~ condition )
#歸一化，因?yàn)闃颖玖砍^了30，因此用vst
vsd <- vst(dds, blind = FALSE)
#返回樣本名和分組
pcaData <- plotPCA(vsd, intgroup=c("condition"),returnData = T)
#這里按照condition排序了，原因見下。
pcaData <- pcaData[order(pcaData$condition,decreasing=F),]
#知道每一個(gè)組有多少樣本
table(pcaData$condition)
# II III IV
# 20 11 10
#根據(jù)上面的結(jié)果，設(shè)置每一個(gè)組的數(shù)量，方便加顏色；這一步是把每一個(gè)樣本根據(jù)分組情況畫出來，效果見圖1
plot(pcaData[,1:2],pch = 19,col= c(rep("red",20),rep("green",11),rep("blue",10)),
cex=1)
#加上樣本名字，效果見圖2
text(pcaData[,1],pcaData[,2],row.names(pcaData),cex=1, font = 1)
#加上圖例，效果見圖3
legend(-70,43,inset = .02,pt.cex= 1.5,title = "Grade",legend = c("II", "III","IV"),
col = c( "red","green","blue"),pch = 19, cex=0.75,bty="n")
下面說一下我畫這個(gè)圖的艱辛的過程。
首先，給PCA上每一個(gè)點(diǎn)加樣本名這個(gè)過程其實(shí)是分步驟進(jìn)行的。也就是我三個(gè)圖的展示，即：畫點(diǎn)、加樣本名，加圖例。分別用了plot(),text(),legend()，三個(gè)函數(shù)。
plot函數(shù)
plot函數(shù)是最基本的畫圖函數(shù)，感覺網(wǎng)上說的挺全的，這里引用https://blog.csdn.net/sl_world/article/details/80864674和https://blog.csdn.net/eric_e/article/details/83243999，自己也收藏一下。
本例中直接用的是PCA中的PC1和PC2，也就是pcaData[,1:2]。col顏色，是根據(jù)分組加的，我感覺是根據(jù)pcaData中的順序，因?yàn)槲抑霸囘^：
plot(pcaData[,1:2],pch = 19,col= c("red","green","blue"),cex=1)
效果如圖4
在這圖中，點(diǎn)的顏色完全是根據(jù)pcaData的順序加的，如下
PC1 PC2 group condition name
X10_T -4.9487193 -30.3912003 III III X10_T
X11_T 23.7381223 -3.9537594 II II X11_T
X12_T -56.8346035 -14.8049233 III III X12_T
X77_T -1.7273395 29.0222667 II II X77_T
X14_T -14.2359870 1.8616294 II II X14_T
X76_T -0.2251326 27.9834964 II II X76_T
可以對(duì)比pcaData和圖4，完全就是紅綠藍(lán)，紅綠藍(lán)這種，而不是按照condition劃分的。
所以在我的代碼中，我先把pcaDatacondition排序了，然后根據(jù)table知道每一個(gè)組有幾個(gè)，再設(shè)置顏色。
text函數(shù)
推薦參考http://www.mamicode.com/info-detail-1774107.html講的還是很詳細(xì)的，也給自己收藏一下。比較重要的參數(shù)是text(x,y,label)是最基本的，x,y分別是確定位置；label是文本的內(nèi)容；另外還可以通過adj,pos對(duì)文字的位置進(jìn)行調(diào)整。
legend函數(shù)
推薦參考http://www.itdecent.cn/p/004f00142d00，其實(shí)一直很想把圖例放在圖的外面，現(xiàn)在還沒有做成。
參考：http://www.itdecent.cn/p/e03af1169e57
http://www.itdecent.cn/p/53258de21108
