Uncovering the spatial landscape of molecular interactions within the tumor microenvironment through latent spaces.

Recent advances in spatial transcriptomics (STs) enable gene expression measurements from a tissue sample while retaining its spatial context. This technology enables unprecedented in situ resolution of the regulatory pathways that underlie the heterogeneity in the tumor as well as the tumor microenvironment (TME). The direct characterization of cellular co-localization with spatial technologies facilities quantification of the molecular changes resulting from direct cell-cell interaction, as it occurs in tumor-immune interactions. We present SpaceMarkers, a bioinformatics algorithm to infer molecular changes from cell-cell interactions from latent space analysis of ST data. We apply this approach to infer the molecular changes from tumor-immune interactions in Visium spatial transcriptomics data of metastasis, invasive and precursor lesions, and immunotherapy treatment. Further transfer learning in matched scRNA-seq data enabled further quantification of the specific cell types in which SpaceMarkers are enriched. Altogether, SpaceMarkers can identify the location and context-specific molecular interactions within the TME from ST data.

Investigation of the Role of PUFA Metabolism in Breast Cancer Using a Rank-Based Random Forest Algorithm.

Polyunsaturated fatty acid (PUFA) metabolism is currently a focus in cancer research due to PUFAs functioning as structural components of the membrane matrix, as fuel sources for energy production, and as sources of secondary messengers, so called oxylipins, important players of inflammatory processes. Although breast cancer (BC) is the leading cause of cancer death among women worldwide, no systematic study of PUFA metabolism as a system of interrelated processes in this disease has been carried out. Here, we implemented a Boruta-based feature selection algorithm to determine the list of most important PUFA metabolism genes altered in breast cancer tissues compared with in normal tissues. A rank-based Random Forest (RF) model was built on the selected gene list (33 genes) and applied to predict the cancer phenotype to ascertain the PUFA genes involved in cancerogenesis. It showed high-performance of dichotomic classification (balanced accuracy of 0.94, ROC AUC 0.99) We also retrieved a list of the important PUFA genes (46 genes) that differed between molecular subtypes at the level of breast cancer molecular subtypes. The balanced accuracy of the classification model built on the specified genes was 0.82, while the ROC AUC for the sensitivity analysis was 0.85. Specific patterns of PUFA metabolic changes were obtained for each molecular subtype of breast cancer. These results show evidence that (1) PUFA metabolism genes are critical for the pathogenesis of breast cancer; (2) BC subtypes differ in PUFA metabolism genes expression; and (3) the lists of genes selected in the models are enriched with genes involved in the metabolism of signaling lipids.

Newly Identified Members of FGFR1 Splice Variants Engage in Cross-talk with AXL/AKT Axis in Salivary Adenoid Cystic Carcinoma.

Adenoid cystic carcinoma (ACC) is the second most common malignancy of the salivary gland. Although characterized as an indolent tumor, ACC often leads to incurable metastatic disease. Patients with ACC respond poorly to currently available therapeutic drugs and factors contributing to the limited response remain unknown. Determining the role of molecular alterations frequently occurring in ACC may clarify ACC tumorigenesis and advance the development of effective treatment strategies. Applying Splice Expression Variant Analysis and outlier statistics on RNA sequencing of primary ACC tumors and matched normal salivary gland tissues, we identified multiple alternative splicing events (ASE) of genes specific to ACC. In ACC cells and patient-derived xenografts, FGFR1 was a uniquely expressed ASE. Detailed PCR analysis identified three novel, truncated, intracellular domain-lacking FGFR1 variants (FGFR1v). Cloning and expression analysis suggest that the three FGFR1v are cell surface proteins, that expression of FGFR1v augmented pAKT activity, and that cells became more resistant to pharmacologic FGFR1 inhibitor. FGFR1v-induced AKT activation was associated with AXL function, and inhibition of AXL activity in FGFR1v knockdown cells led to enhanced cytotoxicity in ACC. Moreover, cell killing effect was increased by dual inhibition of AXL and FGFR1 in ACC cells. This study demonstrates that these previously undescribed FGFR1v cooperate with AXL and desensitize cells to FGFR1 inhibitor, which supports further investigation into combined FGFR1 and AXL inhibition as an effective ACC therapy.This study identifies several FGFR1 variants that function through the AXL/AKT signaling pathway independent of FGF/FGFR1, desensitizing cells to FGFR1 inhibitor suggestive of a potential resistance mechanism in ACC. SIGNIFICANCE: This study identifies several FGFR1 variants that function through the AXL/AKT signaling pathway independent of FGF/FGFR1, desensitizing cells to FGFR1 inhibitor, suggestive of a potential resistance mechanism in ACC.

