Differential expression analysis of proteomics data

This document details a differential expression analysis of total proteome data acquired using a TMT-MS approach with SPS-MS3 on an Orbitrap Fusion Lumos. This document begins with the processing of peptide identification search results that were obtained using the procedure outlined in the document found here.

Setting up the environment

These are packages you will need for this notebook. For exact versions used, please refer to the session info at the bottom of this notebook.

library('tidyverse')
library('ggplot2')
library('RColorBrewer')
library('vroom')
library('DEqMS')

I like to set a base directory that we can use as a link to the directory where we will do most of the work. I use two directories here because the Workspace is what is pushed to my GitHub for version tracking and contains things like scripts, but the Repository is where more of the big data is stored that does not get pushed (raw search results).

baseWorkspace = 'C:/Users/chris/OneDrive/Documents/bccrc/protocolsDryLab/relatedToProteomics/'
baseRepository = 'C:/Users/chris/OneDrive/Documents/bccrc/protocolsRepository/totalProteomeDifferentialExpressionAnalysis'

In order to compile the quantification data later on, I use a function for simplicity. I have defined this function in a separate file that I call into our session using the command below. The file referenced below with the function can be found here.

source(paste(baseWorkspace, '/totalProteomeDifferentialExpressionAnalysisUserDefinedFunctions.R', sep = ''))

The last thing we want to get is our annotated fasta index. If you followed the raw data processing outlined in the file mentioned above, this will have been created in the directory with all of the other search results. If you did not follow that processing, you can use the R script provided in that file to easily generate an index file from a fasta database.

##########################################################################################
proteinAnnotation = readRDS(paste(baseRepository, '/uniprotHumanJun2020.fasta.annotated.rds', sep = ''))

Data processing

I want to process the data using DEqMS for differential expression. I am following their guide as detailed here as it is very nice and quite extensive.

The first thing we need to do is read in the peptide search results. Below there are a few commands, but we are essentially reading our PSM result file (as generated in the raw data processing script above) and doing some basic parsing to make it more manageable. It is all based on the tidyverse package, so if you are unfamiliar with text processing using this package, please look around online as they have extensive documentation. If you have not followed the same raw data processing pipeline as I have above, your search results may look different.

##########################################################################################
psm = vroom(paste(baseRepository, '/n_10Oct2016_ADROIT_PrepCompare_TMT10_hph_1_Default_PSM_Report.txt', sep = '')) %>%
  dplyr::select(`Protein(s)`, Sequence, `Modified Sequence`, `Spectrum File`, `Spectrum Scan Number`) %>%
  mutate(fraction = sub('.*hph_(.*)\\.raw\\.mgf$', '\\1', `Spectrum File`)) %>%
  mutate(accession = `Protein(s)`, scan = `Spectrum Scan Number`, sequence = Sequence, modSequence = `Modified Sequence`) %>%
  left_join(proteinAnnotation) %>%
  dplyr::select(fraction, scan, accession, gene, detectablePeptides, sequence, modSequence)

New names:
* `` -> ...1
Rows: 157,560
Columns: 24
Delimiter: "\t"
chr [16]: Protein(s), Sequence, AAs Before, AAs After, Position, Modified Sequence, Variable Modifications, Fixed Modifications, Spectrum File, Spectrum Title, ...
dbl [ 8]: , Spectrum Scan Number, RT, m/z, Theoretical Mass, Isotope Number, Precursor m/z Error [ppm], Confidence [%]

Use `spec()` to retrieve the guessed column specification
Pass a specification to the `col_types` argument to quiet this message
Joining, by = "accession"

Now we need the quantification data. These are hosted in individual files in the same directory as the above data and were generated as part of the raw data processing described above. Again, if you have not followed the same initial data processing pipeline as me, your files may look different.

##########################################################################################
##get the list of files to be processed
quantFiles = as.list(list.files(paste(baseRepository, '/quantFiles/', sep = ''),
           pattern = '_Matrix.txt', full.names = TRUE))
##process using our defined function, combineQuantFiles
quantDataSet = lapply(quantFiles, combineQuantFiles)

[1] "Processing file for fraction 1."
[1] "Processing file for fraction 10."
[1] "Processing file for fraction 11."
[1] "Processing file for fraction 12."
[1] "Processing file for fraction 2."
[1] "Processing file for fraction 3."
[1] "Processing file for fraction 4."
[1] "Processing file for fraction 5."
[1] "Processing file for fraction 6."
[1] "Processing file for fraction 7."
[1] "Processing file for fraction 8."
[1] "Processing file for fraction 9."

Now combine the quant data into a single data frame and combine with the previously processed PSM data.

