Ostapcuk et al 2018
FragPipe analysis
2023-07-19 20:46:22 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"
fragpipeDir <-
"/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/einprot/extdata/fp_example"
idCol <-
function (df)
combineIds(df, combineCols = c("Gene", "Protein.ID"))
labelCol <-
function (df)
combineIds(df, combineCols = c("Gene", "Protein.ID"))
geneIdCol <-
function (df)
getFirstId(df, colName = "Gene")
proteinIdCol <-
"Protein.ID"
stringIdCol <-
function (df)
combineIds(df, combineCols = c("Gene", "Protein.ID"), combineWhen = "missing",
makeUnique = FALSE)
reportTitle <-
"Ostapcuk et al 2018"
reportAuthor <-
"FragPipe analysis"
iColPattern <-
".MaxLFQ.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 <-
NULL
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 FragPipe. 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 |
FragPipe analysis summary
fp <- readFragPipeInfo(fragpipeDir)
fp <- c(list("Protein file" = file.path(fragpipeDir, "combined_protein.tsv")), fp)
makeTableFromList(fp)| Protein file | /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/einprot/extdata/fp_example/combined_protein.tsv |
| FragPipe version | 19.1 |
| FragPipe parameter file | /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/einprot/extdata/fp_example/fragpipe.workflow |
| FragPipe log file | /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/einprot/extdata/fp_example/log_2023-04-12_20-12-46.txt |
| Search engine | MSFragger-3.7 |
| Raw file location | D:/Data/FUSION |
| Raw files | F_160817_AdnpFB_IP06.raw, F_160817_AdnpFB_IP05.raw, F_160817_RBC_ctrl_IP02.raw, F_160817_AdnpFB_IP04.raw, F_160817_RBC_ctrl_IP01.raw, F_160817_Chd4BF_IP09.raw, F_160817_Chd4BF_IP07.raw, F_160817_Chd4BF_IP08.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/://Data//FASTA//2023-04-12-decoys-contam_MOUSE__190410.fasta.fas |
| Peptides (ranges) | length: 7-50 AA; mass: 500-5000 Da |
| Mass error tolerances | precursor:-20-20 [ppm]; fragment:0.7 [Da] (after optimization:200 PPM) |
| Quantification settings (LFQ) | IonQuant: TRUE, Calculate MaxLFQ intensity: TRUE, Normalization: TRUE, match-between runs (MBR): FALSE, min. ions: 2 |
| Enzymes | stricttrypsin[KR, C-terminal, 2 missed cleavages] |
| Variable modifications | M(15.9949), N-term(42.0106) |
| Fixed modifications | C(57.0215) |
| Database decoy tag | rev_ |
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 | .MaxLFQ.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 | |
| 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 FragPipe output
The input to this workflow is a combined_protein file
from FragPipe (see path in the table above). We read the FragPipe
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 FragPipe file.
## Read FragPipe output
tmp <- importExperiment(inFile = file.path(fragpipeDir, "combined_protein.tsv"),
iColPattern = iColPattern,
includeOnlySamples = includeOnlySamples,
excludeSamples = excludeSamples)
sce <- tmp$sce
aName <- tmp$aName
sce## class: SingleCellExperiment
## dim: 453 9
## metadata(1): colList
## assays(5): MaxLFQ.intensity Unique.spectral.count Total.spectral.count
## Spectral.count Intensity
## rownames(453): 1 2 ... 452 453
## rowData names(15): Protein Protein.ID ... Combined.Total.Spectral.Count
## Indistinguishable.Proteins
## 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## $MaxLFQ.intensity
## [1] "Adnp_IP04.MaxLFQ.Intensity" "Adnp_IP05.MaxLFQ.Intensity"
## [3] "Adnp_IP06.MaxLFQ.Intensity" "Chd4BF_IP07.MaxLFQ.Intensity"
## [5] "Chd4BF_IP08.MaxLFQ.Intensity" "Chd4BF_IP09.MaxLFQ.Intensity"
## [7] "RBC_ctrl_IP01.MaxLFQ.Intensity" "RBC_ctrl_IP02.MaxLFQ.Intensity"
## [9] "RBC_ctrl_IP03.MaxLFQ.Intensity"
##
## $Unique.spectral.count
## [1] "Adnp_IP04.Unique.Spectral.Count" "Adnp_IP05.Unique.Spectral.Count"
## [3] "Adnp_IP06.Unique.Spectral.Count" "Chd4BF_IP07.Unique.Spectral.Count"
## [5] "Chd4BF_IP08.Unique.Spectral.Count" "Chd4BF_IP09.Unique.Spectral.Count"
## [7] "RBC_ctrl_IP01.Unique.Spectral.Count" "RBC_ctrl_IP02.Unique.Spectral.Count"
## [9] "RBC_ctrl_IP03.Unique.Spectral.Count"
##
## $Total.spectral.count
## [1] "Adnp_IP04.Total.Spectral.Count" "Adnp_IP05.Total.Spectral.Count"
## [3] "Adnp_IP06.Total.Spectral.Count" "Chd4BF_IP07.Total.Spectral.Count"
## [5] "Chd4BF_IP08.Total.Spectral.Count" "Chd4BF_IP09.Total.Spectral.Count"
## [7] "RBC_ctrl_IP01.Total.Spectral.Count" "RBC_ctrl_IP02.Total.Spectral.Count"
## [9] "RBC_ctrl_IP03.Total.Spectral.Count"
##
## $Spectral.count
## [1] "Combined.Unique.Spectral.Count" "Combined.Total.Spectral.Count"
## [3] "Adnp_IP04.Spectral.Count" "Adnp_IP05.Spectral.Count"
## [5] "Adnp_IP06.Spectral.Count" "Chd4BF_IP07.Spectral.Count"
## [7] "Chd4BF_IP08.Spectral.Count" "Chd4BF_IP09.Spectral.Count"
## [9] "RBC_ctrl_IP01.Spectral.Count" "RBC_ctrl_IP02.Spectral.Count"
## [11] "RBC_ctrl_IP03.Spectral.Count" "Adnp_IP04.Unique.Spectral.Count"
## [13] "Adnp_IP05.Unique.Spectral.Count" "Adnp_IP06.Unique.Spectral.Count"
## [15] "Chd4BF_IP07.Unique.Spectral.Count" "Chd4BF_IP08.Unique.Spectral.Count"
## [17] "Chd4BF_IP09.Unique.Spectral.Count" "RBC_ctrl_IP01.Unique.Spectral.Count"
## [19] "RBC_ctrl_IP02.Unique.Spectral.Count" "RBC_ctrl_IP03.Unique.Spectral.Count"
## [21] "Adnp_IP04.Total.Spectral.Count" "Adnp_IP05.Total.Spectral.Count"
## [23] "Adnp_IP06.Total.Spectral.Count" "Chd4BF_IP07.Total.Spectral.Count"
## [25] "Chd4BF_IP08.Total.Spectral.Count" "Chd4BF_IP09.Total.Spectral.Count"
## [27] "RBC_ctrl_IP01.Total.Spectral.Count" "RBC_ctrl_IP02.Total.Spectral.Count"
## [29] "RBC_ctrl_IP03.Total.Spectral.Count"
##
## $Intensity
## [1] "Adnp_IP04.Intensity" "Adnp_IP05.Intensity"
## [3] "Adnp_IP06.Intensity" "Chd4BF_IP07.Intensity"
## [5] "Chd4BF_IP08.Intensity" "Chd4BF_IP09.Intensity"
## [7] "RBC_ctrl_IP01.Intensity" "RBC_ctrl_IP02.Intensity"
## [9] "RBC_ctrl_IP03.Intensity" "Adnp_IP04.MaxLFQ.Intensity"
## [11] "Adnp_IP05.MaxLFQ.Intensity" "Adnp_IP06.MaxLFQ.Intensity"
## [13] "Chd4BF_IP07.MaxLFQ.Intensity" "Chd4BF_IP08.MaxLFQ.Intensity"
## [15] "Chd4BF_IP09.MaxLFQ.Intensity" "RBC_ctrl_IP01.MaxLFQ.Intensity"
## [17] "RBC_ctrl_IP02.MaxLFQ.Intensity" "RBC_ctrl_IP03.MaxLFQ.Intensity"Add sample annotations
Next, we compile the sample annotations. The sample IDs have been
extracted from the column names in the FragPipe 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 | MaxLFQ.intensity |
| assayLog2WithNA | log2_MaxLFQ.intensity_withNA |
| assayImputIndic | imputed_MaxLFQ.intensity |
| assayLog2NormWithNA | log2_MaxLFQ.intensity_withNA |
| assayImputed | log2_MaxLFQ.intensity |
| assayBatchCorr | log2_MaxLFQ.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 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, reverse hits, and proteins with low confidence
Next, we filter out any features classified by FragPipe 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 <- filterFragPipe(sce = sce, minPeptides = minPeptides,
plotUpset = TRUE,
revPattern = paste0("^", fp[["Database decoy tag"]]),
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 98 features. The new sce object
has 355 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 einprotLabel IDsForSTRING
## 1 Cul4b Cul4b A2A432-2 Cul4b Cul4b
## 2 Zmym4 Zmym4 A2A791 Zmym4 Zmym4
## 3 Chd5 Chd5 A2A8L1 Chd5 Chd5
## 4 Cep152 Cep152 A2AUM9 Cep152 Cep152
## 5 Odf2 Odf2 A3KGV1-2 Odf2 Odf2
## 6 ... ... ... ... ...
## 7 Sox2 Sox2 Q60I23 Sox2 Sox2
## 8 Bms1 Bms1 Q6PGF5 Bms1 Bms1
## 9 Parp1 Parp1 Q921K2 Parp1 Parp1
## 10 Rbmxl2 Rbmxl2 Q9DAE2 Rbmxl2 Rbmxl2
## 11 Ehmt1.Z4YJZ7 Ehmt1 Z4YJZ7 Ehmt1.Z4YJZ7 Ehmt1Prepare 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 217
## [["mouse: DCS complex (Ptbp1, Ptbp2, Hnrph1, Hnrpf) (+1 alt. ID)"]] Hnrnph1 ...
## [["mouse: Drosha complex"]] Dhx15 Hnrnph1 Ddx5 Hnrnpm
## [["mouse: Gata1-Fog1-MeCP1 complex"]] Hdac1 Hdac2 ... Mta1.F8WHY8 Mta3.Q3UII8
## [["mouse: Gon4l-Hdac1-Sin3a-Yy1 complex"]] Hdac1 Yy1
## [["mouse: Hdac1-Ino80-Kdm1a-Phb2-Rbp1-Taf5 complex"]] Hdac1 Kdm1a
## [["mouse: Hsp110-Hsc70-Hsp25 complex"]] Hspb1 Hspa8
## [["mouse: Ikaros complex"]] Hdac1 Hdac2 Rbbp4 Chd4
## [["mouse: Ikaros-NuRD complex"]] Hdac1 Hdac2 Chd4
## [["mouse: MCM complex (+3 alt. IDs)"]] Mcm5 Mcm6 Mcm7
## [["mouse: Metallothionein-3 complex"]] Hsp90ab1 Actb Ckb
## ...
## <207 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: 355 9
## metadata(1): colList
## assays(7): MaxLFQ.intensity Unique.spectral.count ...
## log2_MaxLFQ.intensity_withNA imputed_MaxLFQ.intensity
## rownames(355): Cul4b Zmym4 ... Rbmxl2 Ehmt1.Z4YJZ7
## rowData names(22): Protein Protein.ID ... 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.7141355plotImputationDistribution(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_MaxLFQ.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 = "FragPipe")
## Add experiment metadata
S4Vectors::metadata(sce) <- c(
S4Vectors::metadata(sce),
list(
fragpipeDir = fragpipeDir,
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: 355 9
## metadata(16): colList iSEE ... testres comparisonList
## assays(8): MaxLFQ.intensity Unique.spectral.count ...
## imputed_MaxLFQ.intensity log2_MaxLFQ.intensity
## rownames(355): Cul4b Zmym4 ... Rbmxl2 Ehmt1.Z4YJZ7
## rowData names(88): Protein Protein.ID ...
## Chd4BF_vs_RBC_ctrl.log2_MaxLFQ.intensity.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 (203 of the 355 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_MaxLFQ.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_fp/ostapcuk_2018_fp_einprot_0.7.4/ostapcuk_2018_fp_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_fp/ostapcuk_2018_fp_einprot_0.7.4/ostapcuk_2018_fp_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_fp/ostapcuk_2018_fp_einprot_0.7.4/ostapcuk_2018_fp_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