Example dataset: Ubiquitin interactors

We analyze a proteomics dataset in which Ubiquitin-protein interactors were characterized (Zhang et al. Mol Cell 2017). The raw mass spectrometry data were first analyzed using MaxQuant (Cox and Mann, Nat Biotech 2007) and the resulting “proteinGroups.txt” file is used as input for the downstream analysis.

Loading of the data

# Loading a package required for data handling
library("dplyr")

# The data is provided with the package
data <- UbiLength

# We filter for contaminant proteins and decoy database hits, which are indicated by "+" in the columns "Potential.contaminants" and "Reverse", respectively. 
data <- filter(data, Reverse != "+", Potential.contaminant != "+")

The data.frame has the following dimensions:

dim(data)

## [1] 2941   23

The data.frame has the following column names:

colnames(data)

##  [1] "Protein.IDs"             "Majority.protein.IDs"   
##  [3] "Protein.names"           "Gene.names"             
##  [5] "Fasta.headers"           "Peptides"               
##  [7] "Razor...unique.peptides" "Unique.peptides"        
##  [9] "LFQ.intensity.Ubi4_1"    "LFQ.intensity.Ubi4_2"   
## [11] "LFQ.intensity.Ubi4_3"    "LFQ.intensity.Ubi6_1"   
## [13] "LFQ.intensity.Ubi6_2"    "LFQ.intensity.Ubi6_3"   
## [15] "LFQ.intensity.Ctrl_1"    "LFQ.intensity.Ctrl_2"   
## [17] "LFQ.intensity.Ctrl_3"    "LFQ.intensity.Ubi1_1"   
## [19] "LFQ.intensity.Ubi1_2"    "LFQ.intensity.Ubi1_3"   
## [21] "Only.identified.by.site" "Reverse"                
## [23] "Potential.contaminant"

The “LFQ.intensity” columns will be used for subsequent analysis.

Data preparation

The dataset has unique Uniprot identifiers, however those are not immediately informative. The associated gene names are informative, however these are not always unique.

# Are there any duplicated gene names?
data$Gene.names %>% duplicated() %>% any()

## [1] TRUE

# Make a table of duplicated gene names
data %>% group_by(Gene.names) %>% summarize(frequency = n()) %>% 
  arrange(desc(frequency)) %>% filter(frequency > 1)

## # A tibble: 51 × 2
##    Gene.names frequency
##    <chr>          <int>
##  1 ""                 7
##  2 "ATXN2"            4
##  3 "ATXN2L"           4
##  4 "SF1"              4
##  5 "HSPA8"            3
##  6 "RBM33"            3
##  7 "UGP2"             3
##  8 "ACTL6A"           2
##  9 "BCLAF1"           2
## 10 "BRAP"             2
## # ℹ 41 more rows

For further analysis these proteins must get unique names. Additionally, some proteins do not have an annotated gene name and for those we will use the Uniprot ID.

# Make unique names using the annotation in the "Gene.names" column as primary names and the annotation in "Protein.IDs" as name for those that do not have an gene name.
data_unique <- make_unique(data, "Gene.names", "Protein.IDs", delim = ";")

# Are there any duplicated names?
data$name %>% duplicated() %>% any()

## [1] FALSE

Generate a SummarizedExperiment object

Many Bioconductor packages use SummarizedExperiment objects as input and/or output. This class of objects contains and coordinates the actual (assay) data, information on the samples as well as feature annotation. We can generate the SummarizedExperiment object from our data using two different approaches. We can extract sample information directly from the column names or we add sample information using an experimental design template. An experimental design is included in the package for our example dataset.

The experimental design must contain ‘label’, ‘condition’ and ‘replicate’ columns. The ‘label’ column contains the identifiers of the different samples and they should correspond to the column names containing the assay data. The ‘condition’ and ‘replicate’ columns contain the annotation of these samples as defined by the user.

# Generate a SummarizedExperiment object using an experimental design
LFQ_columns <- grep("LFQ.", colnames(data_unique)) # get LFQ column numbers
experimental_design <- UbiLength_ExpDesign
data_se <- make_se(data_unique, LFQ_columns, experimental_design)

# Generate a SummarizedExperiment object by parsing condition information from the column names
LFQ_columns <- grep("LFQ.", colnames(data_unique)) # get LFQ column numbers
data_se_parsed <- make_se_parse(data_unique, LFQ_columns)

