从大型数据队列中推断分泌蛋白的信号活性，并计算它们与临床信息的相关性。输入的表达值应来自 RNAseq或MicroArray数据，并需经过log2(x+1) 转换，其中x可以是标准化的数据如FPKM、RPKM或TPM值。演示数据来自一个胰腺癌免疫治疗队列，使用SecAct分析，识别与临床结局相关的分泌蛋白。
推断活性
使用 SecAct.activity.inference 函数在患者样本中推断分泌蛋白的信号活性。如果数据中有对照样本，将它们赋值给inputProfile_control。SecAct将使用对照样本对每个基因的表达值进行标准化，即减去对照样本中该基因的平均表达值。如果没有对照样本，只需设置 inputProfile_control = NULL。SecAct 将对每个基因在所有输入样本中的表达值进行标准化，即将所有样本中该基因的平均表达值视为对照。
library(SecAct)
# prepare expression matrix
dataPath <- file.path(system.file(package="SecAct"), "extdata/")
expr <- read.table(paste0(dataPath,"Pancreatic_Nivolumab_Padron2022.logTPM.gz"), check.names=F)
res <- SecAct.activity.inference(inputProfile = expr, inputProfile_control = NULL)
act <- res$zscore
act[1:3,1:5]
2 3 6 8 9
A1BG -25.291225 -26.150454 -12.0577423 17.828415 -12.3510028
A2M 21.686456 12.635013 36.6278568 6.109273 -18.8522668
A2ML1 -1.085858 -7.682554 -0.7080762 5.670982 0.7868586
推断出的蛋白活性是一个相对指标，正值表示在相应患者中活性较高，负值表示活性较低。
计算临床相关性
接下来，通过将分泌蛋白的活性水平与临床数据关联，计算每个分泌蛋白的信号风险评分。对于临床文件，请确保前两列名称分别为Time和Event。SecAct.coxph.regression 将进行Cox 比例风险 (PH) 回归，风险评分以双侧Wald检验的z-scores (Coef / StdErr) 表示。如果存在临床协变量如年龄、性别、分期，也会被纳入回归分析。
clinical <- read.table(paste0(dataPath, "Pancreatic_Nivolumab_Padron2022.OS_Nivo+Sotiga+Chemo"))
head(clinical, 3)
Time Event Age Gender ECOG
8 834 0 55 1 0
13 28 1 60 0 1
15 273 1 69 0 0
riskMat <- SecAct.coxph.regression(act, clinical)
head(riskMat, 3)
risk score z p value
A1BG -2.550539 0.01075565
A2M -2.375125 0.01754302
A2ML1 1.686931 0.09161658
risk score z正值表示该分泌蛋白与较差的生存相关，负值表示与较好的生存相关。
可视化风险评分
可以使用 SecAct.lollipop.plot 函数可视化感兴趣的分泌蛋白。这里分别选取与较差和较好生存结局相关的风险评分最高和最低的前10个分泌蛋白作为示例。
high.risk.SPs <- names(sort(riskMat[, "risk score z"], decreasing = TRUE))[1:10]
low.risk.SPs <- names(sort(riskMat[, "risk score z"]))[1:10]
SPs <- c(high.risk.SPs, low.risk.SPs)
fg.vec <- riskMat[SPs, "risk score z"]
SecAct.lollipop.plot(fg.vec, title = "risk score z")
绘制生存曲线
可以选择一个感兴趣的分泌蛋白，使用SecAct.survival.plot函数绘制其生存曲线。活性阈值是通过最大化高活性与低活性患者组之间的差异来确定的。
riskMat["WNT7B",]
risk score z p value
3.4497516697 0.0005611024
SecAct.survival.plot(act, clinical, "WNT7B", x.title = "Overall (Days)")
