Ostapcuk et al 2018
MaxQuant/LFQ analysis
2023-07-19 20:41:05 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 <-
"mouse"
mqFile <-
"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/einprot/extdata/mq_example/1356_proteinGroups.txt"
mqParameterFile <-
"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/einprot/extdata/mq_example/1356_mqpar.xml"
idCol <-
function (df)
combineIds(df, combineCols = c("Gene.names", "Majority.protein.IDs"))
labelCol <-
function (df)
combineIds(df, combineCols = c("Gene.names", "Majority.protein.IDs"))
geneIdCol <-
function (df)
getFirstId(df, colName = "Gene.names")
proteinIdCol <-
"Majority.protein.IDs"
stringIdCol <-
function (df)
combineIds(df, combineCols = c("Gene.names", "Majority.protein.IDs"),
combineWhen = "missing", makeUnique = FALSE)
reportTitle <-
"Ostapcuk et al 2018"
reportAuthor <-
"MaxQuant/LFQ analysis"
iColPattern <-
"^LFQ.intensity."
sampleAnnot <-
structure(list(sample = c("Adnp_IP04", "Adnp_IP05", "Adnp_IP06",
"Chd4BF_IP07", "Chd4BF_IP08", "Chd4BF_IP09", "RBC_ctrl_IP01",
"RBC_ctrl_IP02", "RBC_ctrl_IP03"), group = c("Adnp", "Adnp",
"Adnp", "Chd4BF", "Chd4BF", "Chd4BF", "RBC_ctrl", "RBC_ctrl",
"RBC_ctrl")), class = "data.frame", row.names = c(NA, -9L))
includeOnlySamples <-
""
excludeSamples <-
""
minScore <-
10
minPeptides <-
2
imputeMethod <-
"MinProb"
assaysForExport <-
c("iBAQ", "Top3")
mergeGroups <-
list()
comparisons <-
list(c("RBC_ctrl", "Adnp"), c("RBC_ctrl", "Chd4BF"))
ctrlGroup <-
""
allPairwiseComparisons <-
TRUE
singleFit <-
FALSE
subtractBaseline <-
FALSE
baselineGroup <-
""
normMethod <-
"none"
spikeFeatures <-
NULL
stattest <-
"limma"
minNbrValidValues <-
2
minlFC <-
0
samSignificance <-
TRUE
nperm <-
250
volcanoAdjPvalThr <-
0.050000000000000003
volcanoLog2FCThr <-
1
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 <-
Inf
seed <-
42
includeFeatureCollections <-
"complexes"
minSizeToKeepSet <-
2
customComplexes <-
list()
complexSpecies <-
"all"
complexDbPath <-
NULL
stringVersion <-
"11.5"
stringDir <-
NULL
linkTableColumns <-
c("^Unique.peptides", "[^se]\\.logFC$")This report describes a reproducible end-to-end analysis of a
proteomics dataset quantified with MaxQuant (Cox and Mann 2008). 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 additional global overviews of the data are provided via
principal component analysis and summary heatmaps.
## Get species info and define STRINGdb object
speciesInfo <- getSpeciesInfo(species)
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
})
}## *********** STRING - https://string-db.org ***********
## (Search Tool for the Retrieval of Interacting Genes/Proteins)
## version: 11.5
## species: 10090
## ............please wait............
## proteins: 22048
## interactions: 1851626## 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 | Mus musculus |
| Species (common) | mouse |
| Taxonomic ID | 10090 |
MaxQuant analysis summary
mq <- readMaxQuantXML(mqParameterFile)
mq <- c(list("MaxQuant file" = mqFile), mq)
makeTableFromList(mq)| MaxQuant file | /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/einprot/extdata/mq_example/1356_proteinGroups.txt |
| MaxQuant version | 1.5.3.8 |
| Parameter file | /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/einprot/extdata/mq_example/1356_mqpar.xml |
| Search engine | Andromeda |
| Raw file location | D:/Data/MaxQuant/ |
| Raw files | F_160817_AdnpFB_IP04.raw, F_160817_AdnpFB_IP05.raw, F_160817_AdnpFB_IP06.raw, F_160817_Chd4BF_IP07.raw, F_160817_Chd4BF_IP08.raw, F_160817_Chd4BF_IP09.raw, F_160817_RBC_ctrl_IP01.raw, F_160817_RBC_ctrl_IP02.raw, F_160817_RBC_ctrl_IP03.raw |
| Sample names | Adnp_IP04, Adnp_IP05, Adnp_IP06, Chd4BF_IP07, Chd4BF_IP08, Chd4BF_IP09, RBC_ctrl_IP01, RBC_ctrl_IP02, RBC_ctrl_IP03 |
| Databases | D:/Databases/MOUSE__150129.fasta |
| Contaminants | */MaxQuant_1.5.3.8/MaxQuant/bin/conf/contaminants.fasta |
| Quantification settings (LFQ) | LFQ min. ratio count: 1, fastLFQ: false, match-between runs (MBR): true, Intensity based absolute quantification (iBAQ): true |
| Min. razor peptides | 1 |
| Enzymes | Trypsin/P |
| Variable modifications | Oxidation (M), Acetyl (Protein N-term) |
| Fixed modifications | |
| Max peptide mass | 8000 |
| Min peptide length | 7 |
Settings
settingsList <- list(
"Column pattern" = iColPattern,
"Include only samples (if applicable)" = paste(includeOnlySamples,
collapse = ", "),
"Exclude samples" = paste(excludeSamples, collapse = ", "),
"Min. number of peptides" = minPeptides,
"Min. protein score" = minScore,
"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,
"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"))]
}
makeTableFromList(settingsList)| Column pattern | ^LFQ.intensity. |
| Include only samples (if applicable) | |
| Exclude samples | |
| Min. number of peptides | 2 |
| Min. protein score | 10 |
| Imputation method | MinProb |
| Assays(s) to use for exported values | iBAQ, Top3 |
| Min. nbr valid values required for testing | 2 |
| Model fit | Separate model fit for each comparison |
| Groups to merge | |
| Comparisons | RBC_ctrl vs Adnp; RBC_ctrl vs Chd4BF |
| Control group | |
| Do all pairwise comparisons | TRUE |
| Batch correction via baseline subtraction | FALSE |
| Baseline group | |
| Normalization method | none |
| 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 | 1 |
| 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 | complexes |
| Minimal required size for feature sets | 2 |
| Complexes file | complexdb_einprot0.5.0_20220323_orthologs.rds |
| Complexes from species | all |
| Custom complexes | |
| FDR Threshold for complexes | 0.1 |
| Max nbr complexes to plot | Inf |
| Random seed | 42 |
| Columns to add in link table | ^Unique.peptides;[^se].logFC$ |
| Interactive display columns | einprotLabel;adj.P.Val;logFC |
Read MaxQuant output
The input to this workflow is a proteinGroups.txt file
from MaxQuant (see path in the table above). We read the MaxQuant
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 MaxQuant file.
## Read MaxQuant output
tmp <- importExperiment(inFile = mqFile, iColPattern = iColPattern,
includeOnlySamples = includeOnlySamples,
excludeSamples = excludeSamples)
sce <- tmp$sce
aName <- tmp$aName
sce## class: SingleCellExperiment
## dim: 463 9
## metadata(1): colList
## assays(9): LFQ.intensity MS.MS.Count ... iBAQ Identification.type
## rownames(463): 1 2 ... 462 463
## rowData names(35): Protein.IDs Majority.protein.IDs ...
## Oxidation.M.site.IDs Oxidation.M.site.positions
## colnames(9): Adnp_IP04 Adnp_IP05 ... RBC_ctrl_IP02 RBC_ctrl_IP03
## 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## $LFQ.intensity
## [1] "LFQ.intensity.Adnp_IP04" "LFQ.intensity.Adnp_IP05"
## [3] "LFQ.intensity.Adnp_IP06" "LFQ.intensity.Chd4BF_IP07"
## [5] "LFQ.intensity.Chd4BF_IP08" "LFQ.intensity.Chd4BF_IP09"
## [7] "LFQ.intensity.RBC_ctrl_IP01" "LFQ.intensity.RBC_ctrl_IP02"
## [9] "LFQ.intensity.RBC_ctrl_IP03"
##
## $MS.MS.Count
## [1] "MS.MS.Count.Adnp_IP04" "MS.MS.Count.Adnp_IP05"
## [3] "MS.MS.Count.Adnp_IP06" "MS.MS.Count.Chd4BF_IP07"
## [5] "MS.MS.Count.Chd4BF_IP08" "MS.MS.Count.Chd4BF_IP09"
## [7] "MS.MS.Count.RBC_ctrl_IP01" "MS.MS.Count.RBC_ctrl_IP02"
## [9] "MS.MS.Count.RBC_ctrl_IP03"
##
## $Intensity
## [1] "Intensity.Adnp_IP04" "Intensity.Adnp_IP05"
## [3] "Intensity.Adnp_IP06" "Intensity.Chd4BF_IP07"
## [5] "Intensity.Chd4BF_IP08" "Intensity.Chd4BF_IP09"
## [7] "Intensity.RBC_ctrl_IP01" "Intensity.RBC_ctrl_IP02"
## [9] "Intensity.RBC_ctrl_IP03"
##
## $Sequence.coverage
## [1] "Sequence.coverage.Adnp_IP04...." "Sequence.coverage.Adnp_IP05...."
## [3] "Sequence.coverage.Adnp_IP06...." "Sequence.coverage.Chd4BF_IP07...."
## [5] "Sequence.coverage.Chd4BF_IP08...." "Sequence.coverage.Chd4BF_IP09...."
## [7] "Sequence.coverage.RBC_ctrl_IP01...." "Sequence.coverage.RBC_ctrl_IP02...."
## [9] "Sequence.coverage.RBC_ctrl_IP03...."
##
## $Unique.peptides
## [1] "Unique.peptides.Adnp_IP04" "Unique.peptides.Adnp_IP05"
## [3] "Unique.peptides.Adnp_IP06" "Unique.peptides.Chd4BF_IP07"
## [5] "Unique.peptides.Chd4BF_IP08" "Unique.peptides.Chd4BF_IP09"
## [7] "Unique.peptides.RBC_ctrl_IP01" "Unique.peptides.RBC_ctrl_IP02"
## [9] "Unique.peptides.RBC_ctrl_IP03"
##
## $Razor.unique.peptides
## [1] "Razor...unique.peptides.Adnp_IP04"
## [2] "Razor...unique.peptides.Adnp_IP05"
## [3] "Razor...unique.peptides.Adnp_IP06"
## [4] "Razor...unique.peptides.Chd4BF_IP07"
## [5] "Razor...unique.peptides.Chd4BF_IP08"
## [6] "Razor...unique.peptides.Chd4BF_IP09"
## [7] "Razor...unique.peptides.RBC_ctrl_IP01"
## [8] "Razor...unique.peptides.RBC_ctrl_IP02"
## [9] "Razor...unique.peptides.RBC_ctrl_IP03"
##
## $Peptides
## [1] "Peptides.Adnp_IP04" "Peptides.Adnp_IP05" "Peptides.Adnp_IP06"
## [4] "Peptides.Chd4BF_IP07" "Peptides.Chd4BF_IP08" "Peptides.Chd4BF_IP09"
## [7] "Peptides.RBC_ctrl_IP01" "Peptides.RBC_ctrl_IP02" "Peptides.RBC_ctrl_IP03"
##
## $iBAQ
## [1] "iBAQ.Adnp_IP04" "iBAQ.Adnp_IP05" "iBAQ.Adnp_IP06"
## [4] "iBAQ.Chd4BF_IP07" "iBAQ.Chd4BF_IP08" "iBAQ.Chd4BF_IP09"
## [7] "iBAQ.RBC_ctrl_IP01" "iBAQ.RBC_ctrl_IP02" "iBAQ.RBC_ctrl_IP03"
##
## $Identification.type
## [1] "Identification.type.Adnp_IP04" "Identification.type.Adnp_IP05"
## [3] "Identification.type.Adnp_IP06" "Identification.type.Chd4BF_IP07"
## [5] "Identification.type.Chd4BF_IP08" "Identification.type.Chd4BF_IP09"
## [7] "Identification.type.RBC_ctrl_IP01" "Identification.type.RBC_ctrl_IP02"
## [9] "Identification.type.RBC_ctrl_IP03"Add sample annotations
Next, we compile the sample annotations. The sample IDs have been
extracted from the column names in the MaxQuant 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 | LFQ.intensity |
| assayLog2WithNA | log2_LFQ.intensity_withNA |
| assayImputIndic | imputed_LFQ.intensity |
| assayLog2NormWithNA | log2_LFQ.intensity_withNA |
| assayImputed | log2_LFQ.intensity |
| assayBatchCorr | log2_LFQ.intensity |
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 protein intensities
The box plot below displays the distribution of raw protein intensities in each sample, on a log scale (excluding any missing values), as well as the values in the input assay defined above.
if ("Intensity" %in% assayNames(sce)) {
print(makeIntensityBoxplots(sce = sce, assayName = "Intensity",
doLog = TRUE, ylab = "Intensity"))
}
makeIntensityBoxplots(sce = sce, assayName = aNames$assayInput, doLog = TRUE,
ylab = aNames$assayInput)
Filter out contaminants, reverse hits, and proteins with low confidence
Next, we filter out any features classified by MaxQuant as potential contaminants or reverse (decoy) hits, and features identified only by site (which removes proteins that are only identified by peptides carrying one or more modified amino acids). In addition, 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 <- filterMaxQuant(sce = sce, minScore = minScore, minPeptides = minPeptides,
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 232 features. The new sce object
has 231 features.
Modify feature names
The feature ID used when reading the data above are numeric indices.
We replace these IDs with more interpretable ones, defined by the
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
## 1 Dhx9 Dhx9 A0A087WPL5;E9QNN1;O70133;O70133-2;O70133-3;Q3UR42
## 2 Zmynd8 Zmynd8 Q3UH28;Q3U1M7;A2A483;E9Q8D1;A2A482;A2A484;A2A485
## 3 Zmym4 Zmym4 A2A791;A2A791-2;F6VYE2
## 4 Rlf Rlf A2A7F4;E9Q532
## 5 Zfp600 Zfp600 A2A7V0
## 6 ... ... ...
## 7 Mbd3.Q9Z2D8 Mbd3 Q9Z2D8;D3YTR5;Q9Z2D8-2;F6WZH1
## 8 Mbd2 Mbd2 Q9Z2E1
## 9 Zfp292 Zfp292 Q9Z2U2-2;Q9Z2U2
## 10 Hnrnpf Hnrnpf Q9Z2X1-2;Q9Z2X1;J3QM80;J3QMT0
## 11 Zfp280d Zfp280d V9GWW1;Q68FE8;Q68FE8-2;Q68FE8-3
## einprotLabel IDsForSTRING
## 1 Dhx9 Dhx9
## 2 Zmynd8 Zmynd8
## 3 Zmym4 Zmym4
## 4 Rlf Rlf
## 5 Zfp600 Zfp600
## 6 ... ...
## 7 Mbd3.Q9Z2D8 Mbd3
## 8 Mbd2 Mbd2
## 9 Zfp292 Zfp292
## 10 Hnrnpf Hnrnpf
## 11 Zfp280d Zfp280dPrepare feature collections for later testing
In addition to testing individual proteins for differential abundance between groups, we can also test collections of proteins. 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))
}## $complexes
## CharacterList of length 151
## [["mouse: DCS complex (Ptbp1, Ptbp2, Hnrph1, Hnrpf) (+1 alt. ID)"]] Hnrnph1 ...
## [["mouse: Drosha complex"]] Hnrnph1 Ddx5 Hnrnpm
## [["mouse: Gata1-Fog1-MeCP1 complex"]] Mbd3.D3YTR4 Mta1.E9PX23 ... Mbd2
## [["mouse: Hdac1-Ino80-Kdm1a-Phb2-Rbp1-Taf5 complex"]] Kdm1a Hdac1
## [["mouse: Ikaros complex"]] Chd4 Hdac1 Hdac2 Rbbp4
## [["mouse: Ikaros-NuRD complex"]] Chd4 Hdac1 Hdac2
## [["mouse: Metallothionein-3 complex"]] Hsp90ab1 Actb
## [["mouse: Nkx3.2-SMAD1-SMAD4-HDAC-Sin3A complex"]] Hdac1 Rbbp4 Rbbp7
## [["mouse: Nucleolar remodeling complex (NoRC complex) (+1 alt. ID)"]] Baz2a ...
## [["mouse: Parvulin-associated pre-rRNP complex"]] Ncl Rpl7a Rpl7 ... Nop56 Rpl4
## ...
## <141 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: 231 9
## metadata(1): colList
## assays(11): LFQ.intensity MS.MS.Count ... log2_LFQ.intensity_withNA
## imputed_LFQ.intensity
## rownames(231): Dhx9 Zmynd8 ... Hnrnpf Zfp280d
## rowData names(40): Protein.IDs Majority.protein.IDs ... einprotProtein
## IDsForSTRING
## colnames(9): Adnp_IP04 Adnp_IP05 ... RBC_ctrl_IP02 RBC_ctrl_IP03
## colData names(2): sample group
## 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")]))
Finally we show the full missing value structure in the data matrix. The white values correspond to missing values, and the grey ones to observed intensities.
plotMissingValuesHeatmap(sce, assayMissing = aNames$assayImputIndic)
The log2 intensities are not normalized further across samples since ‘normMethod’ is set to ‘none’.
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.9660524plotImputationDistribution(sce, assayToPlot = aNames$assayImputed,
assayImputation = aNames$assayImputIndic,
xlab = aNames$assayImputed)
Overall distribution of log2 protein 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_LFQ.intensity 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## $Adnp_vs_RBC_ctrl
## [1] "RBC_ctrl" "Adnp"
##
## $Chd4BF_vs_RBC_ctrl
## [1] "RBC_ctrl" "Chd4BF"Click to expand list of group compositions
comparisonList$groupComposition## $Adnp
## [1] "Adnp"
##
## $Chd4BF
## [1] "Chd4BF"
##
## $RBC_ctrl
## [1] "RBC_ctrl"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")
}
}Adnp_vs_RBC_ctrl
Click to display design matrix and contrast(s)
(Intercept) fcAdnp
Adnp_IP04 1 1
Adnp_IP05 1 1
Adnp_IP06 1 1
RBC_ctrl_IP01 1 0
RBC_ctrl_IP02 1 0
RBC_ctrl_IP03 1 0
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$fc
[1] "contr.treatment"
Contrast:
[1] 0 1
Sample weights:
NULL
Chd4BF_vs_RBC_ctrl
Click to display design matrix and contrast(s)
(Intercept) fcChd4BF
Chd4BF_IP07 1 1
Chd4BF_IP08 1 1
Chd4BF_IP09 1 1
RBC_ctrl_IP01 1 0
RBC_ctrl_IP02 1 0
RBC_ctrl_IP03 1 0
attr(,"assign")
[1] 0 1
attr(,"contrasts")
attr(,"contrasts")$fc
[1] "contr.treatment"
Contrast:
[1] 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 proteins of each significantly differentially abundant complex in the volcano plot, as well as bar plots of the proteins’ 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
}Adnp_vs_RBC_ctrl
Chd4BF_vs_RBC_ctrl
Interactive volcano plots
Result export
For each comparison, we also save a text file with the “significant” features (defined as those colored in the volcano plots above). The features are ordered by the logFC value. In addition, we save a file with abundance values for all features that are “significant” in at least one comparison, across all the samples used in at least one comparison.
abundanceExport <- makeAbundanceExport(testresList = testres$tests,
abundancePrefix = assaysForExport)
write.table(abundanceExport, file = sub("\\.Rmd$", "_abundance_values_significant.txt",
knitr::current_input()),
row.names = FALSE, col.names = TRUE, quote = FALSE, sep = "\t")## Merge results from all tests and add to rowData(sce)
tests <- testres$tests
for (nm in names(tests)) {
idx <- which(colnames(tests[[nm]]) != "pid")
colnames(tests[[nm]])[idx] <- paste0(nm, ".", colnames(tests[[nm]])[idx])
}
all_tests <- as.data.frame(Reduce(function(...) dplyr::full_join(..., by = "pid"),
tests), optional = TRUE)
rownames(all_tests) <- all_tests$pid
all_tests$pid <- NULL
stopifnot(rownames(sce) == rownames(all_tests))
rowData(sce) <- cbind(rowData(sce), all_tests)Overlap among sets of significant features
The UpSet plot below shows the overlap among the “significant” proteins (defined as the ones that are colored in the volcano plots above, based on the user specifications) from the different comparisons. Note that if there are many comparisons, not all combinations may be displayed in the plot (only the 50 combinations with the largest number of proteins are shown, for interpretability reasons). Moreover, the UpSet plot will only be shown if at least two comparisons have been made, and there are at least two comparisons where any proteins were deemed significant.
tmpsign <- all_tests %>% dplyr::select(contains("showInVolcano")) + 0
tmpsign[is.na(tmpsign)] <- 0
## Only features that are significant in at least one comparison
tmpsign <- tmpsign[rowSums(tmpsign) > 0, , drop = FALSE]
colnames(tmpsign) <- gsub("\\.showInVolcano$", "", colnames(tmpsign))
if (length(tests) > 1 && sum(colSums(tmpsign) > 0) > 1) {
ComplexUpset::upset(tmpsign, intersect = colnames(tmpsign),
sort_intersections_by = "cardinality")
}
Most significant feature sets
Finally, we display the top significant feature sets in each of the tested collections, for each comparison (if applicable). We recommend that the (adjusted) p-values for the feature sets is interpreted with caution, especially in situations where the protein abundances are measured on an isoform level and the feature sets are defined on the protein level, since there will be many (sometimes strongly correlated) features corresponding to a single gene or protein annotated to a feature set.
for (nm in names(testres$topsets)) {
if (length(testres$topsets[[nm]]) > 0) {
plts <- lapply(names(testres$topsets[[nm]]), function(snm) {
df <- testres$topsets[[nm]][[snm]]
if (nrow(df) > maxNbrComplexesToPlot) {
df <- df[seq_len(maxNbrComplexesToPlot), ]
}
if (nrow(df) > 0) {
ggplot(df %>% dplyr::mutate(set = factor(.data$set, levels = rev(.data$set))),
aes(x = .data$set, y = -log10(.data[[paste0(nm, "_FDR")]]),
fill = .data[[paste0(nm, "_Direction")]])) +
geom_bar(stat = "identity") +
coord_flip() + theme_bw() +
labs(x = "", title = snm) +
scale_fill_manual(values = c(Up = scales::muted("blue"),
Down = scales::muted("red")),
name = "Direction") +
theme(axis.text = element_text(size = 6))
} else {
NULL
}
})
if (sum(sapply(plts, is.null)) < length(plts)) {
cat("\n\n### ", nm, " \n\n\n")
print(cowplot::plot_grid(plotlist = plts, ncol = 1, align = "v"))
cat("\n\n")
}
}
}Adnp_vs_RBC_ctrl
Chd4BF_vs_RBC_ctrl
Table with direct database links to sequences, functional information and predicted structures
The table below provides autogenerated links to the UniProt and
AlphaFold pages (as well as selected organism-specific databases) for
the majority protein IDs corresponding to each feature in the data set.
The ‘pid’ column represents the unique feature ID used by
einprot, and the einprotLabel column contains
the user-defined feature labels. UniProt is a resource of protein
sequence and functional information hosted by EMBL-EBI, PIR and SIB. The
AlphaFold Protein Structure Database, developed by DeepMind and
EMBL-EBI, provides open access to protein structure predictions for the
human proteome and other key proteins of interest. Note that (depending
on the species) many proteins are not yet covered in AlphaFold (in this
case, the link below will lead to a non-existent page), and that numeric
values are rounded to four significant digits to increase
readability.
linkTable <- makeDbLinkTable(
df = as.data.frame(rowData(sce)) %>%
rownames_to_column("pid") %>%
dplyr::select("pid", "einprotProtein",
matches(setdiff(c("^einprotLabel$", linkTableColumns), ""), perl = TRUE)),
idCol = "einprotProtein",
speciesCommon = speciesInfo$speciesCommon,
addSpeciesSpecificColumns = TRUE,
convTablePomBase = pbconv,
convTableWormBase = wbconv,
removeSuffix = TRUE, signifDigits = 4
) %>%
dplyr::select(-.data$einprotProtein)
DT::datatable(
as.data.frame(linkTable), escape = FALSE,
filter = list(position = "top", clear = FALSE),
extensions = "Buttons",
options = list(scrollX = TRUE, pageLength = 20,
search = list(regex = FALSE, caseInsensitive = TRUE),
dom = "Bfrltip", buttons =
list(list(extend = "csv",
filename = paste0(sub("\\.Rmd$", "",
knitr::current_input()), "_linktable")),
list(extend = "excel", title = "",
filename = paste0(sub("\\.Rmd$", "",
knitr::current_input()), "_linktable"))))
)Assemble SingleCellExperiment object
We assemble all the information calculated above in a SingleCellExperiment object, which can later be used e.g. for exploration with iSEE (Rue-Albrecht et al. 2018).
sce <- prepareFinalSCE(
sce = sce, baseFileName = sub("\\.Rmd$", "", knitr::current_input()),
featureCollections = testres$featureCollections, expType = "MaxQuant")
## Add experiment metadata
S4Vectors::metadata(sce) <- c(
S4Vectors::metadata(sce),
list(
mqFile = mqFile,
aNames = aNames,
aName = aNames$assayInput,
iColPattern = iColPattern,
imputeMethod = imputeMethod,
ctrlGroup = ctrlGroup,
allPairwiseComparisons = allPairwiseComparisons,
normMethod = normMethod,
stattest = stattest,
minlFC = minlFC,
analysisDate = as.character(Sys.Date()),
rmdFile = knitr::current_input(dir = TRUE),
testres = testres,
comparisonList = comparisonList
)
)
sce## class: SingleCellExperiment
## dim: 231 9
## metadata(16): colList iSEE ... testres comparisonList
## assays(13): LFQ.intensity MS.MS.Count ... log2_LFQ.intensity
## log2_iBAQ_withNA
## rownames(231): Dhx9 Zmynd8 ... Hnrnpf Zfp280d
## rowData names(97): Protein.IDs Majority.protein.IDs ...
## Chd4BF_vs_RBC_ctrl.log2_iBAQ.RBC_ctrl.sd
## Chd4BF_vs_RBC_ctrl.IDsForSTRING
## colnames(9): Adnp_IP04 Adnp_IP05 ... RBC_ctrl_IP02 RBC_ctrl_IP03
## colData names(4): sample group nNA pNA
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):Run PCA
pcafeatures <- which(rowSums(!assay(sce, aNames$assayImputIndic)) >= minNbrValidValues)We run a principal component analysis to obtain a reduced dimensionality representation of the data, in order to visualize the samples in two dimensions. The PCA is based on features with at least 2 non-imputed values (222 of the 231 features). The figure below shows the sample representation in the first two principal components, the fraction of variance explained by each of the principal components, and the features with the highest positive and negative loadings in the displayed components.
interactivePCAs <- htmltools::tagList()
for (a in unique(c(aNames$assayImputed, aNames$assayBatchCorr))) {
pcares <- doPCA(sce = sce, assayName = a, ncomponents = 10, ntop = Inf,
plotpairs = list(c(1, 2)), maxTextWidthBarplot = 1.9,
subset_row = pcafeatures)
sce <- pcares$sce
interactivePCAs[[a]] <- ggplotly(pcares$plotcoord[[1]],
width = 750, height = 500)
for (plc in pcares$plotcombined) {
print(plc)
}
}## Not enough samples - only 8 components will be extracted
Heatmap with hierarchical clustering
For another birds-eye view of the data, we represent it using a heatmap of the (imputed and normalized) log intensities, and cluster the samples and proteins using hierarchical clustering. In the first heatmap below, the values represent the normalized log intensities directly. In the second heatmap, the values for each protein have been centered to mean 0. The latter is also exported to a pdf file with row labels for further exploration.
ht <- makeAbundanceHeatmap(sce, assayToPlot = aNames$assayImputed, doCenter = FALSE,
settings = "report")
draw(ht, merge_legends = TRUE)
## Centered
ht <- makeAbundanceHeatmap(sce, assayToPlot = aNames$assayImputed, doCenter = TRUE,
settings = "report")
draw(ht, merge_legends = TRUE)
## Save to pdf (show row names, but no row dendrogram, order samples by group)
pdf(sub("\\.Rmd$", "_heatmap_centered.pdf", knitr::current_input(dir = TRUE)),
width = ncol(sce)/3.3 + 3.375, height = nrow(sce)/6.33 + 4)
ht <- makeAbundanceHeatmap(sce, assayToPlot = aNames$assayImputed, doCenter = TRUE,
settings = "export")
draw(ht, merge_legends = TRUE, heatmap_legend_side = "bottom",
annotation_legend_side = "bottom")
dev.off()Correlation plot
The plot below shows the pairwise Pearson correlations between all pairs of samples, based on the log2_LFQ.intensity assay.
plotassay <- assay(sce, aNames$assayImputed)
ggplot(data = as.data.frame(cor(plotassay, method = "pearson")) %>%
rownames_to_column("sample1") %>%
tidyr::pivot_longer(names_to = "sample2", values_to = "correlation",
-"sample1"),
aes(x = .data$sample1, y = .data$sample2, fill = .data$correlation)) +
geom_tile(color = "grey95", linewidth = 0.1) +
scale_fill_gradient(low = "grey95", high = "darkblue") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
labs(x = "", y = "") +
scale_x_discrete(expand = c(0, 0)) +
scale_y_discrete(expand = c(0, 0))
Save SingleCellExperiment object
The SingleCellExperiment object created above is saved
in the following location:
sceFile <- sub("\\.Rmd$", paste0("_sce.rds"), knitr::current_input(dir = TRUE))
saveRDS(sce, file = sceFile)
sceFile## [1] "/Users/charlottesoneson/Documents/Rpackages/einprot/einprot_examples/ostapcuk_2018/ostapcuk_2018_einprot_0.7.4/ostapcuk_2018_einprot_0.7.4_sce.rds"In addition, all feature information (the rowData of the
SingleCellExperiment) is written to a text file:
textFile <- sub("\\.Rmd$", paste0("_feature_info.txt"),
knitr::current_input(dir = TRUE))
write.table(as.data.frame(rowData(sce)) %>%
rownames_to_column("FeatureID"),
file = textFile, row.names = FALSE, col.names = TRUE,
quote = FALSE, sep = "\t")
textFile## [1] "/Users/charlottesoneson/Documents/Rpackages/einprot/einprot_examples/ostapcuk_2018/ostapcuk_2018_einprot_0.7.4/ostapcuk_2018_einprot_0.7.4_feature_info.txt"Explore the data interactively
For interactive exploration of your results, we generate a script to launch an adapted iSEE interface. The script can be sourced from an R console:
source('/Users/charlottesoneson/Documents/Rpackages/einprot/einprot_examples/ostapcuk_2018/ostapcuk_2018_einprot_0.7.4/ostapcuk_2018_einprot_0.7.4_iSEE.R')That will open up an iSEE session where you can interactively explore your data.
Session info
This report was compiled with the following package versions:
Click to expand
## R version 4.3.1 (2023-06-16)
## Platform: aarch64-apple-darwin20 (64-bit)
## Running under: macOS Ventura 13.4
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.11.0
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: Europe/Zurich
## tzcode source: internal
##
## attached base packages:
## [1] grid stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] einprot_0.7.4 ComplexHeatmap_2.16.0
## [3] plotly_4.10.2 ggalt_0.4.0
## [5] BiocSingular_1.16.0 scater_1.28.0
## [7] scuttle_1.10.1 SingleCellExperiment_1.22.0
## [9] tibble_3.2.1 ggplot2_3.4.2
## [11] ComplexUpset_1.3.3 dplyr_1.1.2
## [13] htmltools_0.5.5 cowplot_1.1.1
## [15] ExploreModelMatrix_1.12.0 limma_3.56.2
## [17] DT_0.28 SummarizedExperiment_1.30.2
## [19] Biobase_2.60.0 GenomicRanges_1.52.0
## [21] GenomeInfoDb_1.36.1 IRanges_2.34.1
## [23] S4Vectors_0.38.1 BiocGenerics_0.46.0
## [25] MatrixGenerics_1.12.2 matrixStats_1.0.0
## [27] STRINGdb_2.12.0
##
## loaded via a namespace (and not attached):
## [1] ProtGenerics_1.32.0 bitops_1.0-7
## [3] webshot_0.5.5 httr_1.4.6
## [5] ash_1.0-15 RColorBrewer_1.1-3
## [7] doParallel_1.0.17 tools_4.3.1
## [9] utf8_1.2.3 R6_2.5.1
## [11] lazyeval_0.2.2 mgcv_1.8-42
## [13] GetoptLong_1.0.5 withr_2.5.0
## [15] iSEEhex_1.2.0 GGally_2.1.2
## [17] gridExtra_2.3 cli_3.6.1
## [19] Cairo_1.6-0 shinyjs_2.1.0
## [21] sandwich_3.0-2 labeling_0.4.2
## [23] sass_0.4.6 robustbase_0.99-0
## [25] mvtnorm_1.2-2 readr_2.1.4
## [27] genefilter_1.82.1 systemfonts_1.0.4
## [29] svglite_2.1.1 stringdist_0.9.10
## [31] rrcov_1.7-4 plotrix_3.8-2
## [33] maps_3.4.1 rstudioapi_0.15.0
## [35] impute_1.74.1 RSQLite_2.3.1
## [37] generics_0.1.3 shape_1.4.6
## [39] crosstalk_1.2.0 gtools_3.9.4
## [41] Matrix_1.5-4.1 ggbeeswarm_0.7.2
## [43] fansi_1.0.4 imputeLCMD_2.1
## [45] lifecycle_1.0.3 yaml_2.3.7
## [47] iSEEu_1.12.0 gplots_3.1.3
## [49] blob_1.2.4 promises_1.2.0.1
## [51] crayon_1.5.2 shinydashboard_0.7.2
## [53] miniUI_0.1.1.1 lattice_0.21-8
## [55] msigdbr_7.5.1 beachmat_2.16.0
## [57] annotate_1.78.0 KEGGREST_1.40.0
## [59] magick_2.7.4 pillar_1.9.0
## [61] knitr_1.43.1 rjson_0.2.21
## [63] codetools_0.2-19 glue_1.6.2
## [65] ggiraph_0.8.7 pcaMethods_1.92.0
## [67] data.table_1.14.8 MultiAssayExperiment_1.26.0
## [69] vctrs_0.6.3 png_0.1-8
## [71] gtable_0.3.3 gsubfn_0.7
## [73] cachem_1.0.8 xfun_0.39
## [75] S4Arrays_1.0.4 mime_0.12
## [77] pcaPP_2.0-3 survival_3.5-5
## [79] rrcovNA_0.5-0 iterators_1.0.14
## [81] gmm_1.8 iSEE_2.12.0
## [83] ellipsis_0.3.2 nlme_3.1-162
## [85] bit64_4.0.5 bslib_0.5.0
## [87] tmvtnorm_1.5 irlba_2.3.5.1
## [89] vipor_0.4.5 KernSmooth_2.23-21
## [91] colorspace_2.1-0 DBI_1.1.3
## [93] proDA_1.14.0 tidyselect_1.2.0
## [95] bit_4.0.5 compiler_4.3.1
## [97] extrafontdb_1.0 chron_2.3-61
## [99] rvest_1.0.3 BiocNeighbors_1.18.0
## [101] xml2_1.3.5 DelayedArray_0.26.6
## [103] colourpicker_1.2.0 scales_1.2.1
## [105] caTools_1.18.2 DEoptimR_1.0-14
## [107] proj4_1.0-12 hexbin_1.28.3
## [109] stringr_1.5.0 digest_0.6.33
## [111] rmarkdown_2.23 XVector_0.40.0
## [113] pkgconfig_2.0.3 extrafont_0.19
## [115] sparseMatrixStats_1.12.0 highr_0.10
## [117] fastmap_1.1.1 rlang_1.1.1
## [119] GlobalOptions_0.1.2 htmlwidgets_1.6.2
## [121] shiny_1.7.4 DelayedMatrixStats_1.22.0
## [123] farver_2.1.1 jquerylib_0.1.4
## [125] zoo_1.8-12 jsonlite_1.8.7
## [127] mclust_6.0.0 BiocParallel_1.34.2
## [129] RCurl_1.98-1.12 magrittr_2.0.3
## [131] kableExtra_1.3.4 GenomeInfoDbData_1.2.10
## [133] patchwork_1.1.2 munsell_0.5.0
## [135] Rcpp_1.0.11 babelgene_22.9
## [137] viridis_0.6.3 proto_1.0.0
## [139] MsCoreUtils_1.13.1 sqldf_0.4-11
## [141] stringi_1.7.12 rintrojs_0.3.2
## [143] zlibbioc_1.46.0 MASS_7.3-60
## [145] plyr_1.8.8 ggseqlogo_0.1
## [147] parallel_4.3.1 ggrepel_0.9.3
## [149] forcats_1.0.0 Biostrings_2.68.1
## [151] splines_4.3.1 hash_2.2.6.2
## [153] hms_1.1.3 circlize_0.4.15
## [155] igraph_1.5.0 uuid_1.1-0
## [157] QFeatures_1.10.0 ScaledMatrix_1.8.1
## [159] XML_3.99-0.14 evaluate_0.21
## [161] tzdb_0.4.0 foreach_1.5.2
## [163] httpuv_1.6.11 Rttf2pt1_1.3.12
## [165] tidyr_1.3.0 purrr_1.0.1
## [167] reshape_0.8.9 clue_0.3-64
## [169] norm_1.0-11.1 rsvd_1.0.5
## [171] xtable_1.8-4 AnnotationFilter_1.24.0
## [173] later_1.3.1 viridisLite_0.4.2
## [175] memoise_2.0.1 beeswarm_0.4.0
## [177] AnnotationDbi_1.62.2 writexl_1.4.2
## [179] cluster_2.1.4 shinyWidgets_0.7.6
## [181] shinyAce_0.4.2