--- title: "CorNetto COVID-19 Severity Workflow" author: "Bradley Ward" date: "`r Sys.Date()`" output: BiocStyle::html_document: toc: true vignette: > %\VignetteIndexEntry{2. CorNetto COVID-19 Severity Workflow} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set(collapse = TRUE, comment = "#>") library(CorNetto) extdataDir <- system.file("extdata", package = "CorNetto") ``` # Introduction This vignette is an applied companion to `vignette("CorNetto")`, which introduces the package, motivates it against related Bioconductor software, and covers installation and the full function reference. Read that one first. Here the same workflow is run end to end on a small packaged subset of proteomic, transcriptomic, and metabolomic data from moderate and severe COVID-19 patients. The subset exists to demonstrate the workflow during package checks; it is not sized for biological discovery, and its provenance is documented in `inst/scripts/covid-example-data.md`. The candidate edges used below are entirely synthetic. They exercise network import, filtering, differential correlation, and node scoring, but they do not represent biological evidence. Constraining the calculations to an explicit candidate network keeps this example fast and avoids dense all-pairs correlations. # Import synthetic candidate networks CorNetto accepts prior-knowledge networks from external resources when their columns can be mapped to the package schema. For this runnable example, the package instead supplies three deterministic synthetic networks. This avoids redistributing third-party interaction data and keeps the example independent of database versions and licenses. The within-assay and cross-assay files contain invented edges between measured features. The decoy file contains invented edges with one or two unmeasured endpoints, so the filtering step below has observable work to do. Edge weights and directions are illustrative only. `readKnowledgeNetwork()` function takes a delimited table, optionally remaps the input column names to a common convention, and standardizes it for downstream CorNetto functions. The minimum columns required are `fromFeatureIdentifier`, `toFeatureIdentifier`, `fromAssayName`, and `toAssayName`. Other columns are optional but can be useful. The output is a standardized `S4Vectors::DataFrame`. `combineKnowledgeNetworks()` takes these standardized networks and combines their rows while removing duplicate edges (optional). ```{r load-prior-data} # Set column mapping priorColumnMap <- c( fromFeatureIdentifier = "from", toFeatureIdentifier = "to", fromFeatureName = "fromName", toFeatureName = "toName", fromAssayName = "fromOmic", toAssayName = "toOmic", edgeType = "edgeType", edgeDirection = "edgeDirection", evidenceScore = "edgeWeight" ) # Import the three synthetic edge sets withinAssayNetwork <- readKnowledgeNetwork( filePath = file.path( extdataDir, "priorNetworks", "synthetic_within_assay.csv" ), columnMapping = priorColumnMap, knowledgeSource = "CorNetto synthetic within-assay" ) crossAssayNetwork <- readKnowledgeNetwork( filePath = file.path( extdataDir, "priorNetworks", "synthetic_cross_assay.csv" ), columnMapping = priorColumnMap, knowledgeSource = "CorNetto synthetic cross-assay" ) decoyNetwork <- readKnowledgeNetwork( filePath = file.path( extdataDir, "priorNetworks", "synthetic_decoy_edges.csv" ), columnMapping = priorColumnMap, knowledgeSource = "CorNetto synthetic decoys" ) # Combine networks combinedPriorNetwork <- combineKnowledgeNetworks( withinAssayNetwork, crossAssayNetwork, decoyNetwork, removeDuplicates = TRUE ) ``` # Import experimental abundance data into the environment Currently CorNetto is able to accept any normalised and fully preprocessed abundance or count matrix relating to proteomics, transcriptomics, or metabolomics. `createAnalysisData()` creates a MultiAssayExperiment from these normalized assay objects and sample metadata. `filterSamples()` filters the MultiAssayExperiment based on a group column and a range of values. ```{r load-experimental-data} # Import patient data patientInformation <- read.csv(file.path(extdataDir, "covidData", "patientInformation.csv")) rownames(patientInformation) <- patientInformation$Visit_ID # Import omics data dataList <- list(RNA = read.csv(file.path(extdataDir, "covidData", "rnaMatrix.csv"), row.names = 1), Protein = read.csv(file.path(extdataDir, "covidData", "proteinMatrix.csv"), row.names = 1), Metabolite = read.csv(file.path(extdataDir, "covidData", "metaboliteMatrix.csv"), row.names = 1)) # Create analysis data analysisData <- createAnalysisData(assayList = dataList, sampleData = patientInformation) # Filter data for only COVID-19 patients at visit 1 analysisCovid <- filterSamples(analysisData, groupColumn = "Group", groupLevels = c("COVID Severe", "COVID Moderate")) analysisCovid <- filterSamples(analysisCovid, groupColumn = "Visit", groupLevels = "Visit 1") ``` # QC checks of all imported objects Following data importation, objects can be checked for structural conformity and diagnostic summaries. `validateAnalysisData()` checks the CorNetto analysis object for structural conformity, specifically: * is this a MultiAssayExperiment? * does it contain assays? * are sample identifiers valid and unique? * are assay columns matched to colData? * are feature/sample names present and non-duplicated? * are assays extractable as numeric matrices? * are required feature metadata fields present? `validateKnowledgeNetwork()` checks the CorNetto prior-knowledge edge table for structural conformity, specifically: * is this a data frame-like edge table? * are the required edge columns present? * are source and target feature identifiers present? * are source and target assay names present? * are missing optional columns filled with standard defaults? * are edge type, edge direction, and knowledge source represented consistently? * is the network converted to the standardized CorNetto edge-table format? These functions also run quietly after most CorNetto functions to ensure conformity of output. ```{r validateAnalysisData} # Check analysis data analysisCovid <- validateAnalysisData(analysisCovid) # Check knowledge networks combinedPriorNetwork <- validateKnowledgeNetwork(combinedPriorNetwork) ``` `summarizeAnalysisData()` asks "what does this analysis object look like, and are there analysis risks?" rather than only asking whether the object is structurally valid. ```{r summarizeAnalysisData} summarizeAnalysisData(analysisCovid, groupColumn = "Group") ``` The output reports the number of features per assay, the number of samples per assay, and the number of samples in each selected group. In this small packaged example there are 12 moderate COVID-19 samples and 12 severe COVID-19 samples. Zero-variance and near-zero-variance features should be removed before correlation analysis. These can be removed with `removeZeroVariance = TRUE` and `minimumVariance` in `filterFeatures()`. ```{r filterFeatures} analysisCovid <- filterFeatures( analysisCovid, removeZeroVariance = TRUE, minimumVariance = 0.01 ) summarizeAnalysisData(analysisCovid) ``` For the sparse workflow, candidate edges are restricted to edges whose endpoint features remain in the filtered analysis object. The first summary shows that the combined synthetic network includes deliberately unmeasured endpoints. The second summary shows the network after requiring both endpoints to be measured. ```{r restrict-prior-network} measuredNodeKeys <- unlist( lapply( names(MultiAssayExperiment::experiments(analysisCovid)), function(assayName) { assayObject <- MultiAssayExperiment::experiments(analysisCovid)[[assayName]] paste( assayName, rownames(SummarizedExperiment::assay(assayObject)), sep = "::" ) } ), use.names = FALSE ) unfilteredPriorSummary <- summarizeKnowledgeNetwork( combinedPriorNetwork, analysisCovid, quiet = TRUE ) unfilteredPriorSummary$overall measuredPriorNetwork <- filterNetworkByNodes( combinedPriorNetwork, nodes = measuredNodeKeys, mode = "both", nodeIdType = "key" ) measuredPriorSummary <- summarizeKnowledgeNetwork( measuredPriorNetwork, analysisCovid, quiet = TRUE ) measuredPriorSummary$overall ``` # Run prior-guided differential correlations The main function for testing differential correlations between groups is `testDifferentialCorrelation`. Here it tests every measured edge in the synthetic candidate network supplied through `candidateEdgeTable`. The result illustrates the software workflow and is not a prior-supported biological analysis. The differential correlation result is added to `analysisCovid`, then extracted with `differentialCorrelationResults()`. Alternatively, `testDifferentialCorrelation()` can directly return the table with `storeResult = FALSE`. ```{r testDifferentialCorrelation} analysisCovid <- testDifferentialCorrelation( analysisData = analysisCovid, groupColumn = "Group", groupLevels = c("COVID Moderate", "COVID Severe"), candidateEdgeTable = measuredPriorNetwork, correlationMethod = "pearson", minimumAbsoluteCorrelation = 0.3, adjustedPValueThreshold = 0.05, pAdjustMethod = "fdr", resultName = "differentialSparseCorrelations" ) differentialSparseCorrelations <- differentialCorrelationResults(analysisCovid)[["differentialSparseCorrelations"]] as.data.frame(differentialSparseCorrelations) ``` # Build differential network. The next step builds a differential correlation network from the `differentialSparseCorrelations` table. The function returns the network directly unless `storeResult = TRUE`, in which case `analysisData` and `resultName` are supplied. Here the resulting network is not filtered further. In a real analysis, `minimumAbsoluteCorrelation` and `differenceAdjustedPValueThreshold` can be tightened to remove weakly correlated or non-significant edges. ```{r createDifferentialCorrelationNetwork} analysisCovid <- createDifferentialCorrelationNetwork( differentialCorrelationTable = differentialSparseCorrelations, minimumAbsoluteCorrelation = 0, edgeWeightMethod = "signedZScore", analysisData = analysisCovid, resultName = "differentialSparseNetwork", storeResult = TRUE) differentialSparseNetwork <- differentialCorrelationNetworks(analysisCovid)[["differentialSparseNetwork"]] as.data.frame(differentialSparseNetwork) ``` # Calculate rewiring scores. Rewiring scores summarize the amount of differential-correlation signal incident on each node in the differential network. Here the rewiring table is stored in `analysisCovid` and extracted with `rewiringResults()`. ```{r calculateRewiringScores} analysisCovid <- calculateRewiringScores( differentialCorrelationNetwork = differentialSparseNetwork, analysisData = analysisCovid, resultName = "differentialSparseRewiring", storeResult = TRUE ) rewiringScores <- rewiringResults(analysisCovid)[["differentialSparseRewiring"]] rewiringScores <- rewiringScores[ order(rewiringScores$rootMeanSquareRewiringScore, decreasing = TRUE), ] head(rewiringScores, 20) ``` # Run permutation ranking. `permuteRewiringScores()` permutes the supplied group labels, reruns the differential correlation testing, rebuilds the networks, and calculates node-level permutation tail probabilities. This example retains the observed correlation and adjusted-p filters used above, so its output is a conditional ranking, not a randomization p-value. This vignette uses 99 Monte Carlo permutations so that package checks run quickly. Assignments are sampled independently and can repeat; exact enumeration is not implemented. ```{r permuteRewiringScores} analysisCovid <- permuteRewiringScores( analysisData = analysisCovid, groupColumn = "Group", groupLevels = c("COVID Moderate", "COVID Severe"), candidateEdgeTable = measuredPriorNetwork, correlationMethod = "pearson", minimumAbsoluteCorrelation = 0.3, adjustedPValueThreshold = 0.05, pAdjustMethod = "fdr", edgeWeightMethod = "signedZScore", scoreColumn = "rawRewiringScore", nPermutations = 99, seed = 1, keepPermutationScores = TRUE, resultName = "differentialSparseRewiringPermutation", storeResult = TRUE) rewiringPermutation <- validationResults(analysisCovid)[["differentialSparseRewiringPermutation"]] rewiringPermutation$inferenceStatus rankedRewiring <- rewiringPermutation$rewiringTable rankedRewiring <- rankedRewiring[ order( rankedRewiring$adjustedPermutationTailProbability, -rankedRewiring$rawRewiringScore, na.last = TRUE ), ] head(rankedRewiring, 20) ``` # Result export Once results are generated, the network and node rewiring scores can be supplied to `prepareCytoscapeTables()` to create tables that can be imported into Cytoscape for further network analysis. `writeNetworkTables()` writes the node and edge tables to a specified directory. The vignette writes to `tempdir()` so package checks do not write into the user's working directory. ```{r writeNetworkTables} cytoscapeTables <- prepareCytoscapeTables(networkEdgeTable = differentialSparseNetwork, rewiringTable = rankedRewiring) writeNetworkTables(cytoscapeTables, directoryPath = file.path(tempdir(), "cornettoCovidSeverityNetwork"), prefix = "covidModerateVsSevere", fileFormat = "csv") ``` # Session information ```{r session-info} sessionInfo() ```