Re-analysis of PXD014414: Metaplastic breast cancer¶
Sabra I. Djomehri, et al.¶
Nature Communications (2020) 11:1723¶
Data re-analysis performed by Phil Wilmarth, PSR Core OHSU¶
This notebook will be using a two-sample t-test instead of edgeR¶
July 26, 2020¶
- Overview
- Library and data loading
- Clustering and normalization
- 1. Normal versus Triple Negative
- 2. Normal versus Metaplastic
- 3. Triple Negative versus Metaplastic
- Summary
Overview¶
The data were first presented in this publication:
Djomehri, S.I., Gonzalez, M.E., da Veiga Leprevost, F., Tekula, S.R., Chang, H.Y., White, M.J., Cimino-Mathews, A., Burman, B., Basrur, V., Argani, P. and Nesvizhskii, A.I., 2020. Quantitative proteomic landscape of metaplastic breast carcinoma pathological subtypes and their relationship to triple-negative tumors. Nature communications, 11(1), pp.1-15.
The IRS adjusted reporter ion intensities from the PAW pipeline were analyzed with Jupyter notebooks, an R kernel, the Bioconductor package edgeR, and its TMM normalization method. This notebook format allows rich data visualization and quality control checks.
Load the necessary R libraries¶
In [1]:
# library imports
library(tidyverse)
library(scales)
library(limma)
library(edgeR)
library(psych)
Define common functions for notebook use¶
In [2]:
# ============== CV function ===================================================
CV <- function(df) {
# Computes CVs of data frame rows
# df - data frame,
# returns vector of CVs (%)
ave <- rowMeans(df) # compute averages
sd <- apply(df, 1, sd) # compute standard deviations
cv <- 100 * sd / ave # compute CVs in percent (last thing gets returned)
}
# =========== Boxplot with median label ========================================
labeled_boxplot <- function(df, ylim, title) {
# Makes a box plot with the median value labeled
# df - data frame with data to compute CVs of
# ylim - upper limit for y-axis
# title - plot title
cv = CV(df)
boxplot(cv, ylim = c(0, ylim), notch = TRUE, main = title)
text(x = 0.65, y = boxplot.stats(cv)$stats[3],
labels = round(boxplot.stats(cv)$stats[3], 1))
}
# ================== TMM normalization from DGEList object =====================
apply_tmm_factors <- function(y, color = NULL, plot = TRUE) {
# computes the tmm normalized data from the DGEList object
# y - DGEList object
# returns a dataframe with normalized intensities
# compute and print "Sample loading" normalization factors
lib_facs <- mean(y$samples$lib.size) / y$samples$lib.size
cat("\nLibrary size factors:\n",
sprintf("%-5s -> %f\n", colnames(y$counts), lib_facs))
# compute and print TMM normalization factors
tmm_facs <- 1/y$samples$norm.factors
cat("\nTrimmed mean of M-values (TMM) factors:\n",
sprintf("%-5s -> %f\n", colnames(y$counts), tmm_facs))
# compute and print the final correction factors
norm_facs <- lib_facs * tmm_facs
cat("\nCombined (lib size and TMM) normalization factors:\n",
sprintf("%-5s -> %f\n", colnames(y$counts), norm_facs))
# compute the normalized data as a new data frame
tmt_tmm <- as.data.frame(sweep(y$counts, 2, norm_facs, FUN = "*"))
colnames(tmt_tmm) <- str_c(colnames(y$counts), "_tmm")
# visualize results and return data frame
if(plot == TRUE) {
boxplot(log10(tmt_tmm), col = color, notch = TRUE, main = "TMM Normalized data")
}
tmt_tmm
}
# ================= reformat edgeR test results ================================
collect_results <- function(df, tt, x, xlab, y, ylab) {
# Computes new columns and extracts some columns to make results frame
# df - data in data.frame
# tt - top tags table from edgeR test
# x - columns for first condition
# xlab - label for x
# y - columns for second condition
# ylab - label for y
# returns a new dataframe
# condition average vectors
ave_x <- rowMeans(df[x])
ave_y <- rowMeans(df[y])
# FC, direction, candidates
fc <- ifelse(ave_y > ave_x, (ave_y / ave_x), (-1 * ave_x / ave_y))
direction <- ifelse(ave_y > ave_x, "up", "down")
candidate <- cut(tt$FDR, breaks = c(-Inf, 0.01, 0.05, 0.10, 1.0),
labels = c("high", "med", "low", "no"))
# make data frame
temp <- cbind(df[c(x, y)], data.frame(logFC = tt$logFC, FC = fc,
PValue = tt$PValue, FDR = tt$FDR,
ave_x = ave_x, ave_y = ave_y,
direction = direction, candidate = candidate,
Acc = tt$genes))
# fix column headers for averages
names(temp)[names(temp) %in% c("ave_x", "ave_y")] <- str_c("ave_", c(xlab, ylab))
temp # return the data frame
}
# =============== p-value plots ================================================
pvalue_plots <- function(results, ylim, title) {
# Makes p-value distribution plots
# results - results data frame
# ylim - ymax for expanded view
# title - plot title
p_plot <- ggplot(results, aes(PValue)) +
geom_histogram(bins = 100, fill = "white", color = "black") +
geom_hline(yintercept = mean(hist(results$PValue, breaks = 100,
plot = FALSE)$counts[26:100]))
# we will need an expanded plot
p1 <- p_plot + ggtitle(str_c(title, " p-value distribution"))
p2 <- p_plot + coord_cartesian(xlim = c(0, 1.0), ylim = c(0, ylim)) + ggtitle("p-values expanded")
grid.arrange(p1, p2, nrow = 2) # from gridExtra package
}
# ============= log2 fold-change distributions =================================
log2FC_plots <- function(results, range, title) {
# Makes faceted log2FC plots by candidate
# results - results data frame
# range - plus/minus log2 x-axis limits
# title - plot title
ggplot(results, aes(x = logFC, fill = candidate)) +
geom_histogram(binwidth=0.1, color = "black") +
facet_wrap(~candidate) +
ggtitle(title) +
coord_cartesian(xlim = c(-range, range))
}
# ========== Setup for MA and volcano plots ====================================
transform <- function(results, x, y) {
# Make data frame with some transformed columns
# results - results data frame
# x - columns for x condition
# y - columns for y condition
# return new data frame
df <- data.frame(log10((results[x] + results[y])/2),
log2(results[y] / results[x]),
results$candidate,
-log10(results$FDR))
colnames(df) <- c("A", "M", "candidate", "P")
df # return the data frame
}
# ========== MA plots using ggplot =============================================
MA_plots <- function(results, x, y, title) {
# makes MA-plot DE candidate ggplots
# results - data frame with edgeR results and some condition average columns
# x - string for x-axis column
# y - string for y-axis column
# title - title string to use in plots
# returns a list of plots
# uses transformed data
temp <- transform(results, x, y)
# 2-fold change lines
ma_lines <- list(geom_hline(yintercept = 0.0, color = "black"),
geom_hline(yintercept = 1.0, color = "black", linetype = "dotted"),
geom_hline(yintercept = -1.0, color = "black", linetype = "dotted"))
# make main MA plot
ma <- ggplot(temp, aes(x = A, y = M)) +
geom_point(aes(color = candidate, shape = candidate)) +
scale_y_continuous(paste0("logFC (", y, "/", x, ")")) +
scale_x_continuous("Ave_intensity") +
ggtitle(title) +
ma_lines
# make separate MA plots
ma_facet <- ggplot(temp, aes(x = A, y = M)) +
geom_point(aes(color = candidate, shape = candidate)) +
scale_y_continuous(paste0("log2 FC (", y, "/", x, ")")) +
scale_x_continuous("log10 Ave_intensity") +
ma_lines +
facet_wrap(~ candidate) +
ggtitle(str_c(title, " (separated)"))
# make the plots visible
print(ma)
print(ma_facet)
}
# ========== Scatter plots using ggplot ========================================
scatter_plots <- function(results, x, y, title) {
# makes scatter-plot DE candidate ggplots
# results - data frame with edgeR results and some condition average columns
# x - string for x-axis column
# y - string for y-axis column
# title - title string to use in plots
# returns a list of plots
# 2-fold change lines
scatter_lines <- list(geom_abline(intercept = 0.0, slope = 1.0, color = "black"),
geom_abline(intercept = 0.301, slope = 1.0, color = "black", linetype = "dotted"),
geom_abline(intercept = -0.301, slope = 1.0, color = "black", linetype = "dotted"),
scale_y_log10(),
scale_x_log10())
# make main scatter plot
scatter <- ggplot(results, aes_string(x, y)) +
geom_point(aes(color = candidate, shape = candidate)) +
ggtitle(title) +
scatter_lines
# make separate scatter plots
scatter_facet <- ggplot(results, aes_string(x, y)) +
geom_point(aes(color = candidate, shape = candidate)) +
scatter_lines +
facet_wrap(~ candidate) +
ggtitle(str_c(title, " (separated)"))
# make the plots visible
print(scatter)
print(scatter_facet)
}
# ========== Volcano plots using ggplot ========================================
volcano_plot <- function(results, x, y, title) {
# makes a volcano plot
# results - a data frame with edgeR results
# x - string for the x-axis column
# y - string for y-axis column
# title - plot title string
# uses transformed data
temp <- transform(results, x, y)
# build the plot
ggplot(temp, aes(x = M, y = P)) +
geom_point(aes(color = candidate, shape = candidate)) +
xlab("log2 FC") +
ylab("-log10 FDR") +
ggtitle(str_c(title, " Volcano Plot"))
}
# ============== individual protein expression plots ===========================
# function to extract the identifier part of the accesssion
get_identifier <- function(accession) {
identifier <- str_split(accession, "\\|", simplify = TRUE)
identifier[,3]
}
set_plot_dimensions <- function(width_choice, height_choice) {
options(repr.plot.width=width_choice, repr.plot.height=height_choice)
}
plot_top_tags <- function(results, nleft, nright, top_tags) {
# results should have data first, then test results (two condition summary table)
# nleft, nright are number of data points in each condition
# top_tags is number of up and number of down top DE candidates to plot
# get top ipregulated
up <- results %>%
filter(logFC >= 0) %>%
arrange(FDR)
up <- up[1:top_tags, ]
# get top down regulated
down <- results %>%
filter(logFC < 0) %>%
arrange(FDR)
down <- down[1:top_tags, ]
# pack them
proteins <- rbind(up, down)
color = c(rep("red", nleft), rep("blue", nright))
for (row_num in 1:nrow(proteins)) {
row <- proteins[row_num, ]
vec <- as.vector(unlist(row[1:(nleft + nright)]))
names(vec) <- colnames(row[1:(nleft + nright)])
title <- str_c(get_identifier(row$Acc), ", int: ", scientific(mean(vec), 2),
", FDR: ", scientific(row$FDR, digits = 3),
", FC: ", round(row$FC, digits = 1),
", ", row$candidate)
barplot(vec, col = color, main = title,
cex.main = 1.0, cex.names = 0.7, cex.lab = 0.7)
}
}
Functions for the t-test¶
In [3]:
ttest_results <- function(df, x, y, accessions) {
# Runs t-test and makes something like a "tt" dataframe
# df - data in data.frame
# x - columns for first condition
# y - columns for second condition
# accesssions - protein accessions
# returns a new dataframe
# add average ratio columns (non-logged ratios), fold-change column, and row names
logFC <- log2(rowMeans(df[, x]) / rowMeans(df[, y]))
# apply the basic two-sample t-test (we will pool variance)
t.result <- apply(log2(df), 1, function(r) t.test(r[x], r[y], var.equal = TRUE))
# extract the p-value column from the t-test thingy
PValue <- unlist(lapply(t.result, function(x) x$p.value))
# do a Benjamini-Hochberg multiple testing correction
FDR <- p.adjust(PValue, method = "BH")
temp <- data.frame(logFC = logFC, PValue = PValue, FDR = FDR, genes = accessions)
}
count_candidates <- function(results, cutoff) {
# Counts DE candidates give a desired FDR cutoff
all <- dim(results)[1]
up <- dim(results[(results$FDR <= cutoff) & (results$logFC > 0.0), ])[1]
down <- dim(results[(results$FDR <= 0.10) & (results$logFC <= 0.0), ])[1]
cat("This is like decideTest in edgeR:\n")
cat("Up", up, "\n")
cat("not DE:", (all - up - down), "\n")
cat("Down:", down, "\n")
}
Load the IRS-normalized TMT intensity data¶
The pandas Python script that does the IRS normalization arranges the tab-delimited table so that importing into R is relatively straightforward. Note than this experimental design has only a single reference channel per plex instead of two.
We need to drop contaminant and decoy proteins, and proteins with missing sets of reporter ions. We extract the accessions column, the SL and IRS normed data columns, and collect channels by biological condition. We also grab the reference channels to check the IRS procedure.
In [4]:
# load the IRS-normalized data and check the table
data_import <- read_tsv("labeled_grouped_protein_summary_TMT_9_IRS_normalized.txt", guess_max = 5250)
# "Filter" column flags contams and decoys
# "Missing" column flags proteins without reporter ion intensities (full sets missing)
# the table from pandas is sorted so the rows we want come first
data_all <- filter(data_import, is.na(Filter), is.na(Missing))
data_sl <- data_all %>% select(., contains("SLNorm_")) %>%
select(., -contains("_Unused")) %>%
select(., -contains("_Pool"))
data_irs <- data_all %>% select(., contains("IRSNorm_")) %>%
select(., -contains("_Unused")) %>%
select(., -contains("_Pool"))
data_pool <- data_all %>% select(., contains("_Pool_"))
In [5]:
# save a few columns for the results table
all_results <- data_all %>% select(., ProtGroup, Counter, Accession, Description, starts_with("PSMs_Used"))
# save gene names for edgeR so we can double check that results line up
accessions <- data_all$Accession
# see how many rows of data we have
length(accessions)
Group samples into the three types of tissue¶
We will lump all the metaplastic breast cancer (MBC) samples together to contrast to triple negative breast cancer and normal tissue samples. There will be 14 mbc, 6 triple negative, and 7 normal samples.
In [6]:
# all categories of metaplastic breast cancer tissue
mbc_sl <- select(data_sl, contains("_C"), contains("_SP"), contains("_SQ"))
mbc_irs <- select(data_irs, contains("_C"), contains("_SP"), contains("_SQ"))
# triple negative breast cancer tissue
tn_sl <- select(data_sl, contains("_TN"))
tn_irs <- select(data_irs, contains("_TN"))
# Normal tissue
n_sl <- select(data_sl, contains("_N"))
n_irs <- select(data_irs, contains("_N"))
In [7]:
# collect the biological replicate channels
bio_sl <- cbind(n_sl, tn_sl, mbc_sl)
bio_irs <- cbind(n_irs, tn_irs, mbc_irs)
# set a color vector for plots
color <- c(rep("red", 7), rep("blue", 6), rep("green", 14))
set_plot_dimensions(9, 9)
Normalized the data¶
Run TMM normalization on the IRS data so we have a starting point for the t-test
In [8]:
# put groups together into a single data frame
tmt_sl <- bio_sl
tmt_irs <- bio_irs
# define the positions of the groups
N <- 1:7
TN <- 8:13
MP <- 14:27
# set some colors by condition
group <- c(rep("N", 7), rep("TN", 6), rep("MBC", 14))
# get the biological sample data into a DGEList object
y <- DGEList(counts = tmt_irs, group = group, genes = accessions)
# run TMM normalization (also includes a library size factor)
y <- calcNormFactors(y)
tmt_tmm <- apply_tmm_factors(y, color = color)
# check the clustering
plotMDS(y, col = color, main = "all samples after TMM")
Statistical testing starts here¶
We will run pairwise 2-sample 2-tailed t-test with equal variance. The data will be the log2 normalized intensities.
(1) Normal versus Triple Negative¶
Compute averages, flag candidates, and save the pair-wise test results¶
We will add the statistical testing results (logFC, p-values, and FDR), condition intensity averages, and candidate status to the results data frame (which has the TMM-normalized data) and also accumulate all of the three comparisons into all_results.
In [9]:
# the exact test object has columns like fold-change, CPM, and p-values
tt <- ttest_results(tmt_tmm, N, TN, accessions)
# get the results summary
results <- collect_results(tmt_tmm, tt, N, "N", TN, "TN")
# this counts up, down, and unchanged genes (proteins) at 10% FDR
count_candidates(results, cutoff = 0.05)
# check the p-value distribution
ggplot(tt, aes(PValue)) +
geom_histogram(bins = 100, fill = "white", color = "black") +
geom_hline(yintercept = mean(hist(tt$PValue, breaks = 100,
plot = FALSE)$counts[26:100])) +
ggtitle("N versus TN p-value distribution")
We have 580 DE candidates¶
We have a modest number of DE candidates. The p-value distribution suggests that the low p-value candidates are "soft" (not really small p-values) because the "spike" is more of a narrow distribution.
In [10]:
# make column names unique by adding comparison (for the accumulated frame)
results_temp <- results
colnames(results_temp) <- str_c(colnames(results), "_N_TN")
# accumulate the testing results
all_results <- cbind(all_results, results_temp)
Count candidates and look at fold-change distributions¶
We will still make the MA plots, scatter plots, and volcano plot using ggplot2. We will also look at the distribution of log2 expression ratios separated by differential expression (DE) category. The Benjamini-Hochberg corrected edgeR p-values are used to categorize the DE candidates: no > 0.10 > low > 0.05 > med > 0.01 > high.
In [11]:
# see how many candidates by category
results %>% count(candidate)
# plot log2 fold-changes by category
ggplot(results, aes(x = logFC, fill = candidate)) +
geom_histogram(binwidth=0.1, color = "black") +
facet_wrap(~candidate) +
coord_cartesian(xlim = c(-4, 4)) +
ggtitle("N vs TN logFC distributions by candidate")
Main summary plots¶
We have many comparisons to visualize, so we will use functions to generate a series of plots. We will make: an MA plot with candidates highlighted by color, faceted MA plots separated by candidate status, a scatter plot with candidates highlighted by color, faceted scatter plots separated by candidate status, and a volcano plot with candidates highlighted by color. The solid black lines in the MA and scatter plots are the 1-to-1 lines; the dotted lines are 2-fold change lines.
MA plots¶
In [12]:
# make MA plots
MA_plots(results, "ave_N", "ave_TN", "N vs TN")
Pattern of non-DE candidates is qualitatively similar to DE candidates¶
We do not see a distribution of non DE proteins that is more distinct compared to the low and medium DE candidates. Also of note, we do not have any highly significant candidates (the high category is only less than 0.01 FDR).
Scatter plots¶
In [13]:
# make scatter plots
scatter_plots(results, "ave_N", "ave_TN", "N vs TN")
Volcano plot¶
In [14]:
# make a volcano plot
volcano_plot(results, "ave_N", "ave_TN", "N vs TN")
Candidates are not super convincing looking¶
Although we have quite a few DE candidates, they are not very distinct in their expression pattern compared to the non-DE candidates.
Check some individual protein expression¶
We can look at the data for the top 20 (by BH-corrected p-values) up and down regulated proteins.
In [15]:
# look at the top 20 candidates in each direction (up in TN, then down in TN)
set_plot_dimensions(9, 3.5)
plot_top_tags(results, 7, 6, 20)
set_plot_dimensions(9, 9)
(2) Normal versus Metaplastic¶
We will do the same testing to compare normal tissue to metaplastic breast cancer samples.
In [16]:
# the test object has columns like fold-change, CPM, and p-values
tt <- ttest_results(tmt_tmm, N, MP, accessions)
# get the results summary
results <- collect_results(tmt_tmm, tt, N, "N", MP, "MP")
# this counts up, down, and unchanged genes (proteins) at 10% FDR
count_candidates(results, cutoff = 0.05)
# check the p-value distribution
ggplot(tt, aes(PValue)) +
geom_histogram(bins = 100, fill = "white", color = "black") +
geom_hline(yintercept = mean(hist(tt$PValue, breaks = 100,
plot = FALSE)$counts[26:100])) +
ggtitle("N versus MP p-value distribution")
We have more candidates¶
We have around 1,000 candidates in this comparison. The p-value distribution looks better.
Compute averages, flag candidates, and save the pair-wise test results¶
We will add the statistical testing results (logFC, p-values, and FDR), condition intensity averages, and candidate status to the results data frame (which has the TMM-normalized data) and also accumulate all of the three comparisons into all_results.
In [17]:
# make column names unique by adding comparison
results_temp <- results
colnames(results_temp) <- str_c(colnames(results), "_N_MP")
# accumulate the testing results
all_results <- cbind(all_results, results_temp)
Count candidates and look at fold-change distributions¶
In [18]:
# see how many candidates by category
results %>% count(candidate)
# plot log2 fold-changes by category
ggplot(results, aes(x = logFC, fill = candidate)) +
geom_histogram(binwidth=0.1, color = "black") +
facet_wrap(~candidate) +
coord_cartesian(xlim = c(-4, 4)) +
ggtitle("N vs MP logFC distributions by candidate")
MA plots¶
In [19]:
# make MA plots
MA_plots(results, "ave_N", "ave_MP", "N vs MP")
Scatter plots¶
In [20]:
# make scatter plots
scatter_plots(results, "ave_N", "ave_MP", "N vs MP")
Volcano plot¶
In [21]:
# make a volcano plot
volcano_plot(results, "ave_N", "ave_MP", "N vs MP")
In [22]:
# look at the top 20 candidates (up in MP, then down in MP)
set_plot_dimensions(9, 3.5)
plot_top_tags(results, 7, 14, 20)
set_plot_dimensions(9, 9)
The top DE candidates look more convincing¶
There are some sample to sample variations in line with what one expects for human subjects. We do seem to have clear differences in expression levels.
(3) Triple Negative versus Metaplastic¶
We can compare the non-metaplastic breast cancer samples to the metaplastic samples. These samples were not very distinct from each other in the cluster plots.
In [23]:
# the test object has columns like fold-change, CPM, and p-values
tt <- ttest_results(tmt_tmm, TN, MP, accessions)
# get the results summary
results <- collect_results(tmt_tmm, tt, TN, "TN", MP, "MP")
# this counts up, down, and unchanged genes (proteins) at 10% FDR
count_candidates(results, cutoff = 0.05)
# check the p-value distribution
ggplot(tt, aes(PValue)) +
geom_histogram(bins = 100, fill = "white", color = "black") +
geom_hline(yintercept = mean(hist(tt$PValue, breaks = 100,
plot = FALSE)$counts[26:100])) +
ggtitle("TN versus MP p-value distribution")
We have no candidates¶
There is no point in plotting things.
In [24]:
# make column names unique by adding comparison
results_temp <- results
colnames(results_temp) <- str_c(colnames(results), "_TN_MP")
# accumulate the testing results
all_results <- cbind(all_results, results_temp)
MA plots¶
In [25]:
# make MA plots
MA_plots(results, "ave_TN", "ave_MP", "TN vs MP")
Scatter plots¶
In [26]:
# make scatter plots
scatter_plots(results, "ave_TN", "ave_MP", "TN vs MP")
Volcano plot¶
In [27]:
# make a volcano plot
volcano_plot(results, "ave_TN", "ave_MP", "TN vs MP")
Check individual protein expression¶
These will all be non-candidates.
In [28]:
# look at the top 10 candidates (up in MP, then down in MP)
set_plot_dimensions(9, 3.5)
plot_top_tags(results, 6, 14, 10)
set_plot_dimensions(9, 9)
Summary¶
The t-test seems to be a bit messier than the edgeR or limma testing. The moderated testing in those other packages is probably helping.
Save the all_results frame to TSV file¶
In [29]:
# write the results to disk
write.table(all_results, "results_t-test.txt", sep = "\t",
row.names = FALSE, na = " ")
Log the session information¶
In [30]:
# log the session details
sessionInfo()
In [ ]: