基于支持向量机模型的TNBC的分子亚型预测
连续两次求贤令：曾经我给你带来了十万用户，但现在祝你倒闭，以及 生信技能树知识整理实习生招募，让我走大运结识了几位优秀小伙伴！有做ngs实战整理的，也有做临床数据挖掘算法工具介绍的。今天分享的是复旦大学和西北民族大学小伙伴合作的笔记
下面是笔记原文
1. 引入
最近刚好毕业，闲来无事，偶然看见曾老师的知识整理的邀请，故应邀加入知识整理的队伍。TNBC.CMS [1] R包是曾老师安排的第一个整理任务。
TNBC.CMS: prediction of TNBC consensus molecular subtype TNBC共识分子亚型的预测
[1] Kim J, Yu D, Kwon Y, et al. Genomic Characteristics of Triple-Negative Breast Cancer Nominate Molecular Subtypes That Predict Chemotherapy Response. Mol Cancer Res. 2020;18(2):253-263. doi:10.1158/1541-7786.MCR-19-0453
2. TNBC.CMS 开发背景
TNBC.CMS是一个用于三阴性乳腺癌(TNBC)分子亚型分类的包。虽然已有了很多分类策略，但缺乏准确的亚型分类方法仍是患者诊断和TNBC研究的局限。作者的基于机器学习的分类器模型使用957名TNBC患者的基因表达谱。
TNBC.CMS包将患者分为四种常见的共识分子亚型(consensus molecular subtypes，CMS)：mesenchymal-like(MSL)、immunomodulatory(IM)、luminal AR(LAR)和stem-like(SL)。它还提供了基因组和临床特征的概要，包括生存率、风险比、通路活性和药物反应
后续原作者在2021年基于前述的TNBC四种亚型，探索了根据 SMARCA4活性调节的通路及Cisplatin Resistance Score顺铂耐药评分，发表于Cancers[2]。（这个蛮有意思的，值得读）
[2] Kim J, Jang G, Sim SH, Park IH, Kim K, Park C. SMARCA4 Depletion Induces Cisplatin Resistance by Activating YAP1-Mediated Epithelial-to-Mesenchymal Transition in Triple-Negative Breast Cancer. Cancers (Basel). 2021;13(21):5474. Published 2021 Oct 30. doi:10.3390/cancers13215474
在初步了解TNBC.CMS的开发背景后，我也对TNBC.CMS产生了一些兴趣，我好奇TNBC.CMS的分类方法是否具有普及性，能不能够其他肿瘤中进行推广，完成对其他肿瘤的分型。
下面将尝试重复其分类步骤，并尝试加以解读。
3. TNBC.CMS 安装
作者将TNBC.CMS存储在Biocondutor， 因此我们可以使用以下的方式安装：
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("TNBC.CMS")如想直接查看作者给定的说明文档，可以参考以下代码：
browseVignettes("TNBC.CMS")4. TNBC.CMS使用
下面根据作者提供的说明文档，进行试用
4.1 加载R包
library(TNBC.CMS)
## 如不想看到过多加载信息 可使用
## suppressMessages(library(TNBC.CMS))4.2 加载数据
data("GSE25055") #获取R包中提供的示例数据
dim(assays(GSE25055)[[1]]) #查看数据维度 dim == Dimension
## [1] 4746 73 ##表明 有4746个genes 和76个 samples
assays(GSE25055)[[1]][1:5, 1:5] ## 查看数据内容
## 615397 615396 615394 615393 615392
## HBA1 15.42644 16.03556 16.01673 15.33820 12.66187
## HBB 14.91473 15.70969 15.72610 14.80730 12.30998
## B2M 12.91446 13.98487 13.79634 12.62554 11.59076
## FTH1 11.98103 12.57083 12.08568 12.98890 12.16884
## HLA-B 12.49575 14.57368 14.37568 12.32062 10.64294
class(GSE25055) ## 查看数据的类型
## [1] "SummarizedExperiment"
attr(,"package")
## [1] "SummarizedExperiment"
View(GSE25055)
可以看到TNBC.CMS 提供的是一个 “SummarizedExperiment” 格式的数据类型，在assays(GSE25055)[1]中，行名是gene symbol，列名是样本名，数据为表达值。
4.3 示例1：CMS分类（基于共识分子亚型）
`precdictCMS`函数根据输入矩阵或`SummarizedExperiment`对象将共识分子亚型分配给TNBC样本。如果输入是一个`SummarizedExperiment`对象，分析列表中的第一个元素应该是基因表达矩阵。在任何情况下，基因表达谱都不应该被缩放（原文：Scaled ）或对数转换 (原文：log-transformed )。可以通过访问`probabilities`属性来检索分类概率。
predictions <- predictCMS(GSE25055)
## Warning message:
## In class(object) <- "environment" :
## 把class(x)设成"environment"集空属性; 其结果不再是S4目标对象
table(predictions)
## predictions
## MSL IM LAR SL
## 12 14 23 24
head(attr(predictions, "probabilities"))
## MSL IM LAR SL
## 615397 0.004965739 0.003456084 0.006844380 0.98473380
## 615396 0.026430368 0.014144690 0.939662186 0.01976276
## 615394 0.056331795 0.020242679 0.020387143 0.90303838
## 615393 0.001031452 0.001748296 0.956579694 0.04064056
## 615392 0.002081453 0.038822011 0.791781322 0.16731521
## 615390 0.029205717 0.008909284 0.004956605 0.95692839这里四个分类即是前述的几个亚群：
间充质样干细胞样亚群（mesenchymal-like stem-like，MSL）
免疫调节亚群（ immunomodulatory，IM）
管腔雄激素受体亚群（luminal androgen receptor，LAR）
干细胞样亚群（stem-like，SL）
实际上到这里根本看不明白这四个群是怎么来的，只知道，我丢了一个表达矩阵进去，然后这个函数就自动将各个样本分别属于各个群的概率算了出来。按上述结果描述，编号“615397”属于SL的概率为 ‘0.98473380’，因此判定为 SL。
因此为了更好的了解这一步干了什么，有必要查看一下 “predictCMS” 函数的源代码
predictCMS
function (expr)
{
if (is(expr, "SummarizedExperiment")) {
exp.mat <- assays(expr)[[1]]
}
else {
exp.mat <- expr
}
missings <- genelist[!genelist %in% rownames(exp.mat)]
nmissing <- length(missings)
if (nmissing > 0) {
missing.mat <- matrix(NA, nrow = nmissing, ncol = ncol(exp.mat))
rownames(missing.mat) <- missings
colnames(missing.mat) <- colnames(exp.mat)
input <- rbind(exp.mat, missing.mat)
input <- impute(input, what = "median")
}
else {
input <- exp.mat
}
input <- log2(input + 1)
input <- input - apply(input, 1, median)
input <- data.frame(t(input[genelist, ]), check.names = FALSE)
pred <- predict(SVM.model, input, probability = TRUE)
prob <- attr(pred, "probabilities")
attr(pred, "probabilities") <- prob[, c("MSL",
"IM", "LAR", "SL")]
return(pred)
}
<bytecode: 0x000001b9473dc4a8>
<environment: namespace:TNBC.CMS>可以看到，出现“input”变量前面的步骤均是在对数据进行清洗处理，后续的关键的分类步骤是使用SVM.model完成的，然后就直接出现了"MSL", "IM", "LAR", "SL" 几个我们关注的因素， 而前面的 genelist 和 SVM.model 都不是在我们环境下输入的参数。
genelist
## 错误: 找不到对象'genelist'
SVM.model
## 错误: 找不到对象'SVM.model'因此，合理的推测是 作者先将给定genelist 在这个表达矩阵中筛选出来，然后通过已构建好的SVM.model将各个样本进行预测，而这个genelist和SVM.model均被封装在这个R包，在运行这个函数的时候调用了给定的genelist 和 SVM.model。
那下一步我想试试能不能在R包给的文档里面找到这个genelist 和SVM.model。
但……我没找到，先这样吧，就不浪费过多时间了，希望有同学会的可以指导下，想学习这个genelist。
（小编寄语：其实超级简单， 使用3个冒号即可，代替之前的两个冒号的操作符。）
其实就是实习生的编程基础知识不过关，具备基础的计算机知识非常重要，我把它粗略的分成基于R语言的统计可视化，以及基于Linux的NGS数据处理：
但大概意思估摸就是各种给定每组类型匹配一个 genelist，然后根据给定基因的表达量在SVM.model下去预测这个样本分别属于哪个类型的。
4.4 示例2：基因组和临床特征的概要
`TNBC.CMS`包有几个函数用于研究共识分子亚型的基因组和临床特征。在本节中，我们将这些函数应用于GSE25055基因表达和临床特征的数据集。
4.4.1 ComputeGES
`ComputeGES`函数计算以下7个基因表达特征的特征分数：EMT(epithelial-mesencymal transition上皮-间充质转化)、stromal基质、免疫、微环境、stemness干性、hormone激素和CIN(chro- mosomal instability染色体不稳定性)。有关基因表达特征的更多详细信息，请参见`ComputeGES`函数手册页面。
resultGES <- computeGES(expr = GSE25055, pred = predictions,
+ rnaseq = FALSE)
## Warning message:
## In .gsva(expr, mapped.gset.idx.list, method, kcdf, rnaseq, abs.ranking, :
## Some gene sets have size one. Consider setting 'min.sz > 1'.
resultGES[,1:4]
## 615397 615396 615394 615393
## EMT 6.0372043 6.1889247 6.3529614 5.7296889
## Stromal -605.4905027 -642.1497283 -552.4131551 -929.4231951
## Immune 167.0953013 282.4093345 293.8429060 -141.3963687
## Microenvironment 0.3019017 0.2577112 0.3392774 0.2170629
## Stemness 0.5490508 0.3003425 0.3994718 0.6053706
## Hormone 4.9628294 4.9485363 5.5935250 5.4526575
## CIN 7.3907780 6.7571045 7.2265374 7.5212038在运行的过程中，同时生成了如下所示图形
如图1所示，此函数绘制了特征分数的箱线图，并用亚型之间的p值进行比较。MSL亚型、IM亚型、LAR亚型和SL亚型的Stromal基质、immune免疫、hormone激素和stemness干性得分分别显著高于其他亚型。
照例，我们也来看一下这个函数进行了什么操作
computeGES
function (expr, pred, rnaseq = FALSE)
{
if (is(expr, "SummarizedExperiment")) {
exp.mat <- assays(expr)[[1]]
}
else {
exp.mat <- expr
}
CMS.palette <- c(MSL = "brown2", IM = "gold2",
LAR = "yellowgreen", SL = "midnightblue")
pred <- factor(pred, levels = c("MSL", "IM",
"LAR", "SL"))
EMT.score <- colMeans(exp.mat[rownames(exp.mat) %in% EMT.geneset,
])
CIN.score <- colMeans(exp.mat[rownames(exp.mat) %in% CIN.geneset,
])
hormone.score <- colMeans(exp.mat[rownames(exp.mat) %in%
c("AR", "ESR1", "ERBB2", "PGR"),
])
X <- exp.mat[rownames(exp.mat) %in% names(Stemness.geneset),
]
Stemness.geneset = Stemness.geneset[rownames(X)]
s <- apply(X, 2, function(z) {
cor(z, Stemness.geneset, method = "sp", use = "complete.obs")
})
s <- s - min(s)
stemness.score <- s/max(s)
ESTIMATE.score <- computeESTIMATEscore(exp.mat)
stromal.score <- unlist(ESTIMATE.score["StromalSignature",
])
immune.score <- unlist(ESTIMATE.score["ImmuneSignature",
])
microenvironment.score <- computexCellScore(exp.mat, rnaseq)
sig.mat <- rbind(EMT.score, stromal.score, immune.score,
microenvironment.score, stemness.score, hormone.score,
CIN.score)
rownames(sig.mat) <- c("EMT", "Stromal", "Immune",
"Microenvironment", "Stemness", "Hormone",
"CIN")
……
## 后续计算P值及画图部分省略可以发现，这个函数的评分方式很简单粗暴，比如EMT.score评分，实际上就是选取EMT相关基因，然后在各个样本中取ColMeans 即作为他的评分，个人感觉还不如GSVA评分或者用AUCell进行排序，单细胞的addmodule 感觉也比这个靠谱。免疫基质相关评分采用的是 Estimate 的比较经典的方式。但这里有个问题是，在函数运行的时候 出现一个gsva函数的warning，但是函数代码没有看到gsva，也可能是我没有理解到位，持保留意见。
EMT.score <- colMeans(exp.mat[rownames(exp.mat) %in% EMT.geneset,])4.4.2 performGSVA
`performGSVA`函数对基因集进行基因集差异分析，并生成代表GSVA富集分数的热图。如果基因组没有给出，则使用标记通路基因组。用户还可以通过设置`gsva.kcdf`参数来选择用于估计表达式值的累积分布函数的要点，该参数默认设置为"`Gaussian`"。如果表达值是整数，建议使用"`Poisson`"。
resultGSVA <- performGSVA(expr = GSE25055, pred = predictions,
gene.set = NULL)
head(resultGSVA[,1:4])
这部分就不展开了，就是一个常见的GSVA的富集分析。
4.4.3 plotKM 和 plotHR
`TNBC.CMS`包提供了两个生存分析函数：`plotKM`和`plotHR`。在这里，我们使用来自GSE25055数据集的生存数据来研究总生存率和共识分子亚型之间的关系。生存数据也包含在`SummarizedExperiment`对象中，可以使用`colData`函数进行访问。
time <- colData(GSE25055)$DFS.month
event <- colData(GSE25055)$DFS.status
plotKM(pred = predictions, time = time, event = event)
`plotKM`函数为每个共识分子亚型生成Kaplan-Meier曲线，如图3所示。SL组的预后最差，这与作者前述研究是一致的。
> library(survival)
> #Test for difference of survival between low and high expression groups
> surv <- Surv(time, event)
> GSE25055.exprs <- assays(GSE25055)[[1]]
> chisq <- apply(GSE25055.exprs, 1, function(x) survdiff(surv ~ (x > median(x)))$chisq)
> pval <- 1 - pchisq(chisq, 1)
> #Select 10 genes with lowest p-values for the log-rank test
> gs <- names(sort(pval)[1:10])
> gs
[1] "RECK" "RELN" "EHD4" "PRRX2" "FOLR1" "UGCG" "GOT2" "PRPF39" "OPHN1" "CPE"
> ## [1] "RECK" "RELN" "EHD4" "PRRX2" "FOLR1" "UGCG" "GOT2" "PRPF39"
> ## [9] "OPHN1" "CPE"
> plotHR(expr = GSE25055, gene.symbol = gs, pred = predictions, time = time,
+ event = event, by.subtype = FALSE)
也可以设置 “by.subtype= TRUE“ 按亚群进行计算HR
plotHR(expr = GSE25055, gene.symbol = gs[1:4], pred = predictions,
time = time, event = event, by.subtype = TRUE)
4.5 示例3：药物反应研究
`TNBC.CMS` 包还提供了预测药物反应的函数。`computeDS` 函数计算 `MSigDB CGP`（chemical and genetic perturbations化学和遗传扰动）集合中相应基因集的药物特征分数，并绘制特征分数的热图。药物特征分数计算为基因集平均表达值之间的差值，这些基因集与药物反应和耐药性相关的。分数越高，患者越有可能引起反应。用户可以通过`gene.set` 参数提供他们自己的基因集。请注意，基因组的名称必须遵循 `[DRUG NAME]_[RESPONSE/RESISTANCE]_[UP/DN]` 的格式（例如 `CISPLATIN_RESISTANCE_UP`）。
图6显示了每个样本药物特征分数的热图。MSL和SL亚型似乎分别对dasatinib 达沙替尼 和doxorubicin阿霉素耐药。此外，IM和LAR亚型对 androgen agonist雄激素激动剂和SB216763（GSK3B的抑制剂）的得分高于其他亚型。
resultDS <- computeDS(expr = GSE25055, pred = predictions)
head(resultDS[,1:4])
## 615397 615396 615394 615393
## APLIDIN -0.8876518 -0.4680677 -0.8102091 -0.8549382
## CISPLATIN 0.9589450 0.6949340 1.2613627 0.2520613
## DASATINIB 0.5460543 0.7005523 -0.1239162 0.7189546
## FORSKOLIN 6.4843470 6.7434440 6.4589256 6.6777797
## IMATINIB -6.3512657 -5.0914031 -5.4834761 -5.2843840
## TROGLITAZONE 6.8838127 6.7156754 6.7170530 6.6795939
5. 结果保存
以上是TNBC.CMS提供的几个函数的主要用法，最后输出的结果作者也提供了一个保存方式，以供后续使用。
dfCMS <- data.frame(row.names = colnames(GSE25055.exprs), CMS = predictions, t(resultGES),
t(resultDS), stringsAsFactors = FALSE)
head(dfCMS)
## CMS EMT Stromal Immune Microenvironment Stemness Hormone CIN APLIDIN CISPLATIN DASATINIB FORSKOLIN IMATINIB
## 615397 SL 6.037204 -605.4905 167.09530 0.3019017 0.5490508 4.962829 7.390778 -0.8876518 0.9589450 0.5460543 6.484347 -6.351266
## 615396 LAR 6.188925 -642.1497 282.40933 0.2577112 0.3003425 4.948536 6.757104 -0.4680677 0.6949340 0.7005523 6.743444 -5.091403
## 615394 SL 6.352961 -552.4132 293.84291 0.3392774 0.3994718 5.593525 7.226537 -0.8102091 1.2613627 -0.1239162 6.458926 -5.483476
## 615393 LAR 5.729689 -929.4232 -141.39637 0.2170629 0.6053706 5.452658 7.521204 -0.8549382 0.2520613 0.7189546 6.677780 -5.284384
## 615392 LAR 5.576038 -1518.6248 -560.49739 0.1650646 1.0000000 4.224790 8.239136 -1.3667877 0.3644356 0.8224894 6.703497 -4.930696
## 615390 SL 6.233911 -607.0961 31.71176 0.2690540 0.5903492 4.602222 7.540700 -0.8707087 0.7678301 0.1427109 6.795645 -5.007474
## TROGLITAZONE TZD SB216763 ANDROGEN TAMOXIFEN BEXAROTENE DOXORUBICIN
## 615397 6.883813 7.214105 6.142319 3.949582 -6.029159 0.3264693 -6.875380
## 615396 6.715675 5.341876 6.268575 4.039061 -6.283494 0.1351034 -6.320782
## 615394 6.717053 6.388377 5.994192 3.848465 -6.292139 -1.0281784 -6.709809
## 615393 6.679594 6.290939 5.844824 3.895471 -5.917298 0.1482893 -6.984418
## 615392 7.042860 5.494772 5.731897 4.047578 -6.097761 -0.1288242 -7.581800
## 615390 6.856523 7.490999 6.045182 4.173174 -6.205849 -0.8771287 -6.7876716. TNBC.CMS试用总结在最开始拿到这个TNBC.CMS的时候，我以为他是一个无监督聚类生成四个分组，然后作者根据各个分组的基因表达模式分为四个亚群。但是在逐步探索后发现，实际上这个R包还是一个根据既定分组，既定亚型相关基因做出的区分，就相当于，我先划好4个区域，再把这些样本根据不同区域的特征分选出来，分选用的方法就是SVM。因此，虽然R包命名为TNBC.CMS，但实际上和TNBC没啥关系，换个肿瘤，套入这个模式，大概率也能分成四组。
腾讯云开发者