# Let's have a look at the SummarizedExperiment object
data_se

## class: SummarizedExperiment 
## dim: 2941 12 
## metadata(0):
## assays(1): ''
## rownames(2941): RBM47 UBA6 ... ATXN2.3 X6RHB9
## rowData names(13): Protein.IDs Majority.protein.IDs ... name ID
## colnames(12): Ubi4_1 Ubi4_2 ... Ubi1_2 Ubi1_3
## colData names(4): label ID condition replicate

Prerequisites of the SummarizedExperiment object

The make_se and make_se_parse functions generate a SummarizedExperiment object that has a couple of specifications. The assay data is log2-transformed and its rownames depict the protein names. The rowData contains, amongst others, the ‘name’ and ‘ID’ columns that were generated by make_unique. The colData contains the experimental design and thereby the sample annotation. Thereby the colData includes the ‘label’, ‘condition’ and ‘replicate’ columns as well as a newly generated ‘ID’ column. The log2-transformed assay data and the specified rowData and colData columns are prerequisites for the subsequent analysis steps.

Filter on missing values

The dataset contains proteins which are not quantified in all replicates. Some proteins are even only quantified in a single replicate.

# Plot a barplot of the protein identification overlap between samples
plot_frequency(data_se)

This leaves our dataset with missing values, which need to be imputed. However, this should not be done for proteins that contain too many missing values. Therefore, we first filter out proteins that contain too many missing values. This is done by setting the threshold for the allowed number of missing values per condition in the filter_missval function.

# Filter for proteins that are identified in all replicates of at least one condition
data_filt <- filter_missval(data_se, thr = 0)

# Less stringent filtering:
# Filter for proteins that are identified in 2 out of 3 replicates of at least one condition
data_filt2 <- filter_missval(data_se, thr = 1)

After filtering, the number of identified proteins per sample can be plotted as well as the overlap in identifications between samples.

# Plot a barplot of the number of identified proteins per samples
plot_numbers(data_filt)

# Plot a barplot of the protein identification overlap between samples
plot_coverage(data_filt)

Normalization

The data is background corrected and normalized by variance stabilizing transformation (vsn).

# Normalize the data
data_norm <- normalize_vsn(data_filt)

The normalization can be inspected by checking the distributions of the samples before and after normalization.

# Visualize normalization by boxplots for all samples before and after normalization
plot_normalization(data_filt, data_norm)

Impute data for missing values

The remaining missing values in the dataset need to be imputed. The data can be missing at random (MAR), for example if proteins are quantified in some replicates but not in others. Data can also be missing not at random (MNAR), for example if proteins are not quantified in specific conditions (e.g. in the control samples). MNAR can indicate that proteins are below the detection limit in specific samples, which could be very well the case in proteomics experiments. For these different conditions, different imputation methods have to be used, as described in the MSnbase vignette and more specifically in the impute function descriptions.

To explore the pattern of missing values in the data, a heatmap is plotted indicating whether values are missing (0) or not (1). Only proteins with at least one missing value are visualized.

# Plot a heatmap of proteins with missing values
plot_missval(data_filt)

This heatmap indicates that missing values are highly biased to specific samples. The example dataset is an affinity enrichment dataset of ubiquitin interactors, which is likely to have proteins which are below the detection limit in specific samples. These can be proteins that are specifically enriched in the ubiquitin purifications, but are not enriched in the controls samples, or vice versa. To check whether missing values are biased to lower intense proteins, the densities and cumulative fractions are plotted for proteins with and without missing values.

# Plot intensity distributions and cumulative fraction of proteins with and without missing values
plot_detect(data_filt)

Indeed the proteins with missing values have on average low intensities. This data (MNAR and close to the detection limit) should be imputed by a left-censored imputation method, such as the quantile regression-based left-censored function (“QRILC”) or random draws from a left-shifted distribution (“MinProb” and “man”). In contrast, MAR data should be imputed with methods such as k-nearest neighbor (“knn”) or maximum likelihood (“MLE”) functions. See the MSnbase vignette and more specifically the impute function description for more information.

# All possible imputation methods are printed in an error, if an invalid function name is given.
impute(data_norm, fun = "")

## Error in match.arg(fun): 'arg' should be one of "bpca", "knn", "QRILC", "MLE", "MinDet", "MinProb", "man", "min", "zero", "mixed", "nbavg"

