standPep: Standardize peptide results

In qzhang503/proteoQ: Processing and Informatic Analysis of Mass Spectrometrirc Data

Standardize peptide results

Description

standPep standardizes peptide results from mergePep with additional, stand-alone choices in data alignment. The utility is typically applied after the assembly of peptide data via mergePep. It further supports iterative normalization against data under selected sample columns, data rows or both.

Usage

standPep(
  method_align = c("MC", "MGKernel"),
  col_select = NULL,
  range_log2r = c(10, 90),
  range_int = c(5, 95),
  n_comp = NULL,
  seed = NULL,
  plot_log2FC_cv = FALSE,
  ...
)

Arguments

method_align	Character string indicating the method in aligning log2FC across samples. MC: median-centering; MGKernel: the kernel density defined by multiple Gaussian functions (normalmixEM). At the MC default, the ratio profiles of each sample will be aligned in that the medians of the log2FC are zero. At MGKernel, the ratio profiles of each sample will be aligned in that the log2FC at the maximums of kernel density are zero.
col_select	Character string to a column key in expt_smry.xlsx. At the NULL default, the column key of Select in expt_smry.xlsx will be used. In the case of no samples being specified under Select, the column key of Sample_ID will be used. The non-empty entries under the ascribing column will be used in indicated analysis.
range_log2r	Numeric vector at length two. The argument specifies the range of the log2FC for use in the scaling normalization of standard deviation across samples. The default is between the 10th and the 90th quantiles.
range_int	Numeric vector at length two. The argument specifies the range of the intensity of reporter ions (including I000) for use in the scaling normalization of standard deviation across samples. The default is between the 5th and the 95th quantiles.
n_comp	Integer; the number of Gaussian components to be used with method_align = MGKernel. A typical value is 2 or 3. The variable n_comp overwrites the argument k in normalmixEM.
seed	Integer; a seed for reproducible fitting at method_align = MGKernel.
plot_log2FC_cv	Logical; if TRUE, the distributions of the CV of peptide log2FC will be plotted. The default is TRUE.
...	slice_: variable argument statements for the identification of row subsets. The partial data will be taken for parameterizing the alignment of log2FC across samples. The full data set will be updated subsequently with the newly derived parameters. Note that there is no data entry removals from the complete data set with the slice_ procedure. The variable argument statements should be in the following format: each of the statement contains a list of logical expression(s). The lhs needs to start with slice_. The logical condition(s) at the rhs needs to be enclosed in exprs with round parenthesis. For example, pep_len is a column key present in Peptide.txt. The slice_peps_at = exprs(pep_len >= 10, pep_len <= 50) will extract peptide entries with the number of amino acid residues betwen 10 and 50 for log2FC alignment. Shorter or longer peptide sequences will remain in Peptide.txt but not used in the parameterization. See also normPSM for the variable arguments of filter_. Additional parameters from normalmixEM, i.e., maxit, the maximum number of iterations allowed; epsilon, tolerance limit for declaring algorithm convergence.

Details

In the primary output file, "Peptide.txt", values under columns log2_R... are logarithmic ratios at base 2 in relative to the reference(s) within each multiplex TMT set, or to the row means within each plex if no reference(s) are present. Values under columns N_log2_R... are aligned log2_R... according to method_align without scaling normalization. Values under columns Z_log2_R... are N_log2_R... with additional scaling normalization. Values under columns I... are reporter-ion or LFQ intensity before normalization. Values under columns N_I... are normalized I.... Values under columns sd_log2_R... are the standard deviation of the log2FC of proteins from ascribing peptides.

In general, median statistics is applied when summarizing numeric peptide data from different LCMS series. One exception is pep_expect where geometric mean is used.

Description of the column keys in the inputs and outputs:
system.file("extdata", "peptide_keys.txt", package = "proteoQ")

Value

The primary output is in .../Peptide/Peptide.txt.

See Also

Metadata
load_expts for metadata preparation and a reduced working example in data normalization

