--- title: "qc-tidying" author: - name: Oliver M. Crook package: hdxmsqc output: BiocStyle::html_document: toc_float: yes abstract: "This vignette describes how to pefrom quality control for mass-spectrometry based hydrogen deuterium exchange experiment. \n" vignette: | %\VignetteIndexEntry{Qualityt control for differential hydrogen deuterium exchange mass spectrometry data} %\VignetteEngine{knitr::rmarkdown} %\VignetteKeywords{Mass Spectrometry, MS, MSMS, Proteomics, Metabolomics, Infrastructure, Quantitative} %\VignetteEncoding{UTF-8} --- # Introduction The `hdxmsqc` package is a quality control assessment package from hydrogen-deuterium exchange mass-spectrometry (HDX-MS) data. The functions look for outliers in retention time and ion mobility. They also examine missing values, mass errors, intensity based outliers, deviations of the data from monotonicity, the correlation of charge states, whether uptake values are coherent based on overlapping peptides and finally the similarity of the observed to the theoretical spectra observed. This package is designed to help those performing iterative quality control through manual inspection but also a set of metric and visualizations by which practitioners can use to demonstrate they have high quality data. # packages The packages required are the following. ```{r,} suppressMessages(require(hdxmsqc)) require(S4Vectors) suppressMessages(require(dplyr)) require(tidyr) require(QFeatures) require(RColorBrewer) require(ggplot2) require(MASS) require(pheatmap) require(Spectra) require(patchwork) ``` # Data We first load the data, as exported from HDExaminer. ```{r,} BRD4uncurated <- data.frame(read.csv(system.file("extdata", "ELN55049_AllResultsTables_Uncurated.csv", package = "hdxmsqc", mustWork = TRUE))) ``` The following code chunk tidies dataset, which improves the formatting and converts to wide format. It will also note the number of states, timepoints and peptides. ```{r,} BRD4uncurated_wide <- processHDE(HDExaminerFile = BRD4uncurated, proteinStates = c("wt", "iBET")) ``` The next code chunk extracts the columns with the quantitative data. ```{r,} i <- grep(pattern = "X..Deut", x = names(BRD4uncurated_wide)) ``` We now parse the object into an object of class `Qfeatures`. This standardises the formatting of the data. ```{r,} BRD4df <- readQFeatures(assayData = BRD4uncurated_wide, ecol = i, names = "Deuteration", fnames = "fnames") ``` # Visualisation A simple heatmap of our data can give us a sense of it. ```{r,} pheatmap(assay(BRD4df), cluster_rows = FALSE, scale = "row") ``` # Examining missing values Here, we can plot where the missing values are: ```{r,} plotMissing(object = BRD4df) ``` Here, we can filter data that is not missing at random: ```{r,} BRD4df_filtered <- isMissingAtRandom(object = BRD4df) ``` We can then replot missing-ness: ```{r,} plotMissing(object = BRD4df_filtered) ``` The values that are missing are all at the zero time-points where deuterium uptake should be 0, we can simply impute these values. ```{r,} BRD4df_filtered_imputed <- impute(BRD4df_filtered, method = "zero", i = 1) ``` # Empirical vs Theoretical errors ```{r,} massError <- computeMassError(object = BRD4df_filtered_imputed) plotMassError(object = BRD4df_filtered_imputed) ``` # Intensity based outlier detection Using linear-model based outlier detection we see whether there are Spectra that have variable intensity based on their mean intensity. A linear model is fitted to the log-mean and log-variance of the intensities. These should follow a linear trend. Cook's distance is used to determine outliers are consider if their distance is greater than 2/$\sqrt(n)$, where $n$ is the number of peptides. ```{r,} intensityOutlier <- intensityOutliers(object = BRD4df_filtered_imputed) plotIntensityOutliers(object = BRD4df_filtered_imputed) ``` # Retention time analysis Retention time outlier detection looks at the usual variability of retention time search window and the actual left/right windows of the retention time. Outliers are flagged if their retention time falls outside 1.5 * interquartile range. ```{r,} dfrt <- rTimeOutliers(object = BRD4df_filtered_imputed) plotrTimeOutliers(object = BRD4df_filtered_imputed) ``` # Monotonicity statistics This uses a statistic to detect differences from monotonic behavior. First, we need to specify the experimental design and the timepoints used. ```{r,} experiment <- c("wt", "iBET") timepoints <- rep(c(0, 15, 60, 600, 3600, 14000), each = 3) ``` The monotonicity statistic measure the deviation from monotoncity. Whilst some deviation is expected from random fluctuations, it is worth double checking those that are strong deviates compare to the rest of the data. ```{r,} monoStat <- computeMonotoneStats(object = BRD4df_filtered_imputed, experiment = experiment, timepoints = timepoints) out <- plotMonotoneStat(object = BRD4df_filtered_imputed, experiment = experiment, timepoints = timepoints) out ``` # Ion Mobility Time analysis In a similar analysis to the retention time analysis, for ion mobility time we can also see whether there are random deviation in the ion mobility windows. Again, we define outliers that deviate outside the typical 1.5 * IQR. ```{r,} imTimeOut <- imTimeOutlier(object = BRD4df_filtered_imputed) plotImTimeOutlier(object = BRD4df_filtered_imputed) ``` # Charge state correlation We check that charge states are correlated. Whilst we don't expect exactly the same before - low correlation maybe concerning. ```{r,} csCor <- chargeCorrelationHdx(object = BRD4df_filtered_imputed, experiment = experiment, timepoints = timepoints) csCor ``` # Using replicates to determine outliers and variability ```{r,} replicateVar <- replicateCorrelation(object = BRD4df_filtered_imputed, experiment = experiment, timepoints = timepoints) replicateOut <- replicateOutlier(object = BRD4df_filtered_imputed, experiment = experiment, timepoints = timepoints) ``` # Using sequence overlap information are uptake values compatible We can also check whether uptakes are compatible with overlapping peptides. The difference in uptake cannot be more different than the difference in the number of exchangeable amides. The default methodology only checks whether sequence with up-to 5 different exchangeable amides are compatible to keep run-times lower. Larger difference may indicate different chemical changes or back-exchange properties. ```{r,} tocheck <- compatibleUptake(object = BRD4df_filtered_imputed, experiment = experiment, timepoints = timepoints) ``` # Comparison of Spectra In this section, we can directly examine the differences between the theoretical spectra one would expect from the computed deuterium uptake and the actual observed spectra. Deviations observed in the spectra could suggest contamination, false identifications or poor quality spectra. A score is generated using the cosine similarity between the spectra - which is equivalent to the normalized dot product. The spectra pairs can be also be visualized. Load in some Spectra from HDsite which should match those of HDExaminer ```{r,} hdxsite <- data.frame(read.csv(system.file("extdata", "BRD4_RowChecked_20220628_HDsite.csv", package = "hdxmsqc", mustWork = TRUE), header = TRUE, fileEncoding = 'UTF-8-BOM')) BRD4matched <- read.csv(system.file("extdata", "BRD4_RowChecked_20220628_HDE.csv", package = "hdxmsqc", mustWork = TRUE), header = TRUE, fileEncoding = 'UTF-8-BOM') ``` ```{r,} spectraCompare <- spectraSimilarity(peaks = hdxsite, object = BRD4matched, experiment = experiment, numSpectra = NULL) ``` The scores can be accesses as follows: ```{r,} head(spectraCompare$observedSpectra$score) ``` To visualise these spectra we can use the following function ```{r,} plotSpectraMirror(spectraCompare$observedSpectra[1,], spectraCompare$matchedSpectra[1,], ppm = 300) ``` Finally, a summarise quality control table can be produced and saved in a .csv file if desired. ```{r,} qctable <- qualityControl(object = BRD4df_filtered_imputed, massError = massError, intensityOutlier = intensityOutlier, retentionOutlier = dfrt, monotonicityStat = monoStat, mobilityOutlier = imTimeOut, chargeCorrelation = csCor, replicateCorrelation = replicateVar, replicateOutlier = replicateOut, sequenceCheck = tocheck, spectraCheck = spectraCompare, experiment = experiment, timepoints = timepoints ) ``` ```{r,} sessionInfo() ```