DESeq2 PCA 的一些问题

近日,做差异分析的时候,想着看一下样本本身的特征是以什么分类的,除了计算样本之间的距离,还用到的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")
图1 只画点

图2 加了样本名

图3 加上图例

下面说一下我画这个图的艰辛的过程。

首先,给PCA上每一个点加样本名这个过程其实是分步骤进行的。也就是我三个图的展示,即:画点、加样本名,加图例。分别用了plot(),text(),legend(),三个函数。

plot函数

plot函数是最基本的画图函数,感觉网上说的挺全的,这里引用https://blog.csdn.net/sl_world/article/details/80864674https://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


图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

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 205,132评论 6 478
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 87,802评论 2 381
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 151,566评论 0 338
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,858评论 1 277
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,867评论 5 368
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,695评论 1 282
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 38,064评论 3 399
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,705评论 0 258
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 42,915评论 1 300
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,677评论 2 323
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,796评论 1 333
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,432评论 4 322
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 39,041评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,992评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,223评论 1 260
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 45,185评论 2 352
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,535评论 2 343

推荐阅读更多精彩内容