Data normalization
normPSM for extended examples in PSM data normalization
PSM2Pep for extended examples in PSM to peptide summarization
mergePep for extended examples in peptide data merging
standPep for extended examples in peptide data normalization
Pep2Prn for extended examples in peptide to protein summarization
standPrn for extended examples in protein data normalization.
purgePSM and purgePep for extended examples in data purging
pepHist and prnHist for extended examples in histogram visualization.
extract_raws and extract_psm_raws for extracting MS file names

Variable arguments of 'filter_...'
contain_str, contain_chars_in, not_contain_str, not_contain_chars_in, start_with_str, end_with_str, start_with_chars_in and ends_with_chars_in for data subsetting by character strings

Missing values
pepImp and prnImp for missing value imputation

Informatics
pepSig and prnSig for significance tests
pepVol and prnVol for volcano plot visualization
prnGSPA for gene set enrichment analysis by protein significance pVals
gspaMap for mapping GSPA to volcano plot visualization
prnGSPAHM for heat map and network visualization of GSPA results
prnGSVA for gene set variance analysis
prnGSEA for data preparation for online GSEA.
pepMDS and prnMDS for MDS visualization
pepPCA and prnPCA for PCA visualization
pepLDA and prnLDA for LDA visualization
pepHM and prnHM for heat map visualization
pepCorr_logFC, prnCorr_logFC, pepCorr_logInt and prnCorr_logInt for correlation plots
anal_prnTrend and plot_prnTrend for trend analysis and visualization
anal_pepNMF, anal_prnNMF, plot_pepNMFCon, plot_prnNMFCon, plot_pepNMFCoef, plot_prnNMFCoef and plot_metaNMF for NMF analysis and visualization

Custom databases
Uni2Entrez for lookups between UniProt accessions and Entrez IDs
Ref2Entrez for lookups among RefSeq accessions, gene names and Entrez IDs
prepGO for gene ontology
prepMSig for molecular signatures
prepString and anal_prnString for STRING-DB

Column keys in PSM, peptide and protein outputs
system.file("extdata", "psm_keys.txt", package = "proteoQ")
system.file("extdata", "peptide_keys.txt", package = "proteoQ")
system.file("extdata", "protein_keys.txt", package = "proteoQ")

Examples

# ===================================
# Peptide normalization
# ===================================

## !!!require the brief working example in `?load_expts`

# ===================================
# (1) `MGKernel`
# ===================================

# !!! Initial `Peptide.txt` results from `mergePep()` are in median centering.
# !!! The first `MGKernel` normalization will always be applied to all samples.
# !!! If changing `n_comp`, succeeding `MGKernel` normalization(s) will be 
#     applied to all samples (a fresh-start principle, see also section 4)

# (1.1) the first `MGKernel` after `mergePep()`
# (default double trimming by log2FC and intensity percentiles also applied)
standPep(
  method_align = MGKernel, 
  n_comp = 3, 
  seed = 749662, 
  maxit = 200, 
  epsilon = 1e-05, 
)

# histograms for all data rows
# (for samplicity only plot against samples indicated under `expt_smry.xlsx::BI_1`)
pepHist(scale_log2r = TRUE, col_select = BI_1, filename = bi1.png)

# human subset
# (a little off in relative to all data rows)
pepHist(scale_log2r = TRUE, col_select = BI_1, filter_by = exprs(species == "human"), filename = bi1_human.png)

# mouse subset 
# (a lot off in relative to all data rows)
pepHist(scale_log2r = TRUE, col_select = BI_1, filter_by = exprs(species == "mouse"), filename = bi1_mouse.png)

# (1.2) additive step to (1.1): against `expt_smry.xlsx::W2` samples and based on their human subset
standPep(
  method_align = MGKernel, 
  n_comp = 3, 
  seed = 749662, 
  maxit = 200, 
  epsilon = 1e-05, 
  col_select = W2,  
  slice_peps_by = exprs(species == "human"),
)

# human subset 
# (W2 samples are now aligned)
pepHist(scale_log2r = TRUE, col_select = BI_1, filter_by = exprs(species == "human"), filename = bi1_human_slicehuw2.png)

# (1.3) additive to (1.2): against `expt_smry.xlsx::W16` samples and based on their human subset
standPep(
  method_align = MGKernel, 
  n_comp = 3, 
  seed = 749662, 
  maxit = 200, 
  epsilon = 1e-05, 
  col_select = W16,  
  slice_peps_by = exprs(species == "human"),
)

