设为“星标”，精彩可以过

本文主要参考官方先容：https://xlucpu.github.io/MOVICS/MOVICS-VIGNETTE.html

简介

装配

GET Module

准备数据

筛选基因（降维）

详情最好亚型数目

左证单一算法分型

同期进行多种分型算法

整合多种分型截止

稽查分型截止的质料

多组学分型热图

分子分型一直是生信数据挖掘的热点时代，用于分子分型的算法格外多，比如环球常见的非负矩阵解析、一致性聚类、PCA等，一致性聚类咱们在之前也先容过了：免疫浸润截止分子分型

今天给环球先容一个一站式的分子分型R包：MOVICS。

该包与其他分子分型R包最大的不同是它能同期使用多组学的数据，豪放的分子分型R包只可通过一种组学数据进行分析，比如只可通过mRNA的抒发矩阵进行分析。然而这R包它可以同期通过比如说mRNA、lncRNA、甲基化数据、突变数据进行分型。

以外，它还提供了分型之后每个亚型的探索以及每个亚型内的分析。是以说这是一个一站式的包。这个的功能主要分为三个部分，暗意图如下：

图片

第一个部分是左证不同的组学数据进行分型。第二个部分是比拟不同的分型。第三个部分是对每个分型进行探索，以及得回每个分型特异性的分子。

每个部分包含的主要函数如下，底下会先容：

GET Module: get subtypes through multi-omics integrative clustering

COMP Module: compare subtypes from multiple perspectives

RUN Module: run marker identification and verify subtypes

该包已发表，使用时铭记援用：

现在该包在github，只可通过以下款式装配，难得装配时最好先装配依赖包，因为这个包的依赖包格外多，装配经过中格外容易失败。对于入门者来说，这个包的装配不是很友好哦~

# 集会装配devtools::install_github("xlucpu/MOVICS")# 不详下载到土产货装配devtools::install_local("E:/R/R包/MOVICS-master.zip")

咱们先看一下示例数据。

library(MOVICS)## 

使用该包自带数据进行演示，这个自带数据是照旧清洗好的。过几天再专门写一篇推文先容怎样准备这个数据。

# TCGA的乳腺癌数据load(system.file("extdata"， "brca.tcga.RData"， package = "MOVICS"， mustWork = TRUE))load(system.file("extdata"， "brca.yau.RData"，  package = "MOVICS"， mustWork = TRUE))

brca.tcga内部是多个组学的数据，比如mRNA、lncRNA、甲基化、突变数据等，还有临床信息，比如活命时刻和活命景色以及乳腺癌的PAM50分类。

为了演示，这个数据通过MAD筛选了部分数据：

难得，这里最紧迫的少量是：每种组学的数据的样本数目、名字、划定应该透顶一致。环球可以我方看一下这些数据是什么样的。

names(brca.tcga)## [1] "mRNA.expr"   "lncRNA.expr" "meth.beta"   "mut.status"  "count"      ## [6] "fpkm"        "maf"         "segment"     "clin.info"names(brca.yau)## [1] "mRNA.expr" "clin.info"# 索要"mRNA.expr""lncRNA.expr""meth.beta""mut.status"mo.data   <- brca.tcga[1:4]# 索要raw count datacount     <- brca.tcga$count# 索要fpkm datafpkm      <- brca.tcga$fpkm# 索要mafmaf       <- brca.tcga$maf# 索要segmented copy numbersegment   <- brca.tcga$segment# 索要活命信息surv.info <- brca.tcga$clin.info

getElites，顾名想义，找出精英，找出最过劲的，也即是说这个函数可以作念一些预贬责和筛选使命，可以帮你进行数据准备使命。

主要可以作念以下预贬责：

其实这个不是第一步，第一步应该是我方先清洗一下数据，比如抒发矩阵先进行log调解等。

底下是一些功能演示，还好坏常巨大的。

缺失值插补：