Extracellular Vesicles Released by Tumor Endothelial Cells Spread Immunosuppressive and Transforming Signals Through Various Recipient Cells.

Head and neck squamous cell carcinoma (HNSCC) has a high recurrence and metastatic rate with an unknown mechanism of cancer spread. Tumor inflammation is the most critical processes of cancer onset, growth, and metastasis. We hypothesize that the release of extracellular vesicles (EVs) by tumor endothelial cells (TECs) induce reprogramming of immune cells as well as stromal cells to create an immunosuppressive microenvironment that favor tumor spread. We call this mechanism as non-metastatic contagious carcinogenesis. Extracellular vesicles were collected from primary HNSCC-derived endothelial cells (TEC-EV) and were used for stimulation of peripheral blood mononuclear cells (PBMCs) and primary adipose mesenchymal stem cells (ASCs). Regulation of ASC gene expression was investigated by RNA sequencing and protein array. PBMC, stimulated with TEC-EV, were analyzed by enzyme-linked immunosorbent assay and fluorescence-activated cell sorting. We validated in vitro the effects of TEC-EV on ASCs or PBMC by measuring invasion, adhesion, and proliferation. We found and confirmed that TEC-EV were able to change ASC inflammatory gene expression signature within 24-48 h. TEC-EV were also able to enhance the secretion of TGF-ß1 and IL-10 by PBMC and to increase T regulatory cell (Treg) expansion. TEC-EV carry specific proteins and RNAs that are responsible for Treg differentiation and immune suppression. ASCs and PBMC, treated with TEC-EV, enhanced proliferation, adhesion of tumor cells, and their invasion. These data indicate that TEC-EV exhibit a mechanism of non-metastatic contagious carcinogenesis that regulates tumor microenvironment and reprograms immune cells to sustain tumor growth and progression.

HPV E2, E4, E5 drive alternative carcinogenic pathways in HPV positive cancers.

The dominant paradigm for HPV carcinogenesis includes integration into the host genome followed by expression of E6 and E7 (E6/E7). We explored an alternative carcinogenic pathway characterized by episomal E2, E4, and E5 (E2/E4/E5) expression. Half of HPV positive cervical and pharyngeal cancers comprised a subtype with increase in expression of E2/E4/E5, as well as association with lack of integration into the host genome. Models of the E2/E4/E5 carcinogenesis show p53 dependent enhanced proliferation in vitro, as well as increased susceptibility to induction of cancer in vivo. Whole genomic expression analysis of the E2/E4/E5 pharyngeal cancer subtype is defined by activation of the fibroblast growth factor receptor (FGFR) pathway and this subtype is susceptible to combination FGFR and mTOR inhibition, with implications for targeted therapy.

Rational genomic optimization of DNA detection for human papillomavirus type 16 in head and neck squamous cell carcinoma.

BACKGROUND: We aimed to use genomic data for optimizing polymerase chain reaction (PCR) primer/probe sets for detection of human papillomavirus (HPV)-16 in body fluids of patients with HPV-related head and neck squamous cell carcinoma (HPV-HNSCC). METHODS: We used genomic HPV-HNSCC sequencing data from a single institutional and a TCGA cohort. Optimized primer/probe sets were designed and tested for analytical performance in CaSki HPV-16 genome and confirmed in salivary rinse samples from patients with HPV-HNSCC. RESULTS: The highest read density was observed between E5 and L2 regions. The E1 region contained a region that was universally present. Among candidate PCR primer/probe sets created, six reliably detected 30 HPV-16 copy number. In a CLIA certified laboratory setting, the combination of two novel primer/probe with E7 sets improved performance in salivary rinse samples with a sensitivity of 96% and specificity of 100%. CONCLUSIONS: PCR-based detection of HPV-16 DNA in HPV-HNSCC can be improved using rational genomic design.

Differential Variation Analysis Enables Detection of Tumor Heterogeneity Using Single-Cell RNA-Sequencing Data.