# human subset 
# (W16 samples are now also aligned)
pepHist(scale_log2r = TRUE, col_select = BI_1, filter_by = exprs(species == "human"), filename = bi1_human_slicehuw16.png)

# a side effect: to recapitulate the misalignment between human data and human + mouse data
# (this is because density curves are based on the latest `standPep` at `method_align = MGKernel`)
pepHist(scale_log2r = TRUE, col_select = BI_1, filename = bi1_recap.png)


# ===================================
# (2) Mixed-bed
# ===================================
# start over
unlink(file.path(dat_dir, "Peptide"), recursive = TRUE, force = TRUE)
PSM2Pep()
mergePep()

# data in initial `Peptide.txt` from `mergePep()` are aligned by MC
pepHist(scale_log2r = TRUE, col_select = BI_1, filename = mc.png)

# (2.1) the first `MGKernel` alignment against all samples using all data rows
standPep(
  method_align = MGKernel, 
  n_comp = 3, 
  seed = 883, 
  maxit = 200, 
  epsilon = 1e-05, 
)

pepHist(scale_log2r = TRUE, col_select = BI_1, filename = mG.png)

# (2.2) followed by MC against selected samples using all data rows
# (arguments `n_comp`, `seed`, `maxit`, `epsilon` have no effects at `method_align = MC`)
standPep(
  method_align = MC,
  
  n_comp = 3,
  seed = 883,
  maxit = 200,
  epsilon = 1e-05,
  
  col_select = Select_sub, 
)

# the net result is mixed-bed alignment of MGKernel and MC
# (with samples indicated by `Select_sub` aligned by MC and the remaining by MGKernel)
pepHist(scale_log2r = TRUE, col_select = BI_1, filename = mix.png)


# ===================================
# (3) Mixed-bed against data subset
# ===================================
# start over
unlink(file.path(dat_dir, "Peptide"), recursive = TRUE, force = TRUE)
PSM2Pep()
mergePep()

pepHist(scale_log2r = TRUE, col_select = BI_1, filename = mc.png)

# (3.1) MC alignment for all samples, but only selected data rows used for normalization
standPep(
  method_align = MC, 
  slice_peps_by = exprs(prot_n_psm >= 10),
)

pepHist(scale_log2r = TRUE, col_select = BI_1, filename = mc_selrows.png)

# (3.2) the first `MGKernel` for all samples using selected data rows
standPep(
  method_align = MGKernel, 
  n_comp = 3, 
  seed = 400, 
  maxit = 200, 
  epsilon = 1e-05,
  
  # will be forced to all samples since this is the first `MGKernel`
  # col_select = Select_sub, 

  slice_peps_by = exprs(prot_n_psm >= 10),
)

pepHist(scale_log2r = TRUE, col_select = BI_1, filename = mG_selrows.png)

# (3.3) back to MC for selected samples using selected data rows
# (mixed-bed again, but based on data subset by `prot_n_psm >= 10`)
standPep(
  method_align = MC, 
  col_select = Select_sub, 
  slice_peps_by = exprs(prot_n_psm >= 10),
)

# A side effect in comparing `MC` and `MGKernel`
# (density curves are from the preceding `MGKernel` in (3.2))
pepHist(scale_log2r = TRUE, col_select = BI_1, filename = mix_selcols_selrows.png)


# ===================================
# (4) Modified `n_comp`
# ===================================
# start over
unlink(file.path(dat_dir, "Peptide"), recursive = TRUE, force = TRUE)
PSM2Pep()
mergePep()

pepHist(scale_log2r = TRUE, col_select = BI_1, filename = mc.png)

# (4.1) first `MGKernel` at `n_comp = 3` for all samples
standPep(
  method_align = MGKernel, 
  n_comp = 3, 
  seed = 400, 
  maxit = 200, 
  epsilon = 1e-05, 
)

pepHist(scale_log2r = TRUE, col_select = BI_1, filename = mG3.png)