# scenario 1: 贬责缺失值tmp       <- brca.tcga$mRNA.expr # get expression datadim(tmp) # check data dimension## [1] 500 643tmp[1，1]  <- tmp[2，2] <- NA # 添加几个NAtmp[1:3，1:3] # check data##         BRCA-A03L-01A BRCA-A04R-01A BRCA-A075-01A## SCGB2A2            NA          1.42          7.24## SCGB1D2         10.11            NA          5.88## PIP              4.54          2.59          4.35elite.tmp <- getElites(dat       = tmp，                       method    = "mad"，                       na.action = "rm"， # 径直删除                       elite.pct = 1) # 保留100%的数据## --2 features with NA values are removed.## missing elite.num then use elite.pctdim(elite.tmp$elite.dat) ## [1] 498 643elite.tmp <- getElites(dat       = tmp，                       method    = "mad"，                       na.action = "impute"， # 使用knn进行插补                       elite.pct = 1) ## missing elite.num then use elite.pctdim(elite.tmp$elite.dat) ## [1] 500 643elite.tmp$elite.dat[1:3，1:3] # NA values have been imputed ##         BRCA-A03L-01A BRCA-A04R-01A BRCA-A075-01A## SCGB2A2         6.867         1.420          7.24## SCGB1D2        10.110         4.739          5.88## PIP             4.540         2.590          4.35

使用MAD筛选分子：

# scenario 2: 使用MAD筛选，最大中位差tmp       <- brca.tcga$mRNA.expr elite.tmp <- getElites(dat       = tmp，                       method    = "mad"，                       elite.pct = 0.1) # 保留MAD前10%的基因## missing elite.num then use elite.pctdim(elite.tmp$elite.dat) # 500的10%是50## [1]  50 643#> [1]  50 643elite.tmp <- getElites(dat       = tmp，                       method    = "sd"，                       elite.num = 100， # 保留MAD前100的基因                       elite.pct = 0.1) # 此时这个参数就不起作用了## elite.num has been provided then discards elite.pct.dim(elite.tmp$elite.dat) ## [1] 100 643

使用PCA筛选分子，需要了解一些对于PCA的基础常识：R讲话主要素分析

# scenario 3: 使用PCA筛选分子tmp       <- brca.tcga$mRNA.expr # get expression data with 500 featureselite.tmp <- getElites(dat       = tmp，                       method    = "pca"，                       pca.ratio = 0.95) # 主要素的比例## --the ratio used to select principal component is set as 0.95dim(elite.tmp$elite.dat) # get 204 elite (PCs) left## [1] 204 643

使用单因素COX回顾筛选分子，也即是对每个分子作念单因素cox分析，弃取有益旨的留住，需要提供活命信息：

# scenario 4: 使用cox筛选分子tmp       <- brca.tcga$mRNA.expr # get expression data elite.tmp <- getElites(dat       = tmp，                       method    = "cox"，                       surv.info = surv.info， # 活命信息，列名必须有'futime'和'fustat'                       p.cutoff  = 0.05，                       elite.num = 100) # 此时这个参数亦然不起作用的## --all sample matched between omics matrix and survival data.## 5% 10% 15% 20% 25% 30% 35% 40% 45% 50% 55% 60% 65% 70% 75% 80% 85% 90% 95% 100%dim(elite.tmp$elite.dat) # get 125 elites## [1] 125 643table(elite.tmp$unicox$pvalue < 0.05) # 125 genes have nominal pvalue < 0.05 in ## ## FALSE  TRUE ##   375   125tmp       <- brca.tcga$mut.status # get mutation data elite.tmp <- getElites(dat       = tmp，                       method    = "cox"，                       surv.info = surv.info，                        p.cutoff  = 0.05，                       elite.num = 100) ## --all sample matched between omics matrix and survival data.## 7% 13% 20% 27% 33% 40% 47% 53% 60% 67% 73% 80% 87% 93% 100%dim(elite.tmp$elite.dat) # get 3 elites## [1]   3 643table(elite.tmp$unicox$pvalue < 0.05) # 3 mutations have nominal pvalue < 0.05## ## FALSE  TRUE ##    27     3

使用突变频率筛选分子，这个是准们用于0/1矩阵这种二分类数据的：