Tumor heterogeneity provides a complex challenge to cancer treatment and is a critical component of therapeutic response, disease recurrence, and patient survival. Single-cell RNA-sequencing (scRNA-seq) technologies have revealed the prevalence of intratumor and intertumor heterogeneity. Computational techniques are essential to quantify the differences in variation of these profiles between distinct cell types, tumor subtypes, and patients to fully characterize intratumor and intertumor molecular heterogeneity. In this study, we adapted our algorithm for pathway dysregulation, Expression Variation Analysis (EVA), to perform multivariate statistical analyses of differential variation of expression in gene sets for scRNA-seq. EVA has high sensitivity and specificity to detect pathways with true differential heterogeneity in simulated data. EVA was applied to several public domain scRNA-seq tumor datasets to quantify the landscape of tumor heterogeneity in several key applications in cancer genomics such as immunogenicity, metastasis, and cancer subtypes. Immune pathway heterogeneity of hematopoietic cell populations in breast tumors corresponded to the amount of diversity present in the T-cell repertoire of each individual. Cells from head and neck squamous cell carcinoma (HNSCC) primary tumors had significantly more heterogeneity across pathways than cells from metastases, consistent with a model of clonal outgrowth. Moreover, there were dramatic differences in pathway dysregulation across HNSCC basal primary tumors. Within the basal primary tumors, there was increased immune dysregulation in individuals with a high proportion of fibroblasts present in the tumor microenvironment. These results demonstrate the broad utility of EVA to quantify intertumor and intratumor heterogeneity from scRNA-seq data without reliance on low-dimensional visualization. SIGNIFICANCE: This study presents a robust statistical algorithm for evaluating gene expression heterogeneity within pathways or gene sets in single-cell RNA-seq data.

Chromatin dysregulation and DNA methylation at transcription start sites associated with transcriptional repression in cancers.

Although promoter-associated CpG islands have been established as targets of DNA methylation changes in cancer, previous studies suggest that epigenetic dysregulation outside the promoter region may be more closely associated with transcriptional changes. Here we examine DNA methylation, chromatin marks, and transcriptional alterations to define the relationship between transcriptional modulation and spatial changes in chromatin structure. Using human papillomavirus-related oropharyngeal carcinoma as a model, we show aberrant enrichment of repressive H3K9me3 at the transcriptional start site (TSS) with methylation-associated, tumor-specific gene silencing. Further analysis identifies a hypermethylated subtype which shows a functional convergence on MYC targets and association with CREBBP/EP300 mutation. The tumor-specific shift to transcriptional repression associated with DNA methylation at TSSs was confirmed in multiple tumor types. Our data may show a common underlying epigenetic dysregulation in cancer associated with broad enrichment of repressive chromatin marks and aberrant DNA hypermethylation at TSSs in combination with MYC network activation.

Publisher Correction: Chromatin dysregulation and DNA methylation at transcription start sites associated with transcriptional repression in cancers.

The original version of this Article contained an error in the author affiliations. Trey Ideker was incorrectly associated with 'Department of Medicine (Oncology), Stanford University School of Medicine, 875 Blake Wilbur Dr, Palo Alto, CA 94304, USA.' This has now been corrected in both the PDF and HTML versions of the Article.

Functional characterization of alternatively spliced GSN in head and neck squamous cell carcinoma.

We have recently performed the characterization of alternative splicing events (ASEs) in head and neck squamous cell carcinoma, which allows dysregulation of protein expression common for cancer cells. Such analysis demonstrated a high ASE prevalence among tumor samples, including tumor-specific alternative splicing in the GSN gene.In vitro studies confirmed that overall expression of either ASE-GSN or wild-type GSN (WT-GSN) isoform inversely correlated with cell proliferation, whereas the high ratio of ASE-GSN to WT-GSN correlated with increased cellular invasion. Additionally, a change in expression of either isoform caused compensatory changes in expression of the other isoform. Our results suggest that the overall expression and the balance between GSN isoforms are mediating factors in proliferation, while increased overall expression of ASE-GSN is specific to cancer tissues. As a result, we propose ASE-GSN can serve not only as a biomarker of disease and disease progression, but also as a neoantigen for head and neck squamous cell carcinoma treatment, for which only a limited number of disease-specific targeted therapies currently exist.

Integrated time course omics analysis distinguishes immediate therapeutic response from acquired resistance.

