近日,做差异分析的时候,想着看一下样本本身的特征是以什么分类的,除了计算样本之间的距离,还用到的PCA(主成分分析)。在DESeq2包中专门由一个PCA分析的函数,即plotPCA,里面的参数也比较简单。
plotPCA参数
object:对象
在进行PCA之前是要进行归一化的,因为在降维映射的过程中, 存在映射误差, 所有在对高维特征降维之前, 需要做特征归一化(feature normalization), 这个归一化操作包括: (1) feature scaling (让所有的特征拥有相似的尺度, 要不然一个特征特别小, 一个特征特别大会影响降维的效果) (2) zero mean normalization (零均值归一化)。、
但是在DESeq2包中实际上已经有了归一化的方法,rlog和vst,具体是怎么做到归一化的,这个原理不了解。
在使用的需要根据样本量的多少来选择方法。样本量少于30的话,选择rlog,多于30的话,建议选择vst。亲测!样本量多于30的时候,用了rlog,出现以下提示:
rlog() may take a few minutes with 30 or more samples,
vst() is a much faster transformation
所以官方就会建议样本量大于30的时候用vst,而且vst速度是飞起的状态,而rlog是老牛拉慢车。
intgroup:分组
官方解释:interesting groups: a character vector of names in colData(x) to use for grouping。也就是在构建dds的时候的分组,实际上这个分组最后影响的是PCA图中的颜色,但是并不影响PCA图中各个样本的位置。换句话说,PCA降维的是时候并不会考虑这个分组。这个也是我搞了好多天才弄清楚的,所以原理还是要搞懂的,不能只会跑流程。
ntop:用于主成分分析的时候基因数
官方解释:number of top genes to use for principal components, selected by highest row variance。也就是说PCA分析的时候,是根据基因表达的counts值来的。如果给出了ntop,那么意思就是只选出表达量高的这几个基因去计算。没想通作者为什么设置这个,我觉得应该是先筛选出差异的基因,然后再设置这个数会比较好,这样就是看筛选出的基因是不是符合预期。反正感觉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
这个参数的作用就是当你的PCA图需要加样本名的时候,也就是你希望在PCA图上知道哪个点是哪个样本,你就需要导出这个。
其实我的初衷也就是想看在PCA图上哪个点是哪个样本,因此我当然要设置returnData=T。
#构建dds
dds <- DESeqDataSetFromMatrix(countData=indata, colData = state, design= ~ condition )
#归一化,因为样本量超过了30,因此用vst
vsd <- vst(dds, blind = FALSE)
#返回样本名和分组
pcaData <- plotPCA(vsd, intgroup=c("condition"),returnData = T)
#这里按照condition排序了,原因见下。
pcaData <- pcaData[order(pcaData$condition,decreasing=F),]
#知道每一个组有多少样本
table(pcaData$condition)
# II III IV
# 20 11 10
#根据上面的结果,设置每一个组的数量,方便加颜色;这一步是把每一个样本根据分组情况画出来,效果见图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")
下面说一下我画这个图的艰辛的过程。
首先,给PCA上每一个点加样本名这个过程其实是分步骤进行的。也就是我三个图的展示,即:画点、加样本名,加图例。分别用了plot(),text(),legend(),三个函数。
plot函数
plot函数是最基本的画图函数,感觉网上说的挺全的,这里引用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颜色,是根据分组加的,我感觉是根据pcaData中的顺序,因为我之前试过:
plot(pcaData[,1:2],pch = 19,col= c("red","green","blue"),cex=1)
效果如图4
在这图中,点的颜色完全是根据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
可以对比pcaData和图4,完全就是红绿蓝,红绿蓝这种,而不是按照condition划分的。
所以在我的代码中,我先把pcaDatacondition排序了,然后根据table知道每一个组有几个,再设置颜色。
text函数
推荐参考http://www.mamicode.com/info-detail-1774107.html讲的还是很详细的,也给自己收藏一下。比较重要的参数是text(x,y,label)是最基本的,x,y分别是确定位置;label是文本的内容;另外还可以通过adj,pos对文字的位置进行调整。
legend函数
推荐参考https://www.jianshu.com/p/004f00142d00,其实一直很想把图例放在图的外面,现在还没有做成。
参考:https://www.jianshu.com/p/e03af1169e57
https://www.jianshu.com/p/53258de21108