# scenario 5: 使用突变频率筛选tmp       <- brca.tcga$mut.status # get mutation data rowSums(tmp) ## PIK3CA   TP53    TTN   CDH1  GATA3   MLL3  MUC16 MAP3K1  SYNE1  MUC12    DMD ##    208    186    111     83     58     49     48     38     33     32     31 ##  NCOR1    FLG   PTEN   RYR2  USH2A  SPTA1 MAP2K4  MUC5B    NEB   SPEN  MACF1 ##     31     30     29     27     27     25     25     24     24     23     23 ##   RYR3    DST  HUWE1  HMCN1  CSMD1  OBSCN   APOB  SYNE2 ##     23     22     22     22     21     21     21     21elite.tmp <- getElites(dat       = tmp，                       method    = "freq"， # must set as 'freq'                       elite.num = 80， # 这里是指突变频率                       elite.pct = 0.1) # 此时该参数不起作用## --method of 'freq' only supports binary omics data (e.g.， somatic mutation matrix)， and in this manner， elite.pct and elite.num are used to cut frequency.## elite.num has been provided then discards elite.pct.rowSums(elite.tmp$elite.dat) # 只保留在80个及以上样本中突变的基因## PIK3CA   TP53    TTN   CDH1 ##    208    186    111     83elite.tmp <- getElites(dat       = tmp，                       method    = "freq"，                        elite.pct = 0.2) ## --method of 'freq' only supports binary omics data (e.g.， somatic mutation matrix)， and in this manner， elite.pct and elite.num are used to cut frequency.## missing elite.num then use elite.pctrowSums(elite.tmp$elite.dat) # only genes that are mutated in over than 0.2*643=128.6 ## PIK3CA   TP53 ##    208    186

左证分子抒发量对样本进行分型，分子即是上一步得到的mRNA、lncRNA、miRNA、甲基化矩阵等。

先左证CPI和Gaps-statistics详情分红几个亚型：

optk.brca <- getClustNum(data        = mo.data， # 4种组学数据                         is.binary   = c(F，F，F，T)， #前3个不是二分类的，临了一个是                         try.N.clust = 2:8，
江苏海四达电源股份有限公司 # 尝试亚型数目，
首页-科昌宝皮具有限公司从2到8                         fig.name    = "CLUSTER NUMBER OF TCGA-BRCA")#保存的文献名## calculating Cluster Prediction Index...## 5% complete## 5% complete## 10% complete## 10% complete## 15% complete## 15% complete## 20% complete## 25% complete## 25% complete## 30% complete## 30% complete## 35% complete## 35% complete## 40% complete## 45% complete## 45% complete## 50% complete## 50% complete## 55% complete## 55% complete## 60% complete## 65% complete## 65% complete## 70% complete## 70% complete## 75% complete## 75% complete## 80% complete## 85% complete## 85% complete## 90% complete## 90% complete## 95% complete## 95% complete## 100% complete## calculating Gap-statistics...## visualization done...## --the imputed optimal cluster number is 3 arbitrarily，首页-汉昌佳香料有限公司 but it would be better referring to other priori knowledge.

图片

会自动在面前使命目次下产生一个PDF形势的图片。

函数给出的截止是3，然而接洽到乳腺癌的PAM0分类，咱们弃取k=5，也即是分红5个亚型。

是以这个详情最好亚型个数是左证你我方的需要来的哈，纯真诊疗~

详情分红几个亚型之后，可以通过算法进行分型了。提供了格外多的门径，环球常见的非负矩阵解析、异质性聚类等等皆提供了。

比如左证贝叶斯门径进行分型：

# perform iClusterBayes (may take a while)iClusterBayes.res <- getiClusterBayes(data        = mo.data，                                      N.clust     = 5，                                      type        = c("gaussian"，"gaussian"，"gaussian"，"binomial")，                                      n.burnin    = 1800，                                      n.draw      = 1200，                                      prior.gamma = c(0.5， 0.5， 0.5， 0.5)，                                      sdev        = 0.05，                                      thin        = 3)## clustering done...## feature selection done...

不详使用斡旋的函数，我方弃取门径即可，两种门径得到的截止透顶是相同的：