BACKGROUND: Targeted therapies specifically act by blocking the activity of proteins that are encoded by genes critical for tumorigenesis. However, most cancers acquire resistance and long-term disease remission is rarely observed. Understanding the time course of molecular changes responsible for the development of acquired resistance could enable optimization of patients' treatment options. Clinically, acquired therapeutic resistance can only be studied at a single time point in resistant tumors. METHODS: To determine the dynamics of these molecular changes, we obtained high throughput omics data (RNA-sequencing and DNA methylation) weekly during the development of cetuximab resistance in a head and neck cancer in vitro model. The CoGAPS unsupervised algorithm was used to determine the dynamics of the molecular changes associated with resistance during the time course of resistance development. RESULTS: CoGAPS was used to quantify the evolving transcriptional and epigenetic changes. Applying a PatternMarker statistic to the results from CoGAPS enabled novel heatmap-based visualization of the dynamics in these time course omics data. We demonstrate that transcriptional changes result from immediate therapeutic response or resistance, whereas epigenetic alterations only occur with resistance. Integrated analysis demonstrates delayed onset of changes in DNA methylation relative to transcription, suggesting that resistance is stabilized epigenetically. CONCLUSIONS: Genes with epigenetic alterations associated with resistance that have concordant expression changes are hypothesized to stabilize the resistant phenotype. These genes include FGFR1, which was associated with EGFR inhibitors resistance previously. Thus, integrated omics analysis distinguishes the timing of molecular drivers of resistance. This understanding of the time course progression of molecular changes in acquired resistance is important for the development of alternative treatment strategies that would introduce appropriate selection of new drugs to treat cancer before the resistant phenotype develops.

Splice Expression Variation Analysis (SEVA) for inter-tumor heterogeneity of gene isoform usage in cancer.

Motivation: Current bioinformatics methods to detect changes in gene isoform usage in distinct phenotypes compare the relative expected isoform usage in phenotypes. These statistics model differences in isoform usage in normal tissues, which have stable regulation of gene splicing. Pathological conditions, such as cancer, can have broken regulation of splicing that increases the heterogeneity of the expression of splice variants. Inferring events with such differential heterogeneity in gene isoform usage requires new statistical approaches. Results: We introduce Splice Expression Variability Analysis (SEVA) to model increased heterogeneity of splice variant usage between conditions (e.g. tumor and normal samples). SEVA uses a rank-based multivariate statistic that compares the variability of junction expression profiles within one condition to the variability within another. Simulated data show that SEVA is unique in modeling heterogeneity of gene isoform usage, and benchmark SEVA's performance against EBSeq, DiffSplice and rMATS that model differential isoform usage instead of heterogeneity. We confirm the accuracy of SEVA in identifying known splice variants in head and neck cancer and perform cross-study validation of novel splice variants. A novel comparison of splice variant heterogeneity between subtypes of head and neck cancer demonstrated unanticipated similarity between the heterogeneity of gene isoform usage in HPV-positive and HPV-negative subtypes and anticipated increased heterogeneity among HPV-negative samples with mutations in genes that regulate the splice variant machinery. These results show that SEVA accurately models differential heterogeneity of gene isoform usage from RNA-seq data. Availability and implementation: SEVA is implemented in the R/Bioconductor package GSReg. Contact: bahman@jhu.edu or favorov@sensi.org or ejfertig@jhmi.edu. Supplementary information: Supplementary data are available at Bioinformatics online.

Epigenetic regulation of gene expression in cancer: techniques, resources and analysis.

Cancer is a complex disease, driven by aberrant activity in numerous signaling pathways in even individual malignant cells. Epigenetic changes are critical mediators of these functional changes that drive and maintain the malignant phenotype. Changes in DNA methylation, histone acetylation and methylation, noncoding RNAs, posttranslational modifications are all epigenetic drivers in cancer, independent of changes in the DNA sequence. These epigenetic alterations were once thought to be crucial only for the malignant phenotype maintenance. Now, epigenetic alterations are also recognized as critical for disrupting essential pathways that protect the cells from uncontrolled growth, longer survival and establishment in distant sites from the original tissue. In this review, we focus on DNA methylation and chromatin structure in cancer. The precise functional role of these alterations is an area of active research using emerging high-throughput approaches and bioinformatics analysis tools. Therefore, this review also describes these high-throughput measurement technologies, public domain databases for high-throughput epigenetic data in tumors and model systems and bioinformatics algorithms for their analysis. Advances in bioinformatics data that combine these epigenetic data with genomics data are essential to infer the function of specific epigenetic alterations in cancer. These integrative algorithms are also a focus of this review. Future studies using these emerging technologies will elucidate how alterations in the cancer epigenome cooperate with genetic aberrations during tumor initiation and progression. This deeper understanding is essential to future studies with epigenetics biomarkers and precision medicine using emerging epigenetic therapies.