# (4.2) a fresh start since changing `n_comp`
# (e.g. `col_select = Select_sub` ignored; instead apply `MGKernel` to all samples)
standPep(
  method_align = MGKernel, 
  n_comp = 2, 
  seed = 400, 
  maxit = 200, 
  epsilon = 1e-05,
  
  # ignored
  col_select = Select_sub, 
)

pepHist(scale_log2r = TRUE, col_select = BI_1, filename = mG2.png)


# ===================================
# (5) housekeeping normalizers: 
#     suggest `method_align = MC`
# ===================================
# start over
unlink(file.path(dat_dir, "Peptide"), recursive = TRUE, force = TRUE)
PSM2Pep()
mergePep()

# initial `Peptide.txt` from `mergePep()` is aligned by MC
pepHist(scale_log2r = TRUE, col_select = BI_1, filename = mc.png)

# (GAPDH more abundant in W16 according to MC alignment)
pepHist(scale_log2r = TRUE, col_select = BI_1, filter_ = exprs(gene == "GAPDH"), filename = mcGAPDH.png)

# heat map for `BI_1` samples
# (outputs under `Peptide/Heatmap` folder; for help, ?pepHM)
pepHM(
  col_select = BI_1, 
  xmin = -2,
  xmax = 2,
  xmargin = 0.1,
  annot_cols = c("Group", "Color", "Alpha", "Shape"),
  annot_colnames = c("Group", "Lab", "Batch", "WHIM"),
  cluster_rows = TRUE,
  annot_rows = c("gene"),
  show_rownames = TRUE,
  show_colnames = TRUE,
  fontsize_row = 10,
  cellwidth = 12,
  cellheight = 12,
  width = 18,
  height = 12,
  
  filter_by = exprs(gene %in% c("GAPDH")),
  filename = "mcGAPDH.png",
)

# heat map for all samples
# (no GAPDH under `JHU_TMT1` and `PNNL_TMT1`; will be problematic in data alignment next)
pepHM(
  # col_select = BI_1, 
  xmin = -2,
  xmax = 2,
  xmargin = 0.1,
  annot_cols = c("Group", "Color", "Alpha", "Shape"),
  annot_colnames = c("Group", "Lab", "Batch", "WHIM"),
  cluster_rows = TRUE,
  annot_rows = c("gene"),
  show_rownames = TRUE,
  show_colnames = TRUE,
  fontsize_row = 10,
  cellwidth = 12,
  cellheight = 12,
  width = 18,
  height = 12,
  
  filter_by = exprs(gene %in% c("GAPDH")),
  filename = "mcGAPDH_all.png",
)

# first renormalize against GAPDH for all samples
# (not to use `MGkernel` as may not have enough data entries from housekeeper(s))
standPep(
  method_align = MC, 
  slice_hskp = exprs(gene %in% c("GAPDH")),
)

# (now log2FC profiles aligned by GAPDH for `BI_1` samples)
pepHist(scale_log2r = TRUE, col_select = BI_1, filename = wrong_hskp.png)

# (no histograms for `JHU_TMT1` and `PNNL_TMT1` samples as no GAPDH underneath)
pepHist(scale_log2r = TRUE, filename = wrong_hskp_all.png)

# not to keep the above example with no data under `JHU_TMT1` and `PNNL_TMT1`
# (surely need different normalizer(s))
unlink(file.path(dat_dir, "Peptide"), recursive = TRUE, force = TRUE)


## Not run: 
PSM2Pep()
mergePep()

# change to `MGKernel`
standPep(
  method_align = MGKernel, 
  n_comp = 3, 
  seed = 400, 
  maxit = 200, 
  epsilon = 1e-05, 
)

pepHist(scale_log2r = TRUE, col_select = BI_1, filename = mG.png)

# then renormalize against data from GAPDH
# (error: too few entries for fitting with multiple Gaussians)
standPep(
  method_align = MGKernel, 
  n_comp = 3, 
  seed = 400, 
  maxit = 200, 
  epsilon = 1e-05, 
  
  col_select = W2, 
  slice_hskp = exprs(gene %in% c("GAPDH")),
)
## End(Not run)

qzhang503/proteoQ documentation built on Dec. 14, 2024, 12:27 p.m.