iClusterBayes.res <- getMOIC(data        = mo.data，                             N.clust     = 5，                             methodslist = "iClusterBayes"， # 指定算法                             type        = c("gaussian"，"gaussian"，"gaussian"，"binomial")， # data type corresponding to the list                             n.burnin    = 1800，                             n.draw      = 1200，                             prior.gamma = c(0.5， 0.5， 0.5， 0.5)，                             sdev        = 0.05，                             thin        = 3)

复返的截止包含一个clust.res对象，它有两列：clust列携带样本所属的亚型，samID列记载对应的样真称号。对于提供特征弃取经过的算法（如iClusterBayes、CIMLR和MoCluster），截止还包含一个feat.res对象，存储了这种经过的信息。对于触及分层聚类的算法（举例COCA、ConsensusClustering），样本聚类的相应树状图也将行为clust.dend复返，若是用户想要将它们放在热图中会很灵验。

可以同期左证多种算法进行分型，然后整合它们的截止，得到最终的截止，不是一般的巨大：

# perform multi-omics integrative clustering with the rest of 9 algorithmsmoic.res.list <- getMOIC(data        = mo.data，                         methodslist = list("SNF"， "PINSPlus"， "NEMO"， "COCA"， "LRAcluster"， "ConsensusClustering"， "IntNMF"， "CIMLR"， "MoCluster")， # 9种算法                         N.clust     = 5，                         type        = c("gaussian"， "gaussian"， "gaussian"， "binomial"))## --you choose more than 1 algorithm and all of them shall be run with parameters by default.## SNF done...## Clustering method: kmeans## Perturbation method: noise## PINSPlus done...## NEMO done...## COCA done...## LRAcluster done...## end fraction## clustered

## ConsensusClustering done...## IntNMF done...## clustering done...## feature selection done...## CIMLR done...## clustering done...## feature selection done...## MoCluster done...

再把贝叶斯的截止整个加进来，这即是10种算法了：

服装定制 0， 0， 0.55) 0px 2px 10px;line-height: 1.6 !important;">moic.res.list <- append(moic.res.list，                         list("iClusterBayes" = iClusterBayes.res))# 保存下截止save(moic.res.list， file = "moic.res.list.rda")

鉴戒了consensus ensembles的见解，竣事对多个分型算法截止的整合。

可以画出一个一致性热图：

load(file = "moic.res.list.rda")cmoic.brca <- getConsensusMOIC(moic.res.list = moic.res.list，                               fig.name      = "CONSENSUS HEATMAP"，                               distance      = "euclidean"，                               linkage       = "average")

图片

截止会保存在面前使命目次中。

除了通过上头的热图稽查分型截止，还可以使用Silhouette准则判断分型质料。

以下是施展注解，起原于集会：

Silhouette准则是一种用于聚类分析中的评价门径，它通过对每个数据点与其所属簇内其他数据点之间的距离进行比拟，来臆测聚类质料的好坏。Silhouette准则可以匡助咱们详情最好的聚类数目，从而擢升聚类分析的可靠性和准确性。 Silhouette准则的蓄意门径如下：对于每个数据点i，蓄意它与同簇中其他数据点之间的平均距离ai，以及与最近其他簇中数据点之间的平均距离bi。然后，界说每个数据点的Silhouette统统为： s(i) = (bi - ai) / max(ai， bi) Silhouette统统的取值范围在-1到1之间，其中负值线路数据点更容易被分类到空虚的簇中，而恰恰则线路数据点更容易被正确分类。Silhouette统统的平均值可以用来评估整个聚类的质料，因此，Silhouette准则的方向是最大化Silhouette统统的平均值，从而找到最好的聚类数目。 当聚类数目增多时，Silhouette统统的平均值频繁会先增多后减少。因此，咱们需要找到一个聚类数目，使得Silhouette统统的平均值达到最大值。频繁，咱们和会过绘制Silhouette图来弃取最好的聚类数目。Silhouette图是一种以Silhouette统统为纵轴，聚类数目为横轴的图表，它可以匡助咱们直不雅地阐述聚类的质料。 在使用Silhouette准则进行聚类分析时，需要难得以下几点：

截止会保存在面前使命目次中：

