1简介

在本例中，我们将使用临床数据和三种类型的“组学数据”对乳腺肿瘤进行二元分类。我们还使用了一些策略和相似度定义来创建特征。

为此，我们将使用来自癌症基因组图谱的数据，并将整合四种类型的组学数据:

• 来自安捷伦mRNA芯片的基因表达
• DNA甲基化(Illumina HumanMethylation 27K微阵列)
• 反相蛋白质阵列的蛋白质组学测量
• microrna的测序

图1显示了将患者数据转换为相似网络的规则，相似网络作为模型的输入单位(或“特征”)。

• 基因表达:特征是在水平上定义的通路；即，一个特征组对应通路的基因。相似度定义为成对的皮尔森相关
• 临床变量:各变量是将其自身的特征和相似度定义为归一化差异．
• 蛋白质组学和甲基化数据:在整体水平上定义特征数据层；为所有的蛋白质组数据创建一个单一的特征，甲基化也是如此。相似度是两两定义的皮尔森相关

图1:预测设计

构建和测试预测器的整体工作流程如图2所示。

我们从一个包含295个肿瘤的数据集开始。90%的样本用于构建预测器(buildPredictor ())，而10%的样本是随机抽样的，并用于独立验证(subsampleValidationData ())．构建预测器包括:

1. 将样本分成80%的训练样本和20%的测试样本(比例可通过设置参数来改变)buildPredictor ()）
2. 使用训练样本运行特征选择，使特征得分从0到用户指定的最大值。
3. 采用通过截点的特征对20%的测试样本进行分类。

这个过程重复多次随机训练/测试分割，以生成模型在训练数据上的平均性能。选择在所有训练/测试分段中得分一致的特征来构建最终模型。

然后，模型在独立验证集(这里是10%的样本)上进行评估，仅使用始终如一的高分特征(即通过选择的特征)。这一步确定了模型的测试集性能。

图2:工作流

在本例中，我们构建了一个最小预测器，以确保在可行的时间内构建小插图。这里我们使用了两个火车/测试分割，从0-2给特征打分，调用特征在给定的火车/测试分割中为>0个选定的特征打分。为了构建最终的模型，我们选择在至少一半的分割中得分为1或更高的特征。这些是不现实的参数，因为它们没有充分地重新采样训练数据以确保可泛化。

在实践中，合理的值包括10次训练/测试分割，特征得分至少为10分，最终模型的特征得分为7分或更高。尝试几种设计，看看哪一种一致地泛化到一个独立的测试集(超参数调优)。

2设置

加载netDx包中。

suppressWarnings (suppressMessages(需要(netDx)))

3.获取和准备数据

在这个例子中，我们使用了来自BioConductor的癌症基因组图谱(The Cancer Genome Atlas)的多模态数据curatedTCGAData包中。TCGA中所有癌症类型的数据都可通过此包获得;有关详细信息，请参阅本教程．

suppressMessages(库(curatedTCGAData))

让我们来看看乳腺癌的现有数据，不下载任何(集)dry.run = TRUE)．请注意，从BioC 3.13起，用户可以获取两个版本的TCGA数据。

curatedTCGAData (diseaseCode =“BRCA化验=“*”,dry.run = TRUE, version = " 1.1.38 ")

## snapshotDate(): 2022-04-19

##看到'?curatedTCGAData'用于'疾病编码'和'检测'输入