DNA methylation regulates TMEM16A/ANO1 expression through multiple CpG islands in head and neck squamous cell carcinoma.

ANO1 is a calcium-activated chloride channel that is frequently overexpressed in head and neck squamous cell carcinoma (HNSCC) and other cancers. While ANO1 expression negatively correlates with survival in several cancers, its epigenetic regulation is poorly understood. We analyzed HNSCC samples from TCGA and a separate dataset of HPV+ oropharyngeal squamous cell carcinoma (OPSCC) samples to identify differentially methylated regions. E6 and E7 transfected normal oral keratinocytes (NOK) were used to induce hypermethylation of the ANO1 promoter. We found three CpG islands that correlated with ANO1 expression, including two positively correlated with expression. Using two HNSCC datasets with differential expression of ANO1, we showed hypermethylation of positively correlated CpG islands potentiates ANO1 expression. E7 but not E6 transfection of NOK cells led to hypermethylation of a positively correlated CpG island without a change in ANO1 expression. ANO1 promoter methylation was also correlated with patient survival. Our results are the first to show the contribution of positively correlated CpG's for regulating gene expression in HNSCC. Hypermethylation of the ANO1 promoter was strongly correlated with but not sufficient to increase ANO1 expression, suggesting methylation of positively correlated CpG's likely serves as an adjunct to other mechanisms of ANO1 activation.

Integrated Analysis of Whole-Genome ChIP-Seq and RNA-Seq Data of Primary Head and Neck Tumor Samples Associates HPV Integration Sites with Open Chromatin Marks.

Chromatin alterations mediate mutations and gene expression changes in cancer. Chromatin immunoprecipitation followed by sequencing (ChIP-Seq) has been utilized to study genome-wide chromatin structure in human cancer cell lines, yet numerous technical challenges limit comparable analyses in primary tumors. Here we have developed a new whole-genome analytic pipeline to optimize ChIP-Seq protocols on patient-derived xenografts from human papillomavirus-related (HPV+) head and neck squamous cell carcinoma (HNSCC) samples. We further associated chromatin aberrations with gene expression changes from a larger cohort of the tumor and normal samples with RNA-Seq data. We detect differential histone enrichment associated with tumor-specific gene expression variation, sites of HPV integration in the human genome, and HPV-associated histone enrichment sites upstream of cancer driver genes, which play central roles in cancer-associated pathways. These comprehensive analyses enable unprecedented characterization of the complex network of molecular changes resulting from chromatin alterations that drive HPV-related tumorigenesis. Cancer Res; 77(23); 6538-50. ©2017 AACR.

A Novel Functional Splice Variant of AKT3 Defined by Analysis of Alternative Splice Expression in HPV-Positive Oropharyngeal Cancers.

The incidence of HPV-related oropharyngeal squamous cell carcinoma (OPSCC) has increased more than 200% in the past 20 years. Recent genetic sequencing efforts have elucidated relevant genes in head and neck cancer, but HPV-related tumors have consistently shown few DNA mutations. In this study, we sought to analyze alternative splicing events (ASE) that could alter gene function independent of mutations. To identify ASE unique to HPV-related tumors, RNA sequencing was performed on 46 HPV-positive OPSCC and 25 normal tissue samples. A novel algorithm using outlier statistics on RNA-sequencing junction expression identified 109 splicing events, which were confirmed in a validation set from The Cancer Genome Atlas. Because the most common type of splicing event identified was an alternative start site (39%), MBD-seq genome-wide CpG methylation data were analyzed for methylation alterations at promoter regions. ASE in six genes showed significant negative correlation between promoter methylation and expression of an alternative transcriptional start site, including AKT3 The novel AKT3 transcriptional variant and methylation changes were confirmed using qRT-PCR and qMSP methods. In vitro silencing of the novel AKT3 variant resulted in significant growth inhibition of multiple head and neck cell lines, an effect not observed with wild-type AKT3 knockdown. Analysis of ASE in HPV-related OPSCC identified multiple alterations likely involved in carcinogenesis, including a novel, functionally active transcriptional variant of AKT3 Our data indicate that ASEs represent a significant mechanism of oncogenesis with untapped potential for understanding complex genetic changes that result in the development of cancer. Cancer Res; 77(19); 5248-58. ©2017 AACR.