getSilhouette(sil      = cmoic.brca$sil， # a sil object returned by getConsensusMOIC()              fig.path = getwd()，              fig.name = "SILHOUETTE"，              height   = 5.5，              width    = 5)

图片

## png ##   2

分型之后，坚信是要对每个组学数据进行热图展示不同亚型的抒发量情况。

不外需要作念一些准备使命。

# β值矩阵调解为M值矩阵indata <- mo.dataindata$meth.beta <- log2(indata$meth.beta / (1 - indata$meth.beta))# 对数据进行要领化plotdata <- getStdiz(data       = indata，                     halfwidth  = c(2，2，2，NA)， # no truncation for mutation                     centerFlag = c(T，T，T，F)， # no center for mutation                     scaleFlag  = c(T，T，T，F)) # no scale for mutation

咱们这里就用贝叶斯分型的截止进行展示，最初是索要每个组学的截止，然后每个组学中弃取前10个分子进行标注：

feat   <- iClusterBayes.res$feat.resfeat1  <- feat[which(feat$dataset == "mRNA.expr")，][1:10，"feature"] feat2  <- feat[which(feat$dataset == "lncRNA.expr")，][1:10，"feature"]feat3  <- feat[which(feat$dataset == "meth.beta")，][1:10，"feature"]feat4  <- feat[which(feat$dataset == "mut.status")，][1:10，"feature"]annRow <- list(feat1， feat2， feat3， feat4)

底下即是绘画即可，其实亦然借助complexheatmap竣事的，只不外帮你简化了好多经过，截止会自动保存在面前使命目次下，MOVICS的默许出图照旧很好意思不雅的，可能比你我方画的面子~

# 为每个组学的热图自界说颜料，不界说也可mRNA.col   <- c("#00FF00"， "#008000"， "#000000"， "#800000"， "#FF0000")lncRNA.col <- c("#6699CC"， "white"  ， "#FF3C38")meth.col   <- c("#0074FE"， "#96EBF9"， "#FEE900"， "#F00003")mut.col    <- c("grey90" ， "black")col.list   <- list(mRNA.col， lncRNA.col， meth.col， mut.col)# comprehensive heatmap (may take a while)getMoHeatmap(data          = plotdata，             row.title     = c("mRNA"，"lncRNA"，"Methylation"，"Mutation")，             is.binary     = c(F，F，F，T)， # the 4th data is mutation which is binary             legend.name   = c("mRNA.FPKM"，"lncRNA.FPKM"，"M value"，"Mutated")，             clust.res     = iClusterBayes.res$clust.res， # cluster results             clust.dend    = NULL， # no dendrogram             show.rownames = c(F，F，F，F)， # specify for each omics data             show.colnames = FALSE， # show no sample names             annRow        = annRow， # mark selected features             color         = col.list，             annCol        = NULL， # no annotation for samples             annColors     = NULL， # no annotation color             width         = 10， # width of each subheatmap             height        = 5， # height of each subheatmap             fig.name      = "COMPREHENSIVE HEATMAP OF ICLUSTERBAYES")

图片

上头是贝叶斯门径分型截止的展示，你也可以任选一种，毕竟咱们有10种算法。

比如弃取COCA法的截止进行展示，亦然一模相同的用法，截止会自动保存：

# comprehensive heatmap (may take a while)getMoHeatmap(data          = plotdata，             row.title     = c("mRNA"，"lncRNA"，"Methylation"，"Mutation")，             is.binary     = c(F，F，F，T)， # the 4th data is mutation which is binary             legend.name   = c("mRNA.FPKM"，"lncRNA.FPKM"，"M value"，"Mutated")，             clust.res     = moic.res.list$COCA$clust.res， # cluster results             clust.dend    = moic.res.list$COCA$clust.dend， # show dendrogram for samples             color         = col.list，             width         = 10， # width of each subheatmap             height        = 5， # height of each subheatmap             fig.name      = "COMPREHENSIVE HEATMAP OF COCA")

图片

若是你要展示多个临床信息，亦然径直添加即可，难得自界说颜料需要使用circlize竣事：

# extract PAM50， pathologic stage and age for sample annotationannCol    <- surv.info[，c("PAM50"， "pstage"， "age")， drop = FALSE]# generate corresponding colors for sample annotationannColors <- 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"))# comprehensive heatmap (may take a while)getMoHeatmap(data          = plotdata，             row.title     = c("mRNA"，"lncRNA"，"Methylation"，"Mutation")，             is.binary     = c(F，F，F，T)， # the 4th data is mutation which is binary             legend.name   = c("mRNA.FPKM"，"lncRNA.FPKM"，"M value"，"Mutated")，             clust.res     = cmoic.brca$clust.res， # consensusMOIC results             clust.dend    = NULL， # show no dendrogram for samples             show.rownames = c(F，F，F，F)， # specify for each omics data             show.colnames = FALSE， # show no sample names             show.row.dend = c(F，F，F，F)， # show no dendrogram for features             annRow        = NULL， # no selected features             color         = col.list，             annCol        = annCol， # annotation for samples             annColors     = annColors， # annotation color             width         = 10， # width of each subheatmap             height        = 5， # height of each subheatmap             fig.name      = "COMPREHENSIVE HEATMAP OF CONSENSUSMOIC")

图片

是不好坏常过劲？

到这里第一部分的现实就先容收场服装定制，底下即是探索、比拟不同的亚型了。

• 不是
• 3个
• 服装定制
• 分类