# # ah_id标题file_size # # 1 EH584 brca_cnaseq Mb - 20160128 0 # # 2 EH585 brca_cnasnp - 20160128 9.8 Mb # # 3 EH586 brca_cnvsnp - 20160128 2.8 Mb # # 4 EH588 brca_gistic_allbygene - 20160128 1.3 Mb # # 5 EH2121 brca_gistic_peaks Mb - 20160128 0 # # 6 EH589 brca_gistic_thresholdedbygene - 20160128 0.4 Mb # # 7 EH2122 brca_methylation_methyl27 - 20160128 _assays 63.2 Mb # # 8 EH2123 brca_methylation_methyl27 - 20160128 _se 0.4 Mb # # 9 EH2124 brca_methylation_methyl450 - 20160128 _assays 2613.2 Mb # # 10 EH2125brca_methylation_methyl450 - 20160128 _se 6.1 Mb # # 11 EH593 brca_mirnaseqgene - 20160128 0.6 Mb # # 12 EH594 brca_mrnaarray - 20160128 27.3 Mb # # 13 EH595 brca_mutation - 20160128 4.5 Mb # # 14 EH596 brca_rnaseq2genenorm - 20160128 64.5 Mb # # 15 EH597 brca_rnaseqgene 30 Mb - 20160128 # # 16 EH598 brca_rppaarray - 20160128 1.6 Mb # # rdataclass rdatadateadded rdatadateremoved # # 1 RaggedExperiment 2017-10-10 < NA > # # 2 RaggedExperiment 2017-10-10 < NA > # # 3 RaggedExperiment 2017-10-10 < NA > # # 4 SummarizedExperiment2017-10-10 < NA > # # 5 RangedSummarizedExperiment 2017-10-10 < NA > # # 6 SummarizedExperiment 2017-10-10 < NA > # # 7 SummarizedExperiment 2017-10-10 < NA > # # 8 SummarizedExperiment 2019-01-09 < NA > # # 9 RaggedExperiment 2019-01-09 < NA > # # 10 SummarizedExperiment 2019-01-09 < NA > # # 11 SummarizedExperiment 2019-01-09 < NA > # # 12 SummarizedExperiment 2017-10-10 < NA > # # 13 RaggedExperiment 2017-10-10 < NA > # # 14 SummarizedExperiment 2017-10-10 < NA > # # 15 SummarizedExperiment 2017-10-10 < NA > # # 16总结实验2017-10-10  .

在这个调用中，我们只获取基因表达、蛋白质组和甲基化数据;设置dry.run = FALSE开始获取数据。

brca <- suppressMessages(curatedTCGAData(" brca "， c("mRNAArray"，"miRNA*"，"Methylation_methyl27*")， dry.run=FALSE,version="1.1.38"))

这个调用返回一个MultiAssayExperiment对象。回想一下，这是一个容器，用于存储对同一组样品进行的多个分析。参见本教程了解更多。

brca

MultiAssayExperiment对象，包含3个实验，用户自定义名称和各自的类。## [1] BRCA_miRNASeqGene-20160128: summarizeexperiment with 1046 row and 849 columns ## [2] BRCA_mRNAArray-20160128: summarizeexperiment with 17814 row and 590 column ## [3] BRCA_Methylation_methyl27-20160128: summarizeexperiment with 27578 row and 343 column ##功能:## experiments() -获取ExperimentList实例## colData() -主/表型DataFrame ## sampleMap() -样本协调DataFrame ## ' $ '， '['， '[[' -提取colData列，子集，或实验## *格式()-转换为长或宽的DataFrame ## assays() -转换ExperimentList为矩阵的SimpleList ## exportClass() -保存数据到平面文件

总结(化验(brca))

##[1]“长度为3，包含0个元数据列的列表对象”

名称(化验(brca))

"BRCA_miRNASeqGene-20160128" "BRCA_mRNAArray-20160128" "BRCA_Methylation_methyl27-20160128"

mir <- assays(brca)[["BRCA_miRNASeqGene-20160128"]] head(mir[，1:5])

# # # # TCGA-3C-AAAU-01A-11R-A41G-13 TCGA-3C-AALI-01A-11R-A41G-13 hsa-let-7a-1 95618 49201 # # hsa-let-7a-2 189674 98691 # # hsa-let-7a-3 96815 49035 # # hsa-let-7b 264034 148591 # # hsa-let-7c 3641 5095 # # hsa-let-7d 4333 3263 # # TCGA-3C-AALJ-01A-31R-A41G-13 TCGA-3C-AALK-01A-11R-A41G-13 # # hsa-let-7a-1 75342 57278 # # hsa-let-7a-2 150472 114320 # # hsa-let-7a-3 76206 57540 # # hsa-let-7b 99938 164553 # # hsa-let-7c 5799 18464 # # hsa-let-7d 5658 2114 # # # # hsa-let-7a-1 TCGA-4H-AAAK-01A-12R-A41G-13 67196## hsa-let-7a-3 67607 ## hsa-let-7b 136918 ## hsa-let-7c 20429 ## hsa-let-7d 2162