# Impute missing data using random draws from a Gaussian distribution centered around a minimal value (for MNAR)
data_imp <- impute(data_norm, fun = "MinProb", q = 0.01)

# Impute missing data using random draws from a manually defined left-shifted Gaussian distribution (for MNAR)
data_imp_man <- impute(data_norm, fun = "man", shift = 1.8, scale = 0.3)

# Impute missing data using the k-nearest neighbour approach (for MAR)
data_imp_knn <- impute(data_norm, fun = "knn", rowmax = 0.9)

The effect of the imputation on the distributions can be visualized.

# Plot intensity distributions before and after imputation
plot_imputation(data_norm, data_imp)

Differential enrichment analysis

Protein-wise linear models combined with empirical Bayes statistics are used for the differential enrichment analysis (or differential expression analysis). The test_diff function introduced here uses limma and automatically generates the contrasts to be tested. For the contrasts generation, the control sample has to be specified. Additionally, the types of contrasts to be produced need to be indicated, allowing the generation of all possible comparisons (“all”) or the generation of contrasts of every sample versus control (“control”). Alternatively, the user can manually specify the contrasts to be tested (type = “manual”), which need to be specified in the argument test.

# Differential enrichment analysis  based on linear models and empherical Bayes statistics

# Test every sample versus control
data_diff <- test_diff(data_imp, type = "control", control = "Ctrl")

## Tested contrasts: Ubi4_vs_Ctrl, Ubi6_vs_Ctrl, Ubi1_vs_Ctrl

# Test all possible comparisons of samples
data_diff_all_contrasts <- test_diff(data_imp, type = "all")

## Tested contrasts: Ubi4_vs_Ubi6, Ubi4_vs_Ctrl, Ubi4_vs_Ubi1, Ubi6_vs_Ctrl, Ubi6_vs_Ubi1, Ctrl_vs_Ubi1

# Test manually defined comparisons
data_diff_manual <- test_diff(data_imp, type = "manual", 
                              test = c("Ubi4_vs_Ctrl", "Ubi6_vs_Ctrl"))

## Tested contrasts: Ubi4_vs_Ctrl, Ubi6_vs_Ctrl

Finally, significant proteins are defined by user-defined cutoffs using add_rejections.

# Denote significant proteins based on user defined cutoffs
dep <- add_rejections(data_diff, alpha = 0.05, lfc = log2(1.5))

Visualization of the results

The results from the previous analysis can be easily visualized by a number of functions. These visualizations assist in the determination of the optimal cutoffs to be used, highlight the most interesting samples and contrasts, and pinpoint differentially enriched/expressed proteins.

# Plot the first and second principal components
plot_pca(dep, x = 1, y = 2, n = 500, point_size = 4)

## Warning: Use of `pca_df[[indicate[1]]]` is discouraged.
## ℹ Use `.data[[indicate[1]]]` instead.

## Warning: Use of `pca_df[[indicate[2]]]` is discouraged.
## ℹ Use `.data[[indicate[2]]]` instead.

# Plot the Pearson correlation matrix
plot_cor(dep, significant = TRUE, lower = 0, upper = 1, pal = "Reds")

# Plot a heatmap of all significant proteins with the data centered per protein
plot_heatmap(dep, type = "centered", kmeans = TRUE, 
             k = 6, col_limit = 4, show_row_names = FALSE,
             indicate = c("condition", "replicate"))

# Plot a heatmap of all significant proteins (rows) and the tested contrasts (columns)
plot_heatmap(dep, type = "contrast", kmeans = TRUE, 
             k = 6, col_limit = 10, show_row_names = FALSE)

# Plot a volcano plot for the contrast "Ubi6 vs Ctrl""
plot_volcano(dep, contrast = "Ubi6_vs_Ctrl", label_size = 2, add_names = TRUE)

## Warning: ggrepel: 50 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

# Plot a barplot for USP15 and IKBKG
plot_single(dep, proteins = c("USP15", "IKBKG"))

# Plot a barplot for the protein USP15 with the data centered
plot_single(dep, proteins = "USP15", type = "centered")

# Plot a frequency plot of significant proteins for the different conditions
plot_cond(dep)

Results table

To extract a table containing the essential results, the get_results function can be used.

# Generate a results table
data_results <- get_results(dep)

# Number of significant proteins
data_results %>% filter(significant) %>% nrow()

