The HERONSequenceData object

The HERONSequenceDataSet object is the first object needed for running HERON. The HERONSequenceDataSet is a subclass of SummarizedExperiment from base Bioconductor. The components are a sequence data.frame or matrix where rows are amino acid sequences of the peptide binding probe and the columns are samples, a colData DFrame, and a metadata data.frame that maps sequences to probes.

data(heffron2021_wuhan)
knitr::kable(assay(heffron2021_wuhan)[1:5,1:5])

knitr::kable(head(colData(heffron2021_wuhan)))

probe_meta <- metadata(heffron2021_wuhan)$probe_meta
knitr::kable(head(probe_meta[,c("PROBE_SEQUENCE", "PROBE_ID")]))

Creating the HERONSequenceData object

The heffron2021_wuhan data object is already a HERONSequenceDataSet object. We can create the same object using the HERONSequenceDataSet constructor as an example.

seq_mat <- assay(heffron2021_wuhan)
sds <- HERONSequenceDataSet(seq_mat)
colData(sds) <- colData(heffron2021_wuhan)
metadata(sds)$probe_meta <- probe_meta

Pre-process data

The data is pre-processed by first quantile normalizing on sequence-level data.

seq_ds_qn <- quantileNormalize(sds)
knitr::kable(head(assay(seq_ds_qn)[,1:5]))

Calculate probe-level p-values

The quantile normalized data is then used with the colData and probe_meta data frames to calculate probe-level p-values.

The calcCombPValues call will first calculate the probe level p-values using the colData, adjust the probe p-values on a column-by-column basis frame.

The combined p-values can be calculated on the sequence data set, then converted to a probe data set through convertSequenceDSToProbeDS, or calculated on the probe data set. If the data was smoothed prior using consecutive probes on the same protein, then the p-values should calculated on the smoothed probe data set.

seq_pval_res <- calcCombPValues(seq_ds_qn)
probe_pval_res <- convertSequenceDSToProbeDS(seq_pval_res)
knitr::kable(head(assay(probe_pval_res, "padj")[,1:5]))

Obtain probe-level calls

Once the p-values have been calculated, the makeProbeCalls will make the calls using the adjusted p-values calculated. makeProbeCalls also includes a filter to remove one-hit probes, which are probes that were called only in one sample and do not have a consecutive probe called in a single sample. The return value is HERONProbeDataSet with the calls included as an additional assay.

probe_calls_res <- makeProbeCalls(probe_pval_res)
knitr::kable(head(assay(probe_calls_res, "calls")[,1:5]))

K of N, Fraction, and Percent of calls made per sequence, probe, epitope, or protein can be determined using the getKofN function.

probe_calls_k_of_n <- getKofN(probe_calls_res)
knitr::kable(head(probe_calls_k_of_n))

Find Epitope Segments using the unique method

After calculating probe p-values and calls, the HERON can find epitope segments, i.e. blocks of consecutive probes that have been called within a sample or cluster of samples. The unique method finds all consecutive probes for each sample, then returns the unique set of all epitopic regions.

epi_segments_uniq_res <- findEpitopeSegments(
    probe_calls_res,
    segment_method = "unique"
)

knitr::kable(head(epi_segments_uniq_res))

The format of the epitope segments are SEQID_Begin_End, where the SEQID is the protein name, Begin is the starting position of the first probe within the consecutive block, and the End is the starting position of the last probe within the consecutive block.

Calculate Epitope-level p-values

With the epitope segments and the probe p-values, HERON can calculate a significance value for the segments using a meta p-value method. Here, we are calculating epitope p-values using Wilkinson’s max meta p-value method.

epi_padj_uniq <- calcEpitopePValues(
    probe_calls_res,
    epitope_ids = epi_segments_uniq_res,
    metap_method = "wmax1"
)

You can add sequence annotations to the row data of the HERONEpitopeDataSet object.

epi_padj_uniq <- addSequenceAnnotations(epi_padj_uniq)

Obtain Epitope-level calls

The makeEpitopeCalls method will work on the epitope adjusted p-values to make epitope-level sample. getKofN can also get the K of N results.