##########################################################################################
##collapse the data into a single data frame
allQuantData = do.call('rbind', quantDataSet)
##combine with the PSM data from above
psmQuant = psm %>%
  left_join(allQuantData)

Joining, by = c("fraction", "scan")

Now we will filter the expression data to discard any entries where the signal is just too low to be confident in the data. I generally will use a cutoff of a minimum signal of 10 per TMT channel. In this case we have 9 channels, so we will use a sum signal of 90 as a cutoff. Here, we also filter out genes that have no assigned gene name as this will cause problems with DEqMS later on. If you don’t want to do this, I suggest using ‘accession’ instead of ‘gene’ in the DEqMS analysis below.

##########################################################################################
psmQuant$sampleSignal = rowSums(psmQuant[,c(8:16)])
psmQuantFiltered = subset(psmQuant, psmQuant$sampleSignal >= 90 & !is.na(psmQuant$gene) & !grepl('-', psmQuant$gene))
psmQuantFiltered

We will also normalize the data to deal with any large loading differences. We save both the non-normalized data here as well and check the normalization with a plot. In the code below, we only include columns where we have sample data. In this case we didn’t use all 11 TMT reporter ion channels, so we get rid of those where we didn’t have a sample.

##########################################################################################
quantInput = psmQuantFiltered[,c(1:16)] ##make a new data frame object
saveRDS(quantInput, paste(baseRepository, '/dataset_peptideQuantDataPreDEqMS.rds', sep = '')) ##save the data
quantInputLog = quantInput[,c(7,4,8:16)] ##choose only columns where we have samples
quantInputLog[,3:11][quantInputLog[,3:11] == 0] = NA ##set zero values to NA
quantInputLog[,3:11] = log2(quantInputLog[,3:11])
quantInputNormalized = medianSweeping(quantInputLog, group_col = 2)
boxplot(quantInputNormalized, las = 2, ylab = 'log2 ratio', main = 'normalized TMT data')

This plot looks good, Our data are generally mean centered with the middle three samples have some more variation (these are actually treated samples, so this is not surprising). Lets continue with the DEqMS analysis. The next thing we need to do for DeQMS analysis is make our sample table. Our sample layout is ‘A1’,‘A2’,‘A3’,‘B1’,‘B2’,‘B3’,‘C1’,‘C2’,‘C3’, and a ‘pool’ of all samples in the last channel in a TMT10-plex format (TMT131C is not used).

##########################################################################################
cond = as.factor(c('A','A','A','B','B','B','C','C','C'))
design = model.matrix(~0+cond) 
colnames(design) = gsub("cond","",colnames(design))

Now we can build the Limma model and get the different contrasts. Again, as noted above, we are just following the DEqMS pipeline as described on their own page. The limma package also has some great documentation that explains many of the processes that DEqMS is using below.

##########################################################################################
geneMatrix = as.matrix(quantInputNormalized) ##extract just the quantification values to a matrix
limmaFit1 = lmFit(geneMatrix, design) ##generate the initial model fit using limma

Partial NA coefficients for 2 probe(s)

limmaContrasts = c('A-B','A-C','B-C') ##the different conditions we want to compare
limmaContrastDesign =  makeContrasts(contrasts = limmaContrasts, levels = design) ##make a new design matrix with our contrasts
limmaFit2 = eBayes(contrasts.fit(limmaFit1, contrasts = limmaContrastDesign)) ##perform eBayes smoothing using limma

Lastly, we will use DEqMS to summarize the data and extract out the different contrasts to save them to a file for each.

##########################################################################################
psmCountTable = as.data.frame(table(quantInput$gene)) ##count the number of times each gene is appearing
rownames(psmCountTable) = psmCountTable$Var1 ##assign genes as the row names
limmaFit2$count = psmCountTable[rownames(limmaFit2$coefficients),2] ##add the counts to the limma model
limmaFit3 = spectraCounteBayes(limmaFit2) ##do the final processing with deqms
head(limmaFit3$coefficients) ##take a look at the final data

             A-B         A-C        B-C
A2M    1.7989850  0.16655554 -1.6324295
A2ML1  1.5225197  0.44212828 -1.0803914
AAAS   1.5744897  0.25890516 -1.3155845
AACS   0.4925278 -0.04314785 -0.5356756
AAGAB  0.6976049 -0.78425106 -1.4818560
AAK1  -0.6994632  0.17290417  0.8723674

##save the individual contrasts
deqmsResults = outputResult(limmaFit3, coef_col = 1) 
write.table(deqmsResults, 
            paste(baseRepository, '/dataset_deqmsA-B.csv', sep = ''),
            sep = ',', row.names = FALSE, quote = FALSE)