## [1] 258

The resulting table contains the following columns:

# Column names of the results table
colnames(data_results)

##  [1] "name"                     "ID"                      
##  [3] "Ubi1_vs_Ctrl_p.val"       "Ubi4_vs_Ctrl_p.val"      
##  [5] "Ubi6_vs_Ctrl_p.val"       "Ubi1_vs_Ctrl_p.adj"      
##  [7] "Ubi4_vs_Ctrl_p.adj"       "Ubi6_vs_Ctrl_p.adj"      
##  [9] "Ubi1_vs_Ctrl_significant" "Ubi4_vs_Ctrl_significant"
## [11] "Ubi6_vs_Ctrl_significant" "significant"             
## [13] "Ubi1_vs_Ctrl_ratio"       "Ubi4_vs_Ctrl_ratio"      
## [15] "Ubi6_vs_Ctrl_ratio"       "Ctrl_centered"           
## [17] "Ubi1_centered"            "Ubi4_centered"           
## [19] "Ubi6_centered"

Of these columns, the p.val and p.adj columns contain the raw and adjusted p values, respectively, for the contrast as depicted in the column name. The ratio columns contain the average log2 fold changes. The significant columns indicate whether the protein is differentially enriched/expressed, as defined by the chosen cutoffs. The centered columns contain the average log2 fold changes scaled by protein-wise centering.

Generate a data.frame from the resulting SummarizedExperiment object

You might want to obtain an ordinary data.frame of the results. For this purpose, the package provides functions to convert SummarizedExperiment objects to data.frames. get_df_wide will generate a wide table, whereas get_df_long will generate a long table.

# Generate a wide data.frame
df_wide <- get_df_wide(dep)
# Generate a long data.frame
df_long <- get_df_long(dep)

Save your results object for reuse

To facilitate future analysis and/or visualization of your current data, saving your analyzed data is highly recommended. We save the final data object (dep) as well as intermediates of the analysis, i.e. the initial SummarizedExperiment object (data_se), normalized data (data_norm), imputed data (data_imp) and differentially expression analyzed data (data_diff). This allows us to easily change parameters in future analysis.

# Save analyzed data
save(data_se, data_norm, data_imp, data_diff, dep, file = "data.RData")
# These data can be loaded in future R sessions using this command
load("data.RData")

LFQ-based DEP analysis

Differential enrichment analysis of label-free proteomics data can be performed using the LFQ workflow function.

# The data is provided with the package 
data <- UbiLength
experimental_design <- UbiLength_ExpDesign

# The wrapper function performs the full analysis
data_results <- LFQ(data, experimental_design, fun = "MinProb", 
      type = "control", control = "Ctrl", alpha = 0.05, lfc = 1)

## [1] 0.2608395

This wrapper produces a list of objects, which can be used to create a report and/or for further analysis. The report function produces two reports (pdf and html), a results table and a RData object, which are saved in a generated “Report” folder.

# Make a markdown report and save the results
report(data_results)

The results generated by the LFQ function contain, among others, the results object (data.frame object) and the dep object (SummarizedExperiment object).

# See all objects saved within the results object
names(data_results)

## [1] "data"    "se"      "filt"    "norm"    "imputed" "diff"    "dep"    
## [8] "results" "param"

The results object contains the essential results and the dep object contains the full SummarizedExperiment object. The results table can be explored by selecting the $results object

# Extract the results table
results_table <- data_results$results

# Number of significant proteins
results_table %>% filter(significant) %>% nrow()

## [1] 233

The full data (dep object) can be used for the plotting functions as described in the chapter “Visualization of the results”, for example a heatmap.

# Extract the sign object
full_data <- data_results$dep

# Use the full data to generate a heatmap
plot_heatmap(full_data, type = "contrast", kmeans = TRUE, 
             k = 6, col_limit = 4, show_row_names = FALSE)

TMT-based DEP analysis

Differential enrichment analysis of Tandem-Mass-Tag labeled proteomics data is also supported. Protein results tables from IsobarQuant can be directly analyzed using the TMT wrapper function.

# Need example data
TMTdata <- example_data
Exp_Design <- example_Exp_Design

# The wrapper function performs the full analysis
TMTdata_results <- TMT(TMTdata, expdesign = Exp_Design, fun = "MinProb",
    type = "control", control = "Control", alpha = 0.05, lfc = 1)