pheno <- colData(brca) colnames(pheno)[1:20]

# # # #[1]“patientID”[2]# #“years_to_birth”[3]“vital_status”# #[4]“days_to_death”# #[5]“days_to_last_followup”# #[6]“tumor_tissue_site”# #[7]“pathologic_stage”# #[8]“pathology_T_stage”# #[9]“pathology_N_stage”# #[10]“pathology_M_stage”# #[11]“性别”# #[12]“date_of_initial_pathologic_diagnosis”# #[13]“days_to_last_known_alive”# #[14]“radiation_therapy”# #[15]“histological_type”# #[16]“number_of_lymph_nodes”# #[17]“种族”# #[18]“种族”# #[19]”管理。Bcr " ## [20] "admin.day_of_dcc_upload"

(把[1:5])

## 6行5列数据表## patientID年数出生生命期状态日数死亡## <字符>    ## TCGA-A1-A0SD TCGA-A1-A0SD 59 0 NA ## TCGA-A1-A0SE TCGA-A1-A0SE 56 0 NA ## TCGA-A1-A0SF TCGA-A1-A0SF 54 0 NA ## TCGA-A1-A0SG TCGA-A1-A0SG 61 0 NA ## TCGA-A1-A0SH TCGA-A1-A0SH 39 0 NA ## days_to_last_followup ##整数> ## TCGA-A1-A0SB 259 ## TCGA-A1-A0SD 1321 ## TCGA-A1-A0SF 1463 ## TCGA-A1-A0SG 434## TCGA-A1-A0SH 1437

source("prepare_data.R") brca <- prepareData(brca,setBinary=TRUE)

##协调输入:##删除28个带“colname”的sampleMap行，而不是实验的colname

##协调输入:##删除44个带有“colname”的sampleMap行，而不是实验的colname

##协调输入:##删除19个带有“colname”的sampleMap行，而不是实验的colname

pID <- colData(brca)$patientID colData(brca)$ID <- pID

set.seed(123) dsets <- subsampleValidationData(brca,pctValidation=0.1) brca <- dsets$trainMAE holdout <- dsets$validationMAE

4创建特征的规则(患者相似性网络)

现在让我们设置输入到netDx的数据。

netDx允许用户定义如何将数据转换为患者相似性网络(psn)，这是模型中的特征。这是通过告诉模型如何完成的:

• 集团不同类型的数据和
• 定义相似对于每一个(例如皮尔逊相关，归一化差异等)。

相关的输入参数为:

• groupList:对应于单个网络的输入数据集(例如按路径分组的基因)
• 西姆斯:为每个数据层指定相似度指标的列表

mRNA数据中的groupList <- list() #基因按路径pathFile <- sprintf("%s/extdata/pathway_ex3. txt ")进行分组。gmt"， path.package("netDx")) pathList <- suppressMessages(readPathways(pathFile)) groupList[["BRCA_mRNAArray-20160128"]] <- pathList . gmt"， path.package("netDx")

总结(groupList)

长度类模式BRCA_mRNAArray-20160128 3 -none- list

名称(groupList[[“brca_mrnaarray - 20160128]])

##[2] "硬脂酸生物合成动物" ##[2]"色氨酸降解生物" ## [3]"

长度(groupList[[“brca_mrnaarray - 20160128]] [[1]])

## [1]

头(groupList[[“brca_mrnaarray - 20160128]] [[1]])

##[1]“elovl1”“acot7”“acsl1”“acsl5”“elovl6”“acsl4”

groupList[["临床"]]<- list(年龄="病人。age_at_initial_pathologic_diagnosis"， stage=" stage ")

TMP <- list(rownames(实验(brca)[[1]]));names(tmp) <- names(brca)[1] groupList[[names(brca)[[1]]]] <- tmp tmp <- list(rownames(experiments(brca)[[2]]));names(tmp) <- names(brca)[2] groupList[[names(brca)[2]]] <- tmp tmp <- list(rownames(experiments(brca)[[3]]));names(tmp) <- names(brca)[3] groupList[[names(brca)[3]]] <- tmp