##
deqmsResults = outputResult(limmaFit3, coef_col = 2) 
write.table(deqmsResults, 
            paste(baseRepository, '/dataset_deqmsA-C.csv', sep = ''),
            sep = ',', row.names = FALSE, quote = FALSE)
##
deqmsResults = outputResult(limmaFit3, coef_col = 3) 
write.table(deqmsResults, 
            paste(baseRepository, '/dataset_deqmsB-C.csv', sep = ''),
            sep = ',', row.names = FALSE, quote = FALSE)

Generate the variance plot to see the DEqMS model for the data.

##########################################################################################
VarianceBoxplot(limmaFit3, n=20, xlab = 'PSM count', main = 'DEqMS analysis of TMT data')

So we can see variance increasing as the PSM count gets lower, which is what we expect. Now let us make a volcano plot to show the data, as an example.

##########################################################################################
deqmsResults = outputResult(limmaFit3, coef_col = 1)
deqmsResults$plotPValue = -log10(deqmsResults$sca.P.Value)
deqmsResults$pColors = ifelse(deqmsResults$logFC <= -1 & deqmsResults$plotPValue >= 3, brewer.pal(8,'RdBu')[8],
                              ifelse(deqmsResults$logFC >= 1 & deqmsResults$plotPValue >= 3, brewer.pal(8,'RdBu')[1],
                                     brewer.pal(8,'Greys')[6]))
##
ggplot(deqmsResults, aes(logFC, plotPValue)) +
  geom_point(size = 1, color = deqmsResults$pColors, alpha = 0.5) +
  labs(x = 'log2(A fold change to B)', y = '-log10(P-value)', title = 'Sample A versus B') +
  scale_y_continuous(limits = c(0,12), breaks = seq(0,18,3)) +
  geom_vline(xintercept = c(-1,1), linetype = 'dashed') +
  geom_hline(yintercept = 3, linetype = 'dashed') +
  theme_classic()
ggsave(paste(baseRepository, '/scatter_deqmsA-B.pdf', sep = ''),
       height = 4, width = 4, useDingbats = FALSE)

So quite a bit has changed here. If you want to plot other contrasts, simply change the ‘coef_col’ value in the first command in the code above to a different value. You can now investigate these data further to determine biological variance seen in your own data!

Session info

sessionInfo()

R version 3.6.3 (2020-02-29)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 18363)

Matrix products: default

locale:
[1] LC_COLLATE=English_Canada.1252  LC_CTYPE=English_Canada.1252    LC_MONETARY=English_Canada.1252 LC_NUMERIC=C                   
[5] LC_TIME=English_Canada.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] DEqMS_1.4.0        limma_3.42.2       vroom_1.2.1        RColorBrewer_1.1-2 forcats_0.5.0      stringr_1.4.0      dplyr_0.8.5       
 [8] purrr_0.3.4        readr_1.3.1        tidyr_1.0.3        tibble_3.0.1       ggplot2_3.3.1      tidyverse_1.3.0   

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.4.6       lubridate_1.7.8    lattice_0.20-38    assertthat_0.2.1   digest_0.6.25      R6_2.4.1           cellranger_1.1.0  
 [8] plyr_1.8.6         backports_1.1.7    reprex_0.3.0       evaluate_0.14      httr_1.4.1         pillar_1.4.4       rlang_0.4.6       
[15] readxl_1.3.1       rstudioapi_0.11    blob_1.2.1         rmarkdown_2.1      labeling_0.3       bit_1.1-15.2       munsell_0.5.0     
[22] broom_0.5.6        compiler_3.6.3     modelr_0.1.8       xfun_0.14          pkgconfig_2.0.3    base64enc_0.1-3    htmltools_0.4.0   
[29] tidyselect_1.1.0   matrixStats_0.56.0 fansi_0.4.1        crayon_1.3.4       dbplyr_1.4.4       withr_2.2.0        grid_3.6.3        
[36] nlme_3.1-144       jsonlite_1.6.1     gtable_0.3.0       lifecycle_0.2.0    DBI_1.1.0          magrittr_1.5       scales_1.1.1      
[43] cli_2.0.2          stringi_1.4.6      farver_2.0.3       fs_1.4.1           xml2_1.3.2         ellipsis_0.3.1     generics_0.0.2    
[50] vctrs_0.3.0        tools_3.6.3        bit64_0.9-7        glue_1.4.1         hms_0.5.3          parallel_3.6.3     yaml_2.2.1        
[57] colorspace_1.4-1   rvest_0.3.5        knitr_1.28         haven_2.2.0