He et al 2022
Proteome Discoverer / TMT 16-plex
2023-07-19 20:35:56 CEST (einprot v0.7.4)
## The following variables were specified as input arguments when calling the rendering function.
## They will be used in the workflow below.
experimentInfo <-
list()
species <-
"human"
pdOutputFolder <-
"/Users/charlottesoneson/Documents/Rpackages/einprot/einprot_examples/he_2022/pr2c00390_si_002/TMT16_PD3"
pdResultName <-
"TMT16plex_chimerys_intensity_nocontam_220914"
inputLevel <-
"Proteins"
pdAnalysisFile <-
NULL
idCol <-
function (df)
combineIds(df, combineCols = c("Gene.Symbol", "Accession"))
labelCol <-
function (df)
combineIds(df, combineCols = c("Gene.Symbol", "Accession"))
geneIdCol <-
function (df)
getFirstId(df, colName = "Gene.Symbol")
proteinIdCol <-
"Accession"
stringIdCol <-
function (df)
combineIds(df, combineCols = c("Gene.Symbol", "Accession"), combineWhen = "missing",
makeUnique = FALSE)
modificationsCol <-
"Modifications"
excludeUnmodifiedPeptides <-
FALSE
keepModifications <-
NULL
reportTitle <-
"He et al 2022"
reportAuthor <-
"Proteome Discoverer / TMT 16-plex"
iColPattern <-
"^Abundance."
sampleAnnot <-
structure(list(sample = c("F1.127C.Sample", "F1.127N.Sample",
"F1.128C.Sample", "F1.128N.Sample", "F1.129C.Sample", "F1.129N.Sample",
"F1.130C.Sample", "F1.130N.Sample", "F1.131C.Sample", "F1.131N.Sample",
"F1.132C.Sample", "F1.132N.Sample", "F1.133N.Sample", "F2.127C.Sample",
"F2.127N.Sample", "F2.128C.Sample", "F2.128N.Sample", "F2.129C.Sample",
"F2.129N.Sample", "F2.130C.Sample", "F2.130N.Sample", "F2.131C.Sample",
"F2.131N.Sample", "F2.132C.Sample", "F2.132N.Sample", "F2.133N.Sample",
"F3.127C.Sample", "F3.127N.Sample", "F3.128C.Sample", "F3.128N.Sample",
"F3.129C.Sample", "F3.129N.Sample", "F3.130C.Sample", "F3.130N.Sample",
"F3.131C.Sample", "F3.131N.Sample", "F3.132C.Sample", "F3.132N.Sample",
"F3.133C.Sample", "F3.133N.Sample", "F4.127C.Sample", "F4.127N.Sample",
"F4.128C.Sample", "F4.128N.Sample", "F4.129C.Sample", "F4.129N.Sample",
"F4.130C.Sample", "F4.130N.Sample", "F4.131C.Sample", "F4.131N.Sample",
"F4.132C.Sample", "F4.132N.Sample", "F4.133N.Sample"), group = c("Medulla_COVID",
"Cortex_COVID", "Medulla_COVID", "Cortex_COVID", "Medulla_COVID",
"Cortex_COVID", "Medulla_COVID", "Cortex_COVID", "Medulla_Control",
"Cortex_COVID", "Cortex_Control", "Cortex_Control", "Medulla_Control",
"Medulla_COVID", "Cortex_COVID", "Medulla_COVID", "Cortex_COVID",
"Medulla_COVID", "Cortex_COVID", "Medulla_COVID", "Cortex_COVID",
"Medulla_Control", "Cortex_Control", "Medulla_Control", "Cortex_Control",
"Cortex_Control", "Medulla_COVID", "Cortex_COVID", "Medulla_COVID",
"Cortex_COVID", "Medulla_COVID", "Cortex_COVID", "Medulla_COVID",
"Cortex_COVID", "Medulla_COVID", "Cortex_COVID", "Medulla_Control",
"Cortex_Control", "Medulla_Control", "Cortex_Control", "Medulla_COVID",
"Cortex_COVID", "Medulla_COVID", "Cortex_COVID", "Cortex_COVID",
"Cortex_COVID", "Cortex_Control", "Medulla_COVID", "Cortex_Control",
"Medulla_Control", "Cortex_Control", "Medulla_Control", "Medulla_Control"
), batch = c("b1", "b1", "b1", "b1", "b1", "b1", "b1", "b1",
"b1", "b1", "b1", "b1", "b1", "b2", "b2", "b2", "b2", "b2", "b2",
"b2", "b2", "b2", "b2", "b2", "b2", "b2", "b3", "b3", "b3", "b3",
"b3", "b3", "b3", "b3", "b3", "b3", "b3", "b3", "b3", "b3", "b4",
"b4", "b4", "b4", "b4", "b4", "b4", "b4", "b4", "b4", "b4", "b4",
"b4")), class = "data.frame", row.names = c(NA, -53L))
includeOnlySamples <-
c("F1.127C.Sample", "F1.127N.Sample", "F1.128C.Sample", "F1.128N.Sample",
"F1.129C.Sample", "F1.129N.Sample", "F1.130C.Sample", "F1.130N.Sample",
"F1.131C.Sample", "F1.131N.Sample", "F1.132C.Sample", "F1.132N.Sample",
"F1.133N.Sample", "F2.127C.Sample", "F2.127N.Sample", "F2.128C.Sample",
"F2.128N.Sample", "F2.129C.Sample", "F2.129N.Sample", "F2.130C.Sample",
"F2.130N.Sample", "F2.131C.Sample", "F2.131N.Sample", "F2.132C.Sample",
"F2.132N.Sample", "F2.133N.Sample", "F3.127C.Sample", "F3.127N.Sample",
"F3.128C.Sample", "F3.128N.Sample", "F3.129C.Sample", "F3.129N.Sample",
"F3.130C.Sample", "F3.130N.Sample", "F3.131C.Sample", "F3.131N.Sample",
"F3.132C.Sample", "F3.132N.Sample", "F3.133C.Sample", "F3.133N.Sample",
"F4.127C.Sample", "F4.127N.Sample", "F4.128C.Sample", "F4.128N.Sample",
"F4.129C.Sample", "F4.129N.Sample", "F4.130C.Sample", "F4.130N.Sample",
"F4.131C.Sample", "F4.131N.Sample", "F4.132C.Sample", "F4.132N.Sample",
"F4.133N.Sample")
excludeSamples <-
""
minScore <-
2
minDeltaScore <-
0.20000000000000001
minPeptides <-
1
minPSMs <-
2
masterProteinsOnly <-
FALSE
imputeMethod <-
"MinProb"
assaysForExport <-
NULL
mergeGroups <-
list()
comparisons <-
list(c("Medulla_Control", "Medulla_COVID"), c("Cortex_Control",
"Cortex_COVID"))
ctrlGroup <-
""
allPairwiseComparisons <-
TRUE
singleFit <-
FALSE
subtractBaseline <-
FALSE
baselineGroup <-
""
normMethod <-
"center.median"
spikeFeatures <-
NULL
stattest <-
"limma"
minNbrValidValues <-
2
minlFC <-
0
samSignificance <-
FALSE
nperm <-
250
volcanoAdjPvalThr <-
0.050000000000000003
volcanoLog2FCThr <-
0
volcanoMaxFeatures <-
25
volcanoLabelSign <-
"both"
volcanoS0 <-
0.10000000000000001
volcanoFeaturesToLabel <-
""
addInteractiveVolcanos <-
TRUE
interactiveDisplayColumns <-
c(Label = "einprotLabel", adjP = "adj.P.Val", logFC = "logFC"
)
interactiveGroupColumn <-
NULL
complexFDRThr <-
0.10000000000000001
maxNbrComplexesToPlot <-
10
seed <-
42
includeFeatureCollections <-
"GO"
minSizeToKeepSet <-
2
customComplexes <-
list()
complexSpecies <-
"all"
complexDbPath <-
NULL
stringVersion <-
"11.5"
stringDir <-
NULL
linkTableColumns <-
c("Description", "Master", "[^se]\\.logFC$")This report describes a reproducible end-to-end analysis of a
proteomics dataset quantified with Proteome
Discoverer (Orsburn 2021). Most of the
code is hidden by default, but can be displayed by clicking on the
Code buttons (or by selecting
Code -> Show All Code in the top right corner of the
report). Navigation between the different sections can be done via the
table of contents in the left sidebar. In the first part of the report,
the quantified data is read into R and passed through a range of
processing and quality control steps. These are followed by statistical analysis to find and visualize
differentially abundant features. A summary
table provides direct links to external resources, and an additional
global overview of the data is provided via principal
component analysis.
## Get species info and define STRINGdb object
speciesInfo <- getSpeciesInfo(species)
if (inputLevel == "Proteins") {
if (is.null(stringDir)) stringDir <- ""
if (is.null(stringIdCol)) {
## If no STRING IDs are extracted, don't do STRING analysis
string_db <- NULL
} else {
string_db <- tryCatch({
tmp <- STRINGdb$new(version = stringVersion, species = speciesInfo$taxId,
score_threshold = 400, input_directory = stringDir)
if (!exists("tmp")) {
warning("The STRINGdb object can not be created. ",
"No STRING analysis will be performed.",
call. = FALSE)
tmp <- NULL
} else {
print(tmp)
}
tmp
}, error = function(e) {
warning("The STRINGdb object can not be created. ",
"No STRING analysis will be performed.",
call. = FALSE)
NULL
})
}
} else {
string_db <- NULL
}## *********** STRING - https://string-db.org ***********
## (Search Tool for the Retrieval of Interacting Genes/Proteins)
## version: 11.5
## species: 9606
## ............please wait............
## proteins: 19566
## interactions: 1795134## If needed and not provided, define path to complex DB
## (will be added to summary table below)
if ("complexes" %in% includeFeatureCollections && is.null(complexDbPath)) {
complexDbPath <- system.file(EINPROT_COMPLEXES_FILE,
package = "einprot")
}
## Get conversion tables for PomBase and WormBase IDs
pbconv <- readRDS(system.file(EINPROT_POMBASE_CONVTABLE,
package = "einprot"))
wbconv <- readRDS(system.file(EINPROT_WORMBASE_CONVTABLE,
package = "einprot"))Experiment details
makeTableFromList(c(experimentInfo,
list(
"Species" = speciesInfo$species,
"Species (common)" = speciesInfo$speciesCommon,
"Taxonomic ID" = speciesInfo$taxId
)))| Species | Homo sapiens |
| Species (common) | human |
| Taxonomic ID | 9606 |
ProteomeDiscoverer analysis summary
pd <- readProteomeDiscovererInfo(pdOutputFolder, pdResultName,
pdAnalysisFile)## Missing files: /Users/charlottesoneson/Documents/Rpackages/einprot/einprot_examples/he_2022/pr2c00390_si_002/TMT16_PD3/TMT16plex_chimerys_intensity_nocontam_220914_InputFiles.txt, /Users/charlottesoneson/Documents/Rpackages/einprot/einprot_examples/he_2022/pr2c00390_si_002/TMT16_PD3/TMT16plex_chimerys_intensity_nocontam_220914_StudyInformation.txtpdFile <- file.path(pdOutputFolder, paste0(pdResultName, "_", inputLevel, ".txt"))
pd <- c(list("PD quantification file" = pdFile), pd)
makeTableFromList(pd)| PD quantification file | /Users/charlottesoneson/Documents/Rpackages/einprot/einprot_examples/he_2022/pr2c00390_si_002/TMT16_PD3/TMT16plex_chimerys_intensity_nocontam_220914_Proteins.txt |
Settings
settingsList <- list(
"Column pattern" = iColPattern,
"Include only samples (if applicable)" = paste(includeOnlySamples,
collapse = ", "),
"Exclude samples" = paste(excludeSamples, collapse = ", "),
"Input type" = inputLevel,
"Min. number of peptides" = minPeptides,
"Min. number of PSMs" = minPSMs,
"Min. protein score" = minScore,
"Min. delta score" = minDeltaScore,
"Only retain master proteins" = masterProteinsOnly,
"Imputation method" = imputeMethod,
"Assays(s) to use for exported values" = paste(assaysForExport, collapse = ", "),
"Min. nbr valid values required for testing" = minNbrValidValues,
"Model fit" = ifelse(singleFit, "Single (one model fit for all samples)",
"Separate model fit for each comparison"),
"Groups to merge" = paste(unlist(
lapply(names(mergeGroups),
function(nm) paste0(nm, ":", paste(mergeGroups[[nm]], collapse = ",")))),
collapse = "; "),
"Comparisons" = paste(unlist(lapply(comparisons,
function(x) paste(x, collapse = " vs "))),
collapse = "; "),
"Control group" = ctrlGroup,
"Do all pairwise comparisons" = allPairwiseComparisons,
"Batch correction via baseline subtraction" = subtractBaseline,
"Baseline group" = baselineGroup,
"Normalization method" = normMethod,
"Spike features" = paste(spikeFeatures, collapse = ","),
"Statistical test" = stattest,
"Minimal fold change (limma/treat)" = minlFC,
"Adjusted p-value threshold for volcano plots" = volcanoAdjPvalThr,
"Log2 FC threshold for volcano plots" = volcanoLog2FCThr,
"Max nbr features to indicate in volcano plots" = volcanoMaxFeatures,
"Sign of features to indicate in volcano plots" = volcanoLabelSign,
"Use SAM statistic for significance" = samSignificance,
"s0" = volcanoS0,
"Features to always label in volcano plots" = paste(volcanoFeaturesToLabel,
collapse = ", "),
"Feature collections for enrichment testing" = paste(includeFeatureCollections, collapse = "; "),
"Minimal required size for feature sets" = minSizeToKeepSet,
"Complexes file" = gsub(".+\\/(.+.rds)", "\\1", complexDbPath),
"Complexes from species" = complexSpecies,
"Custom complexes" = paste(names(customComplexes), collapse = ";"),
"FDR Threshold for complexes" = complexFDRThr,
"Max nbr complexes to plot" = maxNbrComplexesToPlot,
"Number of permutations" = nperm,
"Random seed" = seed,
"Modifications column" = modificationsCol,
"Exclude unmodified peptides" = excludeUnmodifiedPeptides,
"Modifications to keep" = paste(keepModifications, collapse = ", "),
"Columns to add in link table" = paste(linkTableColumns, collapse = ";"),
"Interactive display columns" = paste(interactiveDisplayColumns, collapse = ";"),
"Interactive group column" = interactiveGroupColumn
)
if (stattest == "ttest") {
settingsList <- settingsList[!(names(settingsList) %in%
c("Minimal fold change (limma/treat)",
"Log2 FC threshold for volcano plots"))]
}
if (stattest %in% c("limma", "proDA")) {
settingsList <- settingsList[!(names(settingsList) %in%
c("s0", "Number of permutations",
"Use SAM statistic for significance"))]
}
if (inputLevel == "Proteins") {
settingsList <- settingsList[!(names(settingsList) %in%
c("Min. number of PSMs",
"Min. delta score"))]
}
if (inputLevel == "PeptideGroups") {
settingsList <- settingsList[!(names(settingsList) %in%
c("Min. number of peptides",
"Min. protein score",
"Only retain master proteins"))]
}
makeTableFromList(settingsList)| Column pattern | ^Abundance. |
| Include only samples (if applicable) | F1.127C.Sample, F1.127N.Sample, F1.128C.Sample, F1.128N.Sample, F1.129C.Sample, F1.129N.Sample, F1.130C.Sample, F1.130N.Sample, F1.131C.Sample, F1.131N.Sample, F1.132C.Sample, F1.132N.Sample, F1.133N.Sample, F2.127C.Sample, F2.127N.Sample, F2.128C.Sample, F2.128N.Sample, F2.129C.Sample, F2.129N.Sample, F2.130C.Sample, F2.130N.Sample, F2.131C.Sample, F2.131N.Sample, F2.132C.Sample, F2.132N.Sample, F2.133N.Sample, F3.127C.Sample, F3.127N.Sample, F3.128C.Sample, F3.128N.Sample, F3.129C.Sample, F3.129N.Sample, F3.130C.Sample, F3.130N.Sample, F3.131C.Sample, F3.131N.Sample, F3.132C.Sample, F3.132N.Sample, F3.133C.Sample, F3.133N.Sample, F4.127C.Sample, F4.127N.Sample, F4.128C.Sample, F4.128N.Sample, F4.129C.Sample, F4.129N.Sample, F4.130C.Sample, F4.130N.Sample, F4.131C.Sample, F4.131N.Sample, F4.132C.Sample, F4.132N.Sample, F4.133N.Sample |
| Exclude samples | |
| Input type | Proteins |
| Min. number of peptides | 1 |
| Min. protein score | 2 |
| Only retain master proteins | FALSE |
| Imputation method | MinProb |
| Assays(s) to use for exported values | |
| Min. nbr valid values required for testing | 2 |
| Model fit | Separate model fit for each comparison |
| Groups to merge | |
| Comparisons | Medulla_Control vs Medulla_COVID; Cortex_Control vs Cortex_COVID |
| Control group | |
| Do all pairwise comparisons | TRUE |
| Batch correction via baseline subtraction | FALSE |
| Baseline group | |
| Normalization method | center.median |
| Spike features | |
| Statistical test | limma |
| Minimal fold change (limma/treat) | 0 |
| Adjusted p-value threshold for volcano plots | 0.05 |
| Log2 FC threshold for volcano plots | 0 |
| Max nbr features to indicate in volcano plots | 25 |
| Sign of features to indicate in volcano plots | both |
| Features to always label in volcano plots | |
| Feature collections for enrichment testing | GO |
| Minimal required size for feature sets | 2 |
| Complexes from species | all |
| Custom complexes | |
| FDR Threshold for complexes | 0.1 |
| Max nbr complexes to plot | 10 |
| Random seed | 42 |
| Modifications column | Modifications |
| Exclude unmodified peptides | FALSE |
| Modifications to keep | |
| Columns to add in link table | Description;Master;[^se].logFC$ |
| Interactive display columns | einprotLabel;adj.P.Val;logFC |
Read PD output
The input to this workflow is a Proteins.txt file from Proteome
Discoverer (see path in the table above). We read the PD intensities
into R and store them in a SingleCellExperiment
object. This object will later be expanded with additional data, such as
transformed and imputed abundances.
At this point, the SingleCellExperiment object holds assays with the different types of intensities and annotations from the ProteomeDiscoverer file.
## Read Proteome Discoverer output
tmp <- importExperiment(inFile = pdFile, iColPattern = iColPattern,
includeOnlySamples = includeOnlySamples,
excludeSamples = excludeSamples)
sce <- tmp$sce
aName <- tmp$aName
sce## class: SingleCellExperiment
## dim: 11724 53
## metadata(1): colList
## assays(1): Abundance
## rownames(11724): 1 2 ... 11723 11724
## rowData names(97): Checked Protein.FDR.Confidence.Combined ...
## TMTpro.N.term.Positions Modifications
## colnames(53): F1.127N.Sample F1.127C.Sample ... F4.132C.Sample
## F4.133N.Sample
## colData names(0):
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):
Click to see which columns from the input file were used to form
each of the assays in the SingleCellExperiment object.
S4Vectors::metadata(sce)$colList## $Abundance
## [1] "Abundance.F1.127N.Sample" "Abundance.F1.127C.Sample"
## [3] "Abundance.F1.128N.Sample" "Abundance.F1.128C.Sample"
## [5] "Abundance.F1.129N.Sample" "Abundance.F1.129C.Sample"
## [7] "Abundance.F1.130N.Sample" "Abundance.F1.130C.Sample"
## [9] "Abundance.F1.131N.Sample" "Abundance.F1.131C.Sample"
## [11] "Abundance.F1.132N.Sample" "Abundance.F1.132C.Sample"
## [13] "Abundance.F1.133N.Sample" "Abundance.F2.127N.Sample"
## [15] "Abundance.F2.127C.Sample" "Abundance.F2.128N.Sample"
## [17] "Abundance.F2.128C.Sample" "Abundance.F2.129N.Sample"
## [19] "Abundance.F2.129C.Sample" "Abundance.F2.130N.Sample"
## [21] "Abundance.F2.130C.Sample" "Abundance.F2.131N.Sample"
## [23] "Abundance.F2.131C.Sample" "Abundance.F2.132N.Sample"
## [25] "Abundance.F2.132C.Sample" "Abundance.F2.133N.Sample"
## [27] "Abundance.F3.127N.Sample" "Abundance.F3.127C.Sample"
## [29] "Abundance.F3.128N.Sample" "Abundance.F3.128C.Sample"
## [31] "Abundance.F3.129N.Sample" "Abundance.F3.129C.Sample"
## [33] "Abundance.F3.130N.Sample" "Abundance.F3.130C.Sample"
## [35] "Abundance.F3.131N.Sample" "Abundance.F3.131C.Sample"
## [37] "Abundance.F3.132N.Sample" "Abundance.F3.132C.Sample"
## [39] "Abundance.F3.133N.Sample" "Abundance.F3.133C.Sample"
## [41] "Abundance.F4.127N.Sample" "Abundance.F4.127C.Sample"
## [43] "Abundance.F4.128N.Sample" "Abundance.F4.128C.Sample"
## [45] "Abundance.F4.129N.Sample" "Abundance.F4.129C.Sample"
## [47] "Abundance.F4.130N.Sample" "Abundance.F4.130C.Sample"
## [49] "Abundance.F4.131N.Sample" "Abundance.F4.131C.Sample"
## [51] "Abundance.F4.132N.Sample" "Abundance.F4.132C.Sample"
## [53] "Abundance.F4.133N.Sample"Add sample annotations
Next, we compile the sample annotations. The sample IDs have been
extracted from the column names in the ProteomeDiscoverer
file, by removing the provided iColPattern from the main
intensity columns. The group column will be used to define
groups for the statistical testing later. If a batch column
is present, that will also be accounted for in the limma
model. See the Comparisons and
design section below for more details about the fitted model(s).
Please check that the table below correspond to your expectations.
sce <- addSampleAnnots(sce, sampleAnnot = sampleAnnot)
DT::datatable(as.data.frame(colData(sce)),
options = list(scrollX = TRUE, pageLength = 20))Overview of the workflow
We already now define the names of the assays we will be generating and using later in the workflow.
aNames <- defineAssayNames(aName = aName, normMethod = normMethod,
doBatchCorr = "batch" %in% colnames(colData(sce)))
makeTableFromList(aNames)| assayInput | Abundance |
| assayLog2WithNA | log2_Abundance_withNA |
| assayImputIndic | imputed_Abundance |
| assayLog2NormWithNA | log2_Abundance_withNA_norm |
| assayImputed | log2_Abundance_norm |
| assayBatchCorr | log2_Abundance_norm_batchCorr |
The figure below provides a high-level overview of the workflow, and where the assays defined above will be generated. Note that depending on the settings specified by the user, not all steps may be performed (this is typically indicated by multiple assay names in the table above being equal).
knitr::include_graphics(system.file("extdata", "einprot_workflow.png",
package = "einprot"))
Overall distribution of feature intensities
The box plot below displays the distribution of values in the input assay defined above in each sample, on a log scale (excluding any missing values).
makeIntensityBoxplots(sce = sce, assayName = aNames$assayInput, doLog = TRUE,
ylab = aNames$assayInput)
Filter out contaminants and features with low confidence
Next, we filter out any proteins classified by PD as potential contaminants, and we may remove protein identifications based on score and number of peptides (i.e. to exclude one-hit wonders).
The excluded features, together with the available annotations, are written to a text file. The UpSet plot below illustrates the overlaps between the sets of features filtered out based on the different criteria (vertical bars).
nbrFeaturesBefore <- nrow(sce)
sce <- filterPDTMT(sce = sce, inputLevel = inputLevel, minScore = minScore,
minPeptides = minPeptides, minDeltaScore = minDeltaScore,
minPSMs = minPSMs, masterProteinsOnly = masterProteinsOnly,
modificationsCol = modificationsCol,
excludeUnmodifiedPeptides = excludeUnmodifiedPeptides,
keepModifications = keepModifications, plotUpset = TRUE,
exclFile = sub("\\.Rmd$", paste0("_filtered_out_features.txt"), knitr::current_input(dir = TRUE)))
nbrFeaturesAfter <- nrow(sce)
if (nbrFeaturesAfter == 0) {
stop("No features left after filtering!")
}This filtering removed 0 features. The new sce object
has 11724 features.
Modify feature names
The feature ID used when reading the data above are numeric indices.
We replace these IDs with more interpretable ones, corresponding to
idCol argument. We also add columns holding gene IDs (if
applicable), protein IDs, IDs for matching to STRING (if applicable) and
IDs for labelling points in plots.
sce <- fixFeatureIds(
sce,
colDefs = list(einprotId = idCol, einprotLabel = labelCol,
einprotGene = geneIdCol, einprotProtein = proteinIdCol,
IDsForSTRING = stringIdCol)
)
if (any(duplicated(rowData(sce)[["einprotId"]]))) {
stop("The 'einprotId' column cannot have duplicated entries.")
}
rownames(sce) <- rowData(sce)[["einprotId"]]Click to see examples of the defined feature identifiers
## einprotId einprotGene einprotProtein einprotLabel IDsForSTRING
## 1 TTN TTN Q8WZ42 TTN TTN
## 2 AHNAK AHNAK Q09666 AHNAK AHNAK
## 3 PLEC PLEC Q15149 PLEC PLEC
## 4 SPTAN1 SPTAN1 Q13813 SPTAN1 SPTAN1
## 5 MYH9 MYH9 P35579 MYH9 MYH9
## 6 ... ... ... ... ...
## 7 FOXL1 FOXL1 Q12952 FOXL1 FOXL1
## 8 TTC23 TTC23 Q5W5X9 TTC23 TTC23
## 9 GH2 GH2 P01242 GH2 GH2
## 10 CCN5 CCN5 O76076 CCN5 CCN5
## 11 Q14654 <NA> Q14654 Q14654 Q14654Prepare feature collections for later testing
In addition to testing individual features for differential abundance between groups, we can also test collections of features. Here we define the collections that will be used.
if (is.null(geneIdCol)) {
## If no gene IDs are extracted, don't compare to feature collections
featureCollections <- list()
} else {
(featureCollections <- prepareFeatureCollections(
sce = sce, idCol = "einprotGene",
includeFeatureCollections = includeFeatureCollections,
complexDbPath = complexDbPath, speciesInfo = speciesInfo,
complexSpecies = complexSpecies, customComplexes = customComplexes,
minSizeToKeep = minSizeToKeepSet))
}## $GO
## CharacterList of length 15200
## [["GOBP_10_FORMYLTETRAHYDROFOLATE_METABOLIC_PROCESS"]] ALDH1L1 ... MTHFD2L
## [["GOBP_2_OXOGLUTARATE_METABOLIC_PROCESS"]] OGDH OGDHL IDH2 ... AADAT TAT
## [["GOBP_2FE_2S_CLUSTER_ASSEMBLY"]] NFS1 GLRX3 GLRX5 HSCB BOLA2
## [["GOBP_3_PHOSPHOADENOSINE_5_PHOSPHOSULFATE_BIOSYNTHETIC_PROCESS"]] PAPSS2 ...
## [["GOBP_3_PHOSPHOADENOSINE_5_PHOSPHOSULFATE_METABOLIC_PROCESS"]] PAPSS2 ...
## [["GOBP_3_UTR_MEDIATED_MRNA_DESTABILIZATION"]] UPF1 KHSRP ... ZFP36 RBM24
## [["GOBP_3_UTR_MEDIATED_MRNA_STABILIZATION"]] HNRNPC RBM10 ... ELAVL4 RBM24
## [["GOBP_5_PHOSPHORIBOSE_1_DIPHOSPHATE_METABOLIC_PROCESS"]] PYGL ... PRPS1L1
## [["GOBP_5S_CLASS_RRNA_TRANSCRIPTION_BY_RNA_POLYMERASE_III"]] GTF3C1 ... GTF3C6
## [["GOBP_ABSCISSION"]] SPART VPS4A IST1 ZFYVE19 CHMP4C AURKB
## ...
## <15190 more elements>## Remove collections without any sets
featureCollections <- featureCollections[lapply(featureCollections, length) > 0]Apply a log2 transformation
Before the downstream analysis, we log2-transform the measured intensities. We also add an additional assay to keep track of the position of the missing values (which will be imputed later).
assay(sce, aNames$assayLog2WithNA) <- log2(assay(sce, aNames$assayInput))
## Add assay indicating missing values, which will be imputed
assay(sce, aNames$assayImputIndic) <- !is.finite(assay(sce, aNames$assayLog2WithNA))
sce## class: SingleCellExperiment
## dim: 11724 53
## metadata(1): colList
## assays(3): Abundance log2_Abundance_withNA imputed_Abundance
## rownames(11724): TTN AHNAK ... CCN5 Q14654
## rowData names(102): Checked Protein.FDR.Confidence.Combined ...
## einprotProtein IDsForSTRING
## colnames(53): F1.127N.Sample F1.127C.Sample ... F4.132C.Sample
## F4.133N.Sample
## colData names(3): sample group batch
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):Visualize missing value patterns
The plot below shows the fraction of the total set of features that are detected (with a non-missing value) in each of the samples.
## Replace zeros/-Inf values by explicit NA values in the assays
assay(sce, aNames$assayInput)[assay(sce, aNames$assayInput) == 0] <- NA
assay(sce, aNames$assayLog2WithNA)[!is.finite(assay(sce, aNames$assayLog2WithNA))] <- NA
## Count number of NA values and add to SCE
colData(sce)$nNA <- colSums(is.na(assay(sce, aNames$assayInput)))
colData(sce)$pNA <- 100 * sce$nNA/nrow(sce)
rowData(sce)$nNA <- rowSums(is.na(assay(sce, aNames$assayInput)))
rowData(sce)$pNA <- 100 * rowData(sce)$nNA/ncol(sce)
plotFractionDetectedPerSample(dfNA = as.data.frame(colData(sce)[, c("sample", "nNA", "pNA")]))
We also plot the number of features that are detected in a given number of samples.
plotDetectedInSamples(dfNA = as.data.frame(rowData(sce)[, c("nNA", "pNA")]))
The log2 intensities are next normalized across samples using the center.median method.
sce <- doNormalization(sce, method = normMethod,
assayName = aNames$assayLog2WithNA,
normalizedAssayName = aNames$assayLog2NormWithNA,
spikeFeatures = spikeFeatures)
makeIntensityBoxplots(sce = sce, assayName = aNames$assayLog2NormWithNA,
doLog = FALSE,
ylab = aNames$assayLog2NormWithNA)## Warning: Removed 99731 rows containing non-finite values (`stat_boxplot()`).
Imputation
Next, we apply the MinProb method to perform imputation of the log2-transformed data, and plot the distribution of imputed and non-imputed values in each sample.
set.seed(seed)
sce <- doImputation(sce, method = imputeMethod,
assayName = aNames$assayLog2NormWithNA,
imputedAssayName = aNames$assayImputed)## [1] 0.6858832plotImputationDistribution(sce, assayToPlot = aNames$assayImputed,
assayImputation = aNames$assayImputIndic,
xlab = aNames$assayImputed)
Overall distribution of log2 feature intensities
Next we consider the overall distribution of log2-intensities among the samples (after imputation).
makeIntensityBoxplots(sce = sce, assayName = aNames$assayImputed,
doLog = FALSE,
ylab = aNames$assayImputed)
if ("batch" %in% colnames(colData(sce))) {
if (subtractBaseline) {
assay(sce, aNames$assayBatchCorr) <-
getMatSubtractedBaseline(sce, assayName = aNames$assayImputed,
baselineGroup = baselineGroup,
sceFull = sce)
} else {
assay(sce, aNames$assayBatchCorr) <-
removeBatchEffect(assay(sce, aNames$assayImputed),
batch = sce$batch,
design = model.matrix(~ sce$group))
}
}Statistical testing
For each feature, we then compare the (possibly imputed) log2 intensities between groups. For this, we use the limma R/Bioconductor package (Ritchie et al. 2015; Phipson et al. 2016). For more information about the df.prior, representing the amount of extra information that is borrowed from the full set of features in order to improve the inference for each feature, see section 13.2 in the limma user guide. In addition to the feature-wise tests, we apply the camera method (Wu and Smyth 2012) to test for significance of each included feature collection. These tests are based on the t-statistics returned from limma.
Comparisons and design
## Set the assay to use for tests later
if (stattest == "proDA") {
assayForTests <- aNames$assayLog2NormWithNA
} else {
assayForTests <- aNames$assayImputed
}The log2_Abundance_norm assay will be used for the tests. The following pairwise comparisons will be performed (in each case, the first listed group will be the ‘baseline’ group):
if (subtractBaseline) {
discardGroup <- baselineGroup
} else {
discardGroup <- NULL
}
if (stattest == "none") {
comparisonList <- list(comparisons = list(),
groupComposition = list())
} else {
comparisonList <- makeListOfComparisons(
allGroups = unique(sce$group), comparisons = comparisons,
mergeGroups = mergeGroups,
allPairwiseComparisons = allPairwiseComparisons,
ctrlGroup = ctrlGroup, discardGroup = discardGroup)
}Click to expand list of comparisons
comparisonList$comparisons## $Medulla_COVID_vs_Medulla_Control
## [1] "Medulla_Control" "Medulla_COVID"
##
## $Cortex_COVID_vs_Cortex_Control
## [1] "Cortex_Control" "Cortex_COVID"Click to expand list of group compositions
comparisonList$groupComposition## $Cortex_COVID
## [1] "Cortex_COVID"
##
## $Medulla_COVID
## [1] "Medulla_COVID"
##
## $Medulla_Control
## [1] "Medulla_Control"
##
## $Cortex_Control
## [1] "Cortex_Control"if (any(assaysForExport %in% assayNames(sce))) {
assaysForExport <- intersect(assaysForExport, assayNames(sce))
} else {
assaysForExport <- aNames$assayInput
}testres <- runTest(sce = sce, comparisons = comparisonList$comparisons,
groupComposition = comparisonList$groupComposition,
testType = stattest,
assayForTests = assayForTests,
assayImputation = aNames$assayImputIndic,
minNbrValidValues = minNbrValidValues,
minlFC = minlFC, featureCollections = featureCollections,
complexFDRThr = complexFDRThr,
volcanoAdjPvalThr = volcanoAdjPvalThr,
volcanoLog2FCThr = volcanoLog2FCThr,
baseFileName = sub("\\.Rmd$", "", knitr::current_input()),
seed = seed, samSignificance = samSignificance,
nperm = nperm, volcanoS0 = volcanoS0,
aName = assaysForExport, addAbundanceValues = TRUE,
singleFit = singleFit,
subtractBaseline = subtractBaseline,
baselineGroup = baselineGroup,
extraColumns = union(interactiveDisplayColumns,
interactiveGroupColumn))
for (cmp in names(comparisonList$comparisons)) {
## Add WormBase/PomBase IDs if applicable
if (speciesInfo$speciesCommon == "roundworm") {
testres$tests[[cmp]]$WormBaseID <-
vapply(testres$tests[[cmp]][["einprotProtein"]],
function(mpds) {
wbids <- unlist(lapply(strsplit(mpds, ";")[[1]], function(mpd) {
wbconv$WormBaseID[wbconv$UniProtKB.ID == mpd |
wbconv$UniProtID == mpd]
}))
if (length(wbids[!is.na(wbids)]) != 0 &&
length(setdiff(wbids[!is.na(wbids)], "")) != 0) {
wbids <- setdiff(wbids, "")
wbids <- wbids[!is.na(wbids)]
paste(wbids, collapse = ";")
} else {
""
}
}, "NA")
} else if (speciesInfo$speciesCommon == "fission yeast") {
testres$tests[[cmp]]$PomBaseID <-
vapply(testres$tests[[cmp]][["einprotProtein"]],
function(mpds) {
pbids <- unlist(lapply(strsplit(mpds, ";")[[1]], function(mpd) {
pbconv$PomBaseID[pbconv$UniProtID == mpd]
}))
if (length(pbids[!is.na(pbids)]) != 0 &&
length(setdiff(pbids[!is.na(pbids)], "")) != 0) {
pbids <- setdiff(pbids, "")
pbids <- pbids[!is.na(pbids)]
paste(pbids, collapse = ";")
} else {
""
}
}, "NA")
}
}The plots below illustrate the experimental design used for the
linear model(s) and contrasts by limma. The plot to the
right shows the number of samples for each combination of factor levels
across the predictors, and is useful for detecting imbalances between
group sizes for different conditions. The plot to the left summarizes
the expected response value for each combination of predictor levels,
expressed in terms of the linear model coefficients. For more details on
how to interpret the plots, we refer to Soneson
et al. (2020) or Law et al. (2020).
Clicking on the arrow below the plots will reveal the design matrix used
by limma, as well as the contrasts that were fit for each
comparison.
if ("design" %in% names(testres$design)) {
cat("\n### Overall design \n")
vd <- VisualizeDesign(
testres$design$sampleData, designFormula = NULL,
designMatrix = testres$design$design)
print(cowplot::plot_grid(
plotlist = c(vd$plotlist, vd$cooccurrenceplots), nrow = 1)
)
cat("\n\n")
cat("<details>\n<summary><b>\nClick to display design matrix and contrast(s)\n</b></summary>\n")
cat("\n````\n")
print(testres$design$design)
cat("\n````\n")
cat("\n````\n")
cat("Contrast(s): \n")
print(testres$design$contrasts)
cat("\n````\n")
cat("\n````\n")
cat("Sample weights: \n")
print(testres$design$sampleWeights)
cat("\n````\n")
cat("\n</details>\n\n")
} else {
for (nm in names(testres$design)) {
cat("\n###", nm, "\n")
vd <- VisualizeDesign(
testres$design[[nm]]$sampleData, designFormula = NULL,
designMatrix = testres$design[[nm]]$design)
print(cowplot::plot_grid(
plotlist = c(vd$plotlist, vd$cooccurrenceplots), nrow = 1)
)
cat("\n\n")
cat("<details>\n<summary><b>\nClick to display design matrix and contrast(s)\n</b></summary>\n")
cat("\n````\n")
print(testres$design[[nm]]$design)
cat("\n````\n")
cat("\n````\n")
cat("Contrast: \n")
print(testres$design[[nm]]$contrast)
cat("\n````\n")
cat("\n````\n")
cat("Sample weights: \n")
print(testres$design[[nm]]$sampleWeights)
cat("\n````\n")
cat("\n</details>\n\n")
}
}Medulla_COVID_vs_Medulla_Control
Click to display design matrix and contrast(s)
(Intercept) bcb2 bcb3 bcb4 fcMedulla_COVID
F1.127C.Sample 1 0 0 0 1
F1.128C.Sample 1 0 0 0 1
F1.129C.Sample 1 0 0 0 1
F1.130C.Sample 1 0 0 0 1
F1.131C.Sample 1 0 0 0 0
F1.133N.Sample 1 0 0 0 0
F2.127C.Sample 1 1 0 0 1
F2.128C.Sample 1 1 0 0 1
F2.129C.Sample 1 1 0 0 1
F2.130C.Sample 1 1 0 0 1
F2.131C.Sample 1 1 0 0 0
F2.132C.Sample 1 1 0 0 0
F3.127C.Sample 1 0 1 0 1
F3.128C.Sample 1 0 1 0 1
F3.129C.Sample 1 0 1 0 1
F3.130C.Sample 1 0 1 0 1
F3.131C.Sample 1 0 1 0 1
F3.132C.Sample 1 0 1 0 0
F3.133C.Sample 1 0 1 0 0
F4.127C.Sample 1 0 0 1 1
F4.128C.Sample 1 0 0 1 1
F4.130N.Sample 1 0 0 1 1
F4.131N.Sample 1 0 0 1 0
F4.132N.Sample 1 0 0 1 0
F4.133N.Sample 1 0 0 1 0
attr(,"assign")
[1] 0 1 1 1 2
attr(,"contrasts")
attr(,"contrasts")$bc
[1] "contr.treatment"
attr(,"contrasts")$fc
[1] "contr.treatment"
Contrast:
[1] 0 0 0 0 1
Sample weights:
NULL
Cortex_COVID_vs_Cortex_Control
Click to display design matrix and contrast(s)
(Intercept) bcb2 bcb3 bcb4 fcCortex_COVID
F1.127N.Sample 1 0 0 0 1
F1.128N.Sample 1 0 0 0 1
F1.129N.Sample 1 0 0 0 1
F1.130N.Sample 1 0 0 0 1
F1.131N.Sample 1 0 0 0 1
F1.132N.Sample 1 0 0 0 0
F1.132C.Sample 1 0 0 0 0
F2.127N.Sample 1 1 0 0 1
F2.128N.Sample 1 1 0 0 1
F2.129N.Sample 1 1 0 0 1
F2.130N.Sample 1 1 0 0 1
F2.131N.Sample 1 1 0 0 0
F2.132N.Sample 1 1 0 0 0
F2.133N.Sample 1 1 0 0 0
F3.127N.Sample 1 0 1 0 1
F3.128N.Sample 1 0 1 0 1
F3.129N.Sample 1 0 1 0 1
F3.130N.Sample 1 0 1 0 1
F3.131N.Sample 1 0 1 0 1
F3.132N.Sample 1 0 1 0 0
F3.133N.Sample 1 0 1 0 0
F4.127N.Sample 1 0 0 1 1
F4.128N.Sample 1 0 0 1 1
F4.129N.Sample 1 0 0 1 1
F4.129C.Sample 1 0 0 1 1
F4.130C.Sample 1 0 0 1 0
F4.131C.Sample 1 0 0 1 0
F4.132C.Sample 1 0 0 1 0
attr(,"assign")
[1] 0 1 1 1 2
attr(,"contrasts")
attr(,"contrasts")$bc
[1] "contr.treatment"
attr(,"contrasts")$fc
[1] "contr.treatment"
Contrast:
[1] 0 0 0 0 1
Sample weights:
NULL
SA plots
We first show a diagnostic plot for each comparison. These plots
display the square root of the residual standard deviation (y-axis)
versus the mean abundance (across all the groups used to perform the
model fit, x-axis). The curve indicated in the plots show the
mean-variance trend inferred by limma.
for (i in seq.int(ceiling(length(testres$tests) / 3))) {
tmplist <- testres$tests[(seq_along(testres$tests) - 1) %/% 3 == (i - 1)]
print(makeSAPlot(tmplist))
}
Volcano plots
Below we display a volcano plot for each comparison. These plots are also saved to pdf files. In each plot, a subset of (up to 25) significant hits are indicated by name (selected as the ones with the largest Manhattan distance to the origin). These proteins are also used to generate STRING networks (Szklarczyk et al. 2021) (separately for the up- and downregulated ones), which are included in the pdf file. Any features explicitly requested (see the table above) are also labeled in the volcano plots. In addition to these pdf files, if “complexes” is specified to be included in the feature collections (and tested for significance using camera), we also generate a multi-page pdf file showing the position of the features of each significantly differentially abundant complex in the volcano plot, as well as bar plots of the features’ abundance values in the compared samples. This pdf file is only generated if there is at least one significant complex (with adjusted p-value below the specified complexFDRThr=0.1).
interactiveVolcanos <- htmltools::tagList()
for (nm in names(testres$tests)) {
plots <- plotVolcano(
sce = sce, res = testres$tests[[nm]], testType = stattest,
xv = NULL, yv = NULL, xvma = NULL, volcind = NULL,
plotnote = testres$plotnotes[[nm]],
plottitle = testres$plottitles[[nm]],
plotsubtitle = testres$plotsubtitles[[nm]],
volcanoFeaturesToLabel = volcanoFeaturesToLabel,
volcanoMaxFeatures = volcanoMaxFeatures,
volcanoLabelSign = volcanoLabelSign,
baseFileName = paste0(sub("\\.Rmd$", "", knitr::current_input()),
"_volcano_", nm),
comparisonString = nm,
groupComposition = comparisonList$groupComposition[comparisonList$comparisons[[nm]]],
stringDb = string_db,
featureCollections = testres$featureCollections,
complexFDRThr = complexFDRThr,
maxNbrComplexesToPlot = maxNbrComplexesToPlot,
curveparam = testres$curveparams[[nm]],
abundanceColPat = assaysForExport,
xlab = "log2(fold change)",
ylab = "-log10(p-value)",
xlabma = "Average abundance",
labelOnlySignificant = TRUE,
interactiveDisplayColumns = interactiveDisplayColumns,
interactiveGroupColumn = interactiveGroupColumn,
maxTextWidthBarplot = 5.1)
if (!is.null(plots$ggma) && !is.null(plots$ggwf) && !is.null(plots$ggbar)) {
cat("\n\n### ", nm, " \n\n\n")
print(plots$gg)
print(plots$ggma)
for (ggb in plots$ggbar) print(ggb)
print(plots$ggwf)
cat("\n\n")
} else if (!is.null(plots$ggma) && !is.null(plots$ggwf)) {
cat("\n\n### ", nm, " \n\n\n")
print(plots$gg)
print(plots$ggma)
print(plots$ggwf)
cat("\n\n")
} else if (!is.null(plots$ggma) && !is.null(plots$ggbar)) {
cat("\n\n### ", nm, " \n\n\n")
print(plots$gg)
print(plots$ggma)
for (ggb in plots$ggbar) print(ggb)
cat("\n\n")
} else if (!is.null(plots$ggbar) && !is.null(plots$ggwf)) {
cat("\n\n### ", nm, " \n\n\n")
print(plots$gg)
for (ggb in plots$ggbar) print(ggb)
print(plots$ggwf)
cat("\n\n")
} else if (!is.null(plots$ggma)) {
cat("\n\n### ", nm, " \n\n\n")
print(plots$gg)
print(plots$ggma)
cat("\n\n")
} else if (!is.null(plots$ggbar)) {
cat("\n\n### ", nm, " \n\n\n")
print(plots$gg)
for (ggb in plots$ggbar) print(ggb)
cat("\n\n")
} else if (!is.null(plots$ggwf)) {
cat("\n\n### ", nm, " \n\n\n")
print(plots$gg)
print(plots$ggwf)
cat("\n\n")
} else {
cat("\n\n### ", nm, " \n\n\n")
print(plots$gg)
cat("\n\n")
}
interactiveVolcanos[[nm]] <- plots$ggint
}Medulla_COVID_vs_Medulla_Control
Warning: we couldn’t map to STRING 4% of your
identifiers
Cortex_COVID_vs_Cortex_Control