模拟人生<- list(a="pearsonCorr"， b="normDiff"， c="pearsonCorr"， d="pearsonCorr") #将图层名称映射到模拟人生名称(sims) <- names(groupList)

5建立预测

现在我们准备训练我们的模型。netDx使用并行处理来加快计算时间。在本例中，让我们使用机器上75%的可用内核。如果提供了一个已经有内容的输出目录，netDx也会抛出一个错误，所以让我们也来清理一下。

nco <- round(parallel::detectCores()*0.75) #使用75%可用内核消息(sprintf("使用%i of %i cores"， nco, parallel::detectCores()))

72核中的54核

outDir <- paste(tempdir()，"pred_output"，sep=getFileSep()) #使用绝对路径if (file.exists(outDir)) unlink(outDir,recursive=TRUE) num< - 2L

最后，我们调用构建预测器的函数。我们提供:

• 病人资料(dataList）
• 分组规则(groupList）
• 函数从数据创建PSN，包括相似度度量的选择(makeNetFunc）
• 收集特征分数和平均表现的训练/测试分段数(numSplits)，
• 在一轮特征选择中特征的最高得分(featScoreMax，设置为10)
• 为每个列车/测试分割调用特征选择网络的阈值(featSelCutoff）;只有得分为这个值或更高的特征才会被用于对测试患者进行分类，
• 用于并行处理的核数(numCores)．

下面的调用运行两个火车/测试分段。在每个分裂中，它:

• 使用默认的80:20分割数据到train/test (trainProp = 0.8）
• 在0到2之间给网络打分(即:featScoreMax l = 2）
• 使用评分为>=1的网络(featSelCutoff l = 1)来对该分裂的测试样本进行分类。

在实践中，一个好的起点是featScoreMax = 10，featSelCutoff = 9而且numSplits = 10 l，但这些参数取决于数据集中的样本大小和样本的异质性。

t0 < - Sys.time模型()< - suppressMessages (buildPredictor (dataList = brca, # #数据groupList = groupList, # #分组策略模拟人生= sims, outDir = outDir, # #输出目录trainProp = 0.8, # # pct的样本用于训练模型# #每一个分割numSplits = 2 l, # #的列车数量/测试分裂featSelCutoff = 1 l,东西称之为# # # #阈值精选featScoreMax = 2 l, # #马克斯得分为特征选择numCores = nco, # #设置更高的并行debugMode = FALSE, keepAllData = FALSE,##设置为TRUE调试日志="none" ##设置为"default"消息))t1 <- Sys.time()打印(t1-t0)

时差14.21651分钟

6检查结果

我们现在用getResults ()获取各种列车/测试分割的模型性能以及特征分数:

results <- getResults(model, unique(colData(brca)$STATUS)， featureSelCutoff=2L, featureSelPct=0.50)

检测到2个分裂和2个职业

## *绘图性能

## *编译特性分数并调用所选特性

结果包含性能，selectedFeatures为每个病人标签，并表特征分数．

总结(结果)

长度类模式##选择功能2 -none- list ##功能评分2 -none- list ##性能4 -none- list

看看表现:

结果美元性能

## [1] 0.8402542 ## ## $splitAccuracy ## [1] 0.8500000 0.8305085 ## ## $splitAUROC ## [1] 1 1 ## ## $splitAUPR ## [1] 0.9767442 0.9761905

看看所有标签的特征分数，在所有的火车测试中:

结果featureScores美元

# #腔的美元。# # 1 # #特性Split1 Split2 brca_methylation_methyl27 - 20160128 1 1 # # 2 brca_mrnaarray - 20160128 2 2 # # 3 brca_mirnaseqgene - 20160128 NA NA 1 # # 4年龄2 # # 5第1阶段NA其他# # # # # # $特性Split1 Split2 # # 1 brca_methylation_methyl27 - 20160128 2 2 # # 2 brca_mrnaarray - 20160128 2 2 # # 3 brca_mirnaseqgene - 20160128 1岁1 # # 4 2 2 # # 5第二阶段2

让我们来看看我们的困惑矩阵:

confMat <- confusionMatrix(model)

注意:这个矩阵的行相加不等于100%，因为这个矩阵是来自所有训练/测试分割的混淆矩阵的平均值。

这里是精选的特征，是那些在至少一半的分割中得到2分的特征。这个阈值只是为了说明。在实践中，我们将运行至少10次训练/测试分割(理想情况是100+)，并寻找在>70%分割中得分为7+的特征。

结果selectedFeatures美元

# #腔的美元。##[1]“BRCA_mRNAArray-20160128”“年龄”## ## $other #[1]“BRCA_Methylation_methyl27-20160128”“BRCA_mRNAArray-20160128”##[3]“年龄”“阶段”

7在独立样本上验证

现在我们用预测()对独立数据集中的样本进行分类。为模型提供特征设计规则groupList，要在其中使用的选定功能的列表featSelNet，将数据转换为患者相似度网络的功能makeNets，以及中原始和验证的数据集brca而且坚持分别。

需要提供训练数据，因为netDx使用训练数据和测试数据创建了单个患者相似度网络。然后，它使用标签传播将患者标签从训练样本“扩散”到测试样本，并根据它们与哪个类别最相似来标记后者。

outDir <- paste(tempdir()， randAlphanumString()， sep = getFileSep()) if (file.exists(outDir)) unlink(outDir，递归=TRUE) dir.create(outDir) predModel <- suppressMessages(predict(trainMAE=brca, testMAE=holdout, groupList=groupList, selectedFeatures=results$selectedFeatures, sims=sims, outDir=outDir, verbose = FALSE))

测试列车运行状态。其他16类状态管腔。另一个## Luminal。A 1 ##其他0 16

8验证图结果

最后，我们检查了我们的模型执行得如何getPerformance ()．

计算性能:

- getPerformance(predModel, unique(colData(brca)$STATUS)) summary(perf)

##长度类模式## rocCurve 1性能S4 ## prCurve 1性能S4 ## auroc 1 -none-数值## auroc 1 -none-数值##精度1 -none-数值

我们绘制AUROC和AUPR曲线使用plotPerf_multi ()．在这个例子中，我们得到了两个类的完美分离。

plotPerf_multi(list(perf$rocCurve)， plotTitle = sprintf("BRCA验证:%i个样本"，nrow(colData(holdout))))

plotPerf_multi(list(perf$prCurve)， plotType = "PR"， plotTitle = sprintf("BRCA验证:%i样本"，nrow(colData(holdout))))

##这个调用在Rstudio中不起作用;现在我们已经把这个注释掉并保存了PSN文件。psn <- suppressMessages(getPSN(brca, groupList, sims=sims, selectedFeatures=results$selectedFeatures))

##目录下的警告。create(paste(netDir， "profiles"， sep = fsep)): '/tmp/Rtmp52iamk/ ## profiles' already exists

##目录下的警告。create(paste(netDir， "INTERACTIONS"， sep = fsep)): '/tmp/ ## Rtmp52iamk/INTERACTIONS'已经存在

- tSNEPlotter(psn$patientSimNetwork_unpruned, colData(brca))

## *制作对称矩阵

## *运行tSNE

## *绘图

9sessionInfo

sessionInfo ()

## R版本4.2.0 RC (2022-04-19 r82224) ##平台:x86_64-pc-linux-gnu(64位)##运行在Ubuntu 20.04.4 LTS ## ##矩阵产品:默认## BLAS: /home/biocbuild/bbs-3.15-bioc/R/lib/libRblas。/home/biocbuild/bbs-3.15-bioc/R/lib/libRlapack。so ## ## locale: ## [1] LC_CTYPE=C LC_NUMERIC= c# # [3] LC_TIME=C LC_COLLATE= c# # [5] LC_MONETARY=C LC_MESSAGES=en_US。UTF-8 ## [7] LC_PAPER=en_US。UTF-8 LC_NAME= c# # [9] LC_ADDRESS=C lc_phone = c# # [11] LC_MEASUREMENT=en_US。UTF-8 LC_IDENTIFICATION=C ## ##附加的基本包:## [1]stats4 stats graphics grDevices utils datasets methods ##[8]基础## ##其他附加包:[7] MultiAssayExperiment_1.22.0 SummarizedExperiment_1.26.0 ## [9] Biobase_2.56.0 GenomicRanges_1.48.0 ## [11] GenomeInfoDb_1.32.0 IRanges_2.30.0 ## [13] S4Vectors_0.34.0 BiocGenerics_0.42.0 ## [15] MatrixGenerics_1.8.0 matrixStats_0.62.0 ## [17] netDx_1.8.0 BiocStyle_2.24.0 ## ##通过命名空间加载(并且没有附加):# # # # [1] bigmemory_4.5.36 bigmemory.sri_0.1.3 [3] colorspace_2.0-3 ellipsis_0.3.2 # # [5] XVector_0.36.0 farver_2.1.0 # # [7] bit64_4.0.5 interactiveDisplayBase_1.34.0 # # [9] AnnotationDbi_1.58.0 fansi_1.0.3 # # [11] codetools_0.2-18 splines_4.2.0 # # [13] doParallel_1.0.17 cachem_1.0.6 # # [15] knitr_1.38 jsonlite_1.8.0 # # [17] dbplyr_2.1.1 png_0.1-7 # # [19] shiny_1.7.1 BiocManager_1.30.17 # # [21] compiler_4.2.0 httr_1.4.2 # # [23] assertthat_0.2.1 fastmap_1.1.0 # # [25] cli_3.3.0 later_1.3.0 # # [27]htmltools_0.5.2 tools_4.2.0 ## [29] igraph_1.3.1 gtable_0.3.0 ## [31] glue_1.6.2 GenomeInfoDbData_1.2.8 ## [33] reshape2_1.4.4 dplyr_1.0.8 ## [35] rappdirs_0.3.3 Rcpp_1.0.8.3 ## [37] jquerylib_0.1.4 rhdf5filters_1.8.0 ## [39] vctrs_0.4.1 Biostrings_2.64.0 ## [41] ExperimentHub_2.4.0 iterators_1.0.14 ## [43] xfun_0.30 stringr_1.4.0 ## [45] mime_0.12 lifecycle_1.0.1 ## [47] AnnotationHub_3.4.0 zlibbioc_1.42.0 ## [49] scales_1.2.0 promises_1.2.0.1 ## [51] parallel_4.2.0 RColorBrewer_1.1-3 ## [53] yaml_2.3.5 curl_4.3.2 ## [55] memoise_2.0.1 ggplot2_3.3.5 ## [57] sass_0.4.1 stringi_1.7.6 ## [59] RSQLite_2.2.12 BiocVersion_3.15.2 ## [61] highr_0.9 foreach_1.5.2 ## [63] plotrix_3.8-2 filelock_1.0.2 ## [65] shape_1.4.6 rlang_1.0.2 ## [67] pkgconfig_2.0.3 bitops_1.0-7 ## [69] pracma_2.3.8 evaluate_0.15 ## [71] lattice_0.20-45 ROCR_1.0-11 ## [73] purrr_0.3.4 Rhdf5lib_1.18.0 ## [75] labeling_0.4.2 bit_4.0.4 ## [77] tidyselect_1.1.2 plyr_1.8.7 ## [79] magrittr_2.0.3 bookdown_0.26 ## [81] R6_2.5.1 generics_0.1.2 ## [83] combinat_0.0-8 DBI_1.1.2 ## [85] pillar_1.7.0 withr_2.5.0 ## [87] survival_3.3-1 KEGGREST_1.36.0 ## [89] RCurl_1.98-1.6 tibble_3.1.6 ## [91] crayon_1.5.1 utf8_1.2.2 ## [93] BiocFileCache_2.4.0 rmarkdown_2.14 ## [95] grid_4.2.0 blob_1.2.3 ## [97] digest_0.6.29 xtable_1.8-4 ## [99] httpuv_1.6.5 munsell_0.5.0 ## [101] glmnet_4.1-4 bslib_0.3.1