意外之喜之不完全的喜-MOVICS

发表于2022-01-26更新于2022-10-21

本来以为要自己把脚本都搞一遍，结果这两天出现了MOVICS这个包，简单看了一下，联合IOBR，基本所有需要的脚本都能直接用了，省下不少时间，当然，也有需要改造的部分，老规矩先走一遍示例数据。

注意，在这个包的示例里面也说到了这个问题：

“Notably, all these omics data used for clustering share exactly the same samples with exactly the same order, and please make sure of this when using MOVICS on your own data.”

所有数据都统一了sample以及完全相同的“顺序”，尤其是顺序问题，在实际分析中非常容易搞错，特别是涉及到cbind/rbind/left_join/right_join/inner_join数据框之后，如果结果明显感觉不对，八成就是因为合并的时候顺序出了问题，这里建议统一使用固定主键的inner_join，基本上这个是不会出问题的。

结果在研究了一上午之后……

我他喵直接不写这个的教程了，因为这玩意只能用在多组学上面，而且这个包对比分析的时候用的是PAM50，明明这玩意是多组学的，用单组学的结果进行对比，但是居然不支持单组学的分析，感觉有点离谱，就不是很想写。

但是这个包的vignette里写了一个compSurv的函数，可以用来很方便地写多组学的生存分析，有值得借鉴的地方。

原本这个教程里面写的是这样：

1
2
3
4
5
6
7
8
9
10
11
# survival comparison
surv.brca <- compSurv(moic.res         = cmoic.brca,
                      surv.info        = surv.info,
                      convt.time       = "m", # convert day unit to month
                      surv.median.line = "h", # draw horizontal line at median survival
                      xyrs.est         = c(5,10), # estimate 5 and 10-year survival
                      fig.name         = "KAPLAN-MEIER CURVE OF CONSENSUSMOIC")
#> --a total of 643 samples are identified.
#> --removed missing values.
#> --leaving 642 observations.
#> --cut survival curve up to 10 years.

但是我测试了一下，提取list的subset也能进行分析，能得到一样的结果，只需要提供cmoic.brca中的clust.res这个就行，这是一个data.frame，rownames是样本名，然后还有一个就是cluster。理论上就只需要这两个就行了。

但是在我自己构造数据的时候发现不行，会报错，经过查看，初步锁定是因为没有获取到list里面的样本名，无法和生存矩阵匹配上。

最后在源码里面发现了这个：

1
2
comsam <- intersect(rownames(surv.info), rownames(moic.res$clust.res))
mosurv.res <- cbind.data.frame(surv.info[comsam, c("futime", "fustat")], moic.res$clust.res[comsam, "clust", drop = FALSE])

嗯，虽然但是，这其实是个反面的代码案例就是了，直接将列名写死了，虽然对整个包的功能有帮助，但是灵活性就失去了，而且这里既然已经用上了data.frame的一些使用方法，不过为什么必须要调原本的list对象呢….这他喵感觉就是个典型的反面代码案例了……

因此，如果想在自己的数据里用这个函数，只需要构造名字为clust.res的list就行了。

1
testlist <- list("clust.res" = clust.data.frame)

so…….看起来MOVICS的本身流程的确适用面比较窄了。但是既然它接收我们自己改造的函数，那么就有一些能纳入我们自己的分析流了。

同样的，下面一些功能函数也可以用类似的方法进行替代执行，只需要更换对应的数据即可。

比如这个跑药物敏感性的，就能直接拿来用。

1
2
3
4
5
6
drug.brca <- compDrugsen(moic.res    = iClusterBayes.res,
                         norm.expr   = fpkm[,iClusterBayes.res$clust.res$samID], # double guarantee sample order
                         drugs       = c("Cisplatin","Paclitaxel"),
                         tissueType  = "breast",
                         test.method = "nonparametric",
                         prefix      = "BOXVIOLIN OF ESTIMATED IC50")

还有差异分析的。但是颜色居然要手动设定，简直有点hp，太麻烦了。因为调用的是complexheatmap，所以颜色只能用circlize::colorRamp2传进去。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
runDEA(dea.method = "deseq2",
       expr       = count,
       moic.res   = cmoic.brca,
       prefix     = "TCGA-BRCA")

# extract PAM50, pathologic stage and age for sample annotation
annCol    <- surv.info[,c("PAM50", "pstage", "age"), drop = FALSE]

# generate corresponding colors for sample annotation
annColors <- list(age    = circlize::colorRamp2(breaks = c(min(annCol$age),
                                                           median(annCol$age),
                                                           max(annCol$age)), 
                                                colors = c("#0000AA", "#555555", "#AAAA00")),
                  PAM50  = c("Basal" = "blue",
                            "Her2"   = "red",
                            "LumA"   = "yellow",
                            "LumB"   = "green",
                            "Normal" = "black"),
                  pstage = c("T1"    = "green",
                             "T2"    = "blue",
                             "T3"    = "red",
                             "T4"    = "yellow", 
                             "TX"    = "black"))


# choose edgeR result to identify subtype-specific up-regulated biomarkers
marker.up <- runMarker(moic.res      = cmoic.brca,
                       dea.method    = "deseq2", # name of DEA method
                       prefix        = "TCGA-BRCA", # MUST be the same of argument in runDEA()
                       dat.path      = getwd(), # path of DEA files
                       res.path      = getwd(), # path to save marker files
                       p.cutoff      = 0.05, # p cutoff to identify significant DEGs
                       p.adj.cutoff  = 0.05, # padj cutoff to identify significant DEGs
                       dirct         = "up", # direction of dysregulation in expression
                       n.marker      = 100, # number of biomarkers for each subtype
                       doplot        = TRUE, # generate diagonal heatmap
                       norm.expr     = fpkm, # use normalized expression as heatmap input
                       annCol        = annCol, # sample annotation in heatmap
                       annColors     = annColors, # colors for sample annotation
                       show_rownames = FALSE, # show no rownames (biomarker name)
                       fig.name      = "UPREGULATED BIOMARKER HEATMAP")
#> --all samples matched.
#> --log2 transformation done for expression data.

所以大概总结一下这里能用在自己的常规分析流里面的：

compClinvar:创建样本表并进行统计分析，用于得到三线表，类似tableone

compSurv：subtype的生存分析，输出总体p值和两两之间的生存差异

compDrugsen：subtype药物敏感性预测，用的是GCP2016的数据，也就是pRRophetic

runDEA：subtype的差异分析，可选deseq2，limma和edger，是下游功能分析的前置

runMarker：subtype的marker选择，有这none overlap的特性

runGSEA：subtype的GSEA分析，none overlap且可选up和down

后续会出把这些脚本纳入其中的代码，大概…..下下次.jpg

在做了在做了（咕咕咕）。