epi_calls_uniq <- makeEpitopeCalls(epi_padj_uniq, one_hit_filter = TRUE)
epi_calls_k_of_n_uniq <- getKofN(epi_calls_uniq)
knitr::kable(head(epi_calls_k_of_n_uniq))

Calculate Protein-level p-values

Calculate protein p-values using Tippett’s (Wilkinson’s Min) method.

prot_padj_uniq <- calcProteinPValues(
    epi_padj_uniq,
    metap_method = "tippetts"
)

Obtain Protein-level calls

prot_calls_uniq <- makeProteinCalls(prot_padj_uniq)
prot_calls_k_of_n_uniq <- getKofN(prot_calls_uniq)
knitr::kable(head(prot_calls_k_of_n_uniq))

Find Epitope Segments using the hclust method

epi_segments_hclust_res <- findEpitopeSegments(
    probe_calls_res,
    segment_method = "hclust",
    segment_score_type = "binary",
    segment_dist_method = "hamming",
    segment_cutoff = "silhouette"
)

epi_segments_hclust_res2 <- findEpitopeSegments(
    probe_calls_res,
    segment_method = "hclust",
    segment_score_type = "zscore",
    segment_dist_method = "euclidean",
    segment_cutoff = "silhouette"
)

Find Epitope Segments using the skater method

epi_segments_skater_res <- findEpitopeSegments(
    probe_calls_res,
    segment_method = "skater",
    segment_score_type = "binary",
    segment_dist_method = "hamming",
    segment_cutoff = "silhouette"
)

epi_segments_skater_res <- findEpitopeSegments(
    probe_calls_res,
    segment_method = "skater",
    segment_score_type = "zscore",
    segment_dist_method = "euclidean",
    segment_cutoff = "silhouette"
)

Other meta p-value methods

HERON also allows a choice for a number of meta p-value methods for combining p-values for epitopes (calcEpitopePValuesPDS) and proteins (calcProteinPValuesEDS). The methods supported by HERON are : min_bonf, min, max, fischer/sumlog, hmp/harmonicmeanp, wilkinsons_min1/tippets, wilkinsons_min2/wmin2, wilkinsons_min3, wilkinsons_min4, wilkinsons_min5, wilkinsons_max1/wmax1, wilkinsons_max2/wmax2, and cct.

When choosing a p-value method, keep in mind that the epitope p-value should be one that requires most of the probe p-values to be small (e.g. wmax1) and the protein should be one that requires at least one of the epitope p-values to be small (e.g. wmax1). Other p-value methods such as the cct and the hmp have been shown to be more accurate with p-value that have dependencies.

Making z-score cutoff calls

Some investigators would like to make z-score global level calls rather than using adjusted p-values. HERON is flexible to allow for such a setup. For example, the user wanted to make probe-level calls using a z-score cutoff > 2.

seq_pval_res_z <- calcCombPValues(
    obj = seq_ds_qn, 
    use = "z", 
    p_adjust_method = "none"
)


p_cutoff <- pnorm(2, lower.tail = FALSE)

probe_pval_res_z <- convertSequenceDSToProbeDS(seq_pval_res_z, probe_meta)

probe_calls_z2 <- makeProbeCalls(probe_pval_res_z, padj_cutoff = p_cutoff)
probe_k_of_n_z2 <- getKofN(probe_calls_z2)

knitr::kable(head(assay(probe_calls_z2,"calls")[,1:5]))

knitr::kable(head(probe_k_of_n_z2[probe_k_of_n_z2$K > 0,]))

The calls can also be used to find epitopes segments, p-values, and calls. Also, can be used to make protein level calls.

epi_segments_uniq_z2_res <- findEpitopeSegments(
    probe_calls_z2,
    segment_method = "unique"
)

If we want to find regions where the z > 2 for all of the consecutive probes, using the max meta p-value method will ensure that.

epi_pval_uniq_z2 <- calcEpitopePValues(
    probe_pds = probe_pval_res_z,
    epitope_ids = epi_segments_uniq_z2_res,
    metap_method = "max",
    p_adjust_method = "none"
)

Again, makeEpitopeCalls will be used in order to find significant regions.