CoGAPS matrix factorization algorithm identifies transcriptional changes in AP-2alpha target genes in feedback from therapeutic inhibition of the EGFR network.

Patients with oncogene driven tumors are treated with targeted therapeutics including EGFR inhibitors. Genomic data from The Cancer Genome Atlas (TCGA) demonstrates molecular alterations to EGFR, MAPK, and PI3K pathways in previously untreated tumors. Therefore, this study uses bioinformatics algorithms to delineate interactions resulting from EGFR inhibitor use in cancer cells with these genetic alterations. We modify the HaCaT keratinocyte cell line model to simulate cancer cells with constitutive activation of EGFR, HRAS, and PI3K in a controlled genetic background. We then measure gene expression after treating modified HaCaT cells with gefitinib, afatinib, and cetuximab. The CoGAPS algorithm distinguishes a gene expression signature associated with the anticipated silencing of the EGFR network. It also infers a feedback signature with EGFR gene expression itself increasing in cells that are responsive to EGFR inhibitors. This feedback signature has increased expression of several growth factor receptors regulated by the AP-2 family of transcription factors. The gene expression signatures for AP-2alpha are further correlated with sensitivity to cetuximab treatment in HNSCC cell lines and changes in EGFR expression in HNSCC tumors with low CDKN2A gene expression. In addition, the AP-2alpha gene expression signatures are also associated with inhibition of MEK, PI3K, and mTOR pathways in the Library of Integrated Network-Based Cellular Signatures (LINCS) data. These results suggest that AP-2 transcription factors are activated as feedback from EGFR network inhibition and may mediate EGFR inhibitor resistance.

Inferring causal molecular networks: empirical assessment through a community-based effort.

It remains unclear whether causal, rather than merely correlational, relationships in molecular networks can be inferred in complex biological settings. Here we describe the HPN-DREAM network inference challenge, which focused on learning causal influences in signaling networks. We used phosphoprotein data from cancer cell lines as well as in silico data from a nonlinear dynamical model. Using the phosphoprotein data, we scored more than 2,000 networks submitted by challenge participants. The networks spanned 32 biological contexts and were scored in terms of causal validity with respect to unseen interventional data. A number of approaches were effective, and incorporating known biology was generally advantageous. Additional sub-challenges considered time-course prediction and visualization. Our results suggest that learning causal relationships may be feasible in complex settings such as disease states. Furthermore, our scoring approach provides a practical way to empirically assess inferred molecular networks in a causal sense.

Preserving biological heterogeneity with a permuted surrogate variable analysis for genomics batch correction.

MOTIVATION: Sample source, procurement process and other technical variations introduce batch effects into genomics data. Algorithms to remove these artifacts enhance differences between known biological covariates, but also carry potential concern of removing intragroup biological heterogeneity and thus any personalized genomic signatures. As a result, accurate identification of novel subtypes from batch-corrected genomics data is challenging using standard algorithms designed to remove batch effects for class comparison analyses. Nor can batch effects be corrected reliably in future applications of genomics-based clinical tests, in which the biological groups are by definition unknown a priori. RESULTS: Therefore, we assess the extent to which various batch correction algorithms remove true biological heterogeneity. We also introduce an algorithm, permuted-SVA (pSVA), using a new statistical model that is blind to biological covariates to correct for technical artifacts while retaining biological heterogeneity in genomic data. This algorithm facilitated accurate subtype identification in head and neck cancer from gene expression data in both formalin-fixed and frozen samples. When applied to predict Human Papillomavirus (HPV) status, pSVA improved cross-study validation even if the sample batches were highly confounded with HPV status in the training set. AVAILABILITY AND IMPLEMENTATION: All analyses were performed using R version 2.15.0. The code and data used to generate the results of this manuscript is available from https://sourceforge.net/projects/psva.