epi_calls_uniq_z2 <- makeEpitopeCalls(
    epi_ds = epi_pval_uniq_z2, 
    padj_cutoff = p_cutoff
)

Other segmentation methods can be used with the called regions through the binary score clustering methods.

epi_segments_skater_z2_res <- findEpitopeSegments(
    probe_calls_z2,
    segment_method = "skater",
    segment_score_type = "binary",
    segment_dist_method = "hamming",
    segment_cutoff = "silhouette")

Smoothing across probes

The pepStat paper (https://pubmed.ncbi.nlm.nih.gov/23770318/) mentions that smoothing can sometimes help reduce the high-variability of the peptide array data. HERON has implemented a running average smoothing algorithm across protein probes that can handle missing values. The function smoothProbeMat will take as input a probe matrix and will return a matrix that has been smoothed. The smoothed matrix can then be processed through the workflow using the calcProbePValuesProbeMat call instead of the calcProbePValuesSeqMat call since now the probe signal is no longer a direct copy from the sequence.

probe_ds_qn <- convertSequenceDSToProbeDS(seq_ds_qn, probe_meta )

smooth_ds <- smoothProbeDS(probe_ds_qn)

After you smoothed the data using the probes, the probe p-value will be calculated on the probe-level rather than the sequence-level.

probe_sm_pval <- calcCombPValues(smooth_ds)

The probe calls can then be made as before.

probe_sm_calls <- makeProbeCalls(probe_sm_pval)
probe_sm_k_of_n <- getKofN(probe_sm_calls)
knitr::kable(assay(probe_sm_calls,"calls")[1:3,1:3])

knitr::kable(head(probe_sm_k_of_n[probe_sm_k_of_n$K > 0,]))

Calculate paired t-tests

Paired t-tests can be utilized and combined with the z-tests as well. Here is an example where we pair the five COVID- samples to the COVID+ samples and run calcProbePValuesSeqMat with the d_paired parameter set.

data(heffron2021_wuhan)
probe_meta <- attr(heffron2021_wuhan, "probe_meta")
colData_paired <- colData(heffron2021_wuhan)

## Make some samples paired by setting the sample ids
pre_idx <- which(colData_paired$visit == "pre")
colData_post <- colData_paired[colData_paired$visit == "post",]
new_ids <- colData_post$SampleName[seq_len(5)]
colData_paired$ptid[pre_idx[seq_len(5)]] = new_ids

paired_ds <- heffron2021_wuhan
colData(paired_ds) <- colData_paired

## calculate p-values
pval_res <- calcCombPValues(
    obj = paired_ds,
    d_paired = TRUE
)

knitr::kable(assay(pval_res[1:3,],"pvalue"))

Use the wilcox test for probe-level p-values

To account for non-normality, HERON implements the two sample wilcox test in placed of the t-test used for the differential p-values. Below is an example of how to use it.

seq_pval_res_w <- calcCombPValues(
    obj = seq_ds_qn, 
    use = "w", d_abs_shift = 1
)
#> exact:
probe_pval_res_w <- convertSequenceDSToProbeDS(seq_pval_res_w)
knitr::kable(assay(probe_pval_res_w[1:3, 1:5], "padj"))

Use of the condition column

You can use the condition column in the colData dataframe. It is used in the getKofN function. The data type of the column can be either a character or factor. When the condition column exists in the colData, three columns (number of samples (K), frequency (F), and percent (P)) are provided for each unique factor/character value.

col_data <- colData(heffron2021_wuhan)
covid <- which(col_data$visit == "post")

col_data$condition[covid[1:10]] <- "COVID2"

seq_ds <- heffron2021_wuhan
colData(seq_ds) <- col_data

seq_ds_qn <- quantileNormalize(seq_ds)
seq_pval_res <- calcCombPValues(seq_ds_qn)
probe_pval_res <- convertSequenceDSToProbeDS(seq_pval_res)
probe_calls_res <- makeProbeCalls(probe_pval_res)
probe_calls_k_of_n <- getKofN(probe_calls_res)
probe_calls_k_of_n <- 
    probe_calls_k_of_n[
        order(probe_calls_k_of_n$K, decreasing = TRUE),]
knitr::kable(head(probe_calls_k_of_n))