A Matched Molecular and Clinical Analysis of the Epithelioid Haemangioendothelioma Cohort in the Stafford Fox Rare Cancer Program and Contextual Literature Review.

BACKGROUND: Epithelioid haemangioendothelioma (EHE) is an ultra-rare malignant vascular tumour with a prevalence of 1 per 1,000,000. It is typically molecularly characterised by a WWTR1::CAMTA1 gene fusion in approximately 90% of cases, or a YAP1::TFE3 gene fusion in approximately 10% of cases. EHE cases are typically refractory to therapies, and no anticancer agents are reimbursed for EHE in Australia. METHODS: We report a cohort of nine EHE cases with comprehensive histologic and molecular profiling from the Walter and Eliza Hall Institute of Medical Research Stafford Fox Rare Cancer Program (WEHI-SFRCP) collated via nation-wide referral to the Australian Rare Cancer (ARC) Portal. The diagnoses of EHE were confirmed by histopathological and immunohistochemical (IHC) examination. Molecular profiling was performed using the TruSight Oncology 500 assay, the TruSight RNA fusion panel, whole genome sequencing (WGS), or whole exome sequencing (WES). RESULTS: Molecular analysis of RNA, DNA or both was possible in seven of nine cases. The WWTR1::CAMTA1 fusion was identified in five cases. The YAP1::TFE3 fusion was identified in one case, demonstrating unique morphology compared to cases with the more common WWTR1::CAMTA1 fusion. All tumours expressed typical endothelial markers CD31, ERG, and CD34 and were negative for pan-cytokeratin. Cases with a WWTR1::CAMTA1 fusion displayed high expression of CAMTA1 and the single case with a YAP1::TFE3 fusion displayed high expression of TFE3. Survival was highly variable and unrelated to molecular profile. CONCLUSIONS: This cohort of EHE cases provides molecular and histopathological characterisation and matching clinical information that emphasises the molecular patterns and variable clinical outcomes and adds to our knowledge of this ultra-rare cancer. Such information from multiple studies will advance our understanding, potentially improving treatment options.

Targeting homologous recombination deficiency in uterine leiomyosarcoma.

BACKGROUND: Uterine leiomyosarcoma (uLMS) is a rare and aggressive gynaecological malignancy, with individuals with advanced uLMS having a five-year survival of < 10%. Mutations in the homologous recombination (HR) DNA repair pathway have been observed in ~ 10% of uLMS cases, with reports of some individuals benefiting from poly (ADP-ribose) polymerase (PARP) inhibitor (PARPi) therapy, which targets this DNA repair defect. In this report, we screened individuals with uLMS, accrued nationally, for mutations in the HR repair pathway and explored new approaches to therapeutic targeting. METHODS: A cohort of 58 individuals with uLMS were screened for HR Deficiency (HRD) using whole genome sequencing (WGS), whole exome sequencing (WES) or NGS panel testing. Individuals identified to have HRD uLMS were offered PARPi therapy and clinical outcome details collected. Patient-derived xenografts (PDX) were generated for therapeutic targeting. RESULTS: All 13 uLMS samples analysed by WGS had a dominant COSMIC mutational signature 3; 11 of these had high genome-wide loss of heterozygosity (LOH) (> 0.2) but only two samples had a CHORD score > 50%, one of which had a homozygous pathogenic alteration in an HR gene (deletion in BRCA2). A further three samples harboured homozygous HRD alterations (all deletions in BRCA2), detected by WES or panel sequencing, with 5/58 (9%) individuals having HRD uLMS. All five individuals gained access to PARPi therapy. Two of three individuals with mature clinical follow up achieved a complete response or durable partial response (PR) with the subsequent addition of platinum to PARPi upon minor progression during initial PR on PARPi. Corresponding PDX responses were most rapid, complete and sustained with the PARP1-specific PARPi, AZD5305, compared with either olaparib alone or olaparib plus cisplatin, even in a paired sample of a BRCA2-deleted PDX, derived following PARPi therapy in the patient, which had developed PARPi-resistance mutations in PRKDC, encoding DNA-PKcs. CONCLUSIONS: Our work demonstrates the value of identifying HRD for therapeutic targeting by PARPi and platinum in individuals with the aggressive rare malignancy, uLMS and suggests that individuals with HRD uLMS should be included in trials of PARP1-specific PARPi.

BRCA1 secondary splice-site mutations drive exon-skipping and PARP inhibitor resistance.

BRCA1 splice isoforms Δ11 and Δ11q can contribute to PARP inhibitor (PARPi) resistance by splicing-out the mutation-containing exon, producing truncated, partially-functional proteins. However, the clinical impact and underlying drivers of BRCA1 exon skipping remain undetermined. We analyzed nine ovarian and breast cancer patient derived xenografts (PDX) with BRCA1 exon 11 frameshift mutations for exon skipping and therapy response, including a matched PDX pair derived from a patient pre- and post-chemotherapy/PARPi. BRCA1 exon 11 skipping was elevated in PARPi resistant PDX tumors. Two independent PDX models acquired secondary BRCA1 splice site mutations (SSMs), predicted in silico to drive exon skipping. Predictions were confirmed using qRT-PCR, RNA sequencing, western blots and BRCA1 minigene modelling. SSMs were also enriched in post-PARPi ovarian cancer patient cohorts from the ARIEL2 and ARIEL4 clinical trials. We demonstrate that SSMs drive BRCA1 exon 11 skipping and PARPi resistance, and should be clinically monitored, along with frame-restoring secondary mutations.

CIS controls the functional polarization of GM-CSF-derived macrophages.

The cytokine granulocyte-macrophage-colony stimulating factor (GM-CSF) possesses the capacity to differentiate monocytes into macrophages (MØs) with opposing functions, namely, proinflammatory M1-like MØs and immunosuppressive M2-like MØs. Despite the importance of these opposing biological outcomes, the intrinsic mechanism that regulates the functional polarization of MØs under GM-CSF signaling remains elusive. Here, we showed that GM-CSF-induced MØ polarization resulted in the expression of cytokine-inducible SH2-containing protein (CIS) and that CIS deficiency skewed the differentiation of monocytes toward immunosuppressive M2-like MØs. CIS deficiency resulted in hyperactivation of the JAK-STAT5 signaling pathway, consequently promoting downregulation of the transcription factor Interferon Regulatory Factor 8 (IRF8). Loss- and gain-of-function approaches highlighted IRF8 as a critical regulator of the M1-like polarization program. In vivo, CIS deficiency induced the differentiation of M2-like macrophages, which promoted strong Th2 immune responses characterized by the development of severe experimental asthma. Collectively, our results reveal a CIS-modulated mechanism that clarifies the opposing actions of GM-CSF in MØ differentiation and uncovers the role of GM-CSF in controlling allergic inflammation.

Epithelial-to-Mesenchymal Transition Supports Ovarian Carcinosarcoma Tumorigenesis and Confers Sensitivity to Microtubule Targeting with Eribulin.

Ovarian carcinosarcoma (OCS) is an aggressive and rare tumor type with limited treatment options. OCS is hypothesized to develop via the combination theory, with a single progenitor resulting in carcinomatous and sarcomatous components, or alternatively via the conversion theory, with the sarcomatous component developing from the carcinomatous component through epithelial-to-mesenchymal transition (EMT). In this study, we analyzed DNA variants from isolated carcinoma and sarcoma components to show that OCS from 18 women is monoclonal. RNA sequencing indicated that the carcinoma components were more mesenchymal when compared with pure epithelial ovarian carcinomas, supporting the conversion theory and suggesting that EMT is important in the formation of these tumors. Preclinical OCS models were used to test the efficacy of microtubule-targeting drugs, including eribulin, which has previously been shown to reverse EMT characteristics in breast cancers and induce differentiation in sarcomas. Vinorelbine and eribulin more effectively inhibited OCS growth than standard-of-care platinum-based chemotherapy, and treatment with eribulin reduced mesenchymal characteristics and N-MYC expression in OCS patient-derived xenografts. Eribulin treatment resulted in an accumulation of intracellular cholesterol in OCS cells, which triggered a downregulation of the mevalonate pathway and prevented further cholesterol biosynthesis. Finally, eribulin increased expression of genes related to immune activation and increased the intratumoral accumulation of CD8+ T cells, supporting exploration of immunotherapy combinations in the clinic. Together, these data indicate that EMT plays a key role in OCS tumorigenesis and support the conversion theory for OCS histogenesis. Targeting EMT using eribulin could help improve OCS patient outcomes. SIGNIFICANCE: Genomic analyses and preclinical models of ovarian carcinosarcoma support the conversion theory for disease development and indicate that microtubule inhibitors could be used to suppress EMT and stimulate antitumor immunity.

From a large-scale genomic analysis of insertion sequences to insights into their regulatory roles in prokaryotes.

BACKGROUND: Insertion sequences (ISs) are mobile repeat sequences and most of them can copy themselves to new host genome locations, leading to genome plasticity and gene regulation in prokaryotes. In this study, we present functional and evolutionary relationships between IS and neighboring genes in a large-scale comparative genomic analysis. RESULTS: IS families were located in all prokaryotic phyla, with preferential occurrence of IS3, IS4, IS481, and IS5 families in Alpha-, Beta-, and Gammaproteobacteria, Actinobacteria and Firmicutes as well as in eukaryote host-associated organisms and autotrophic opportunistic pathogens. We defined the concept of the IS-Gene couple (IG), which allowed to highlight the functional and regulatory impacts of an IS on the closest gene. Genes involved in transcriptional regulation and transport activities were found overrepresented in IG. In particular, major facilitator superfamily (MFS) transporters, ATP-binding proteins and transposases raised as favorite neighboring gene functions of IS hotspots. Then, evolutionary conserved IS-Gene sets across taxonomic lineages enabled the classification of IS-gene couples into phylum, class-to-genus, and species syntenic IS-Gene couples. The IS5, IS21, IS4, IS607, IS91, ISL3 and IS200 families displayed two to four times more ISs in the phylum and/or class-to-genus syntenic IGs compared to other IS families. This indicates that those families were probably inserted earlier than others and then subjected to horizontal transfer, transposition and deletion events over time. In phylum syntenic IG category, Betaproteobacteria, Crenarchaeota, Calditrichae, Planctomycetes, Acidithiobacillia and Cyanobacteria phyla act as IS reservoirs for other phyla, and neighboring gene functions are mostly related to transcriptional regulators. Comparison of IS occurrences with predicted regulatory motifs led to ~ 26.5% of motif-containing ISs with 2 motifs per IS in average. These results, concomitantly with short IS-Gene distances, suggest that those ISs would interfere with the expression of neighboring genes and thus form strong candidates for an adaptive pairing. CONCLUSIONS: All together, our large-scale study provide new insights into the IS genetic context and strongly suggest their regulatory roles.

Lung Imaging Reveals Stroke-Induced Impairment in Pulmonary Intravascular Neutrophil Function, a Response Exacerbated with Aging.

In stroke patients, infection is a significant contributor to morbidity and mortality. Moreover, older stroke patients show an increased risk of developing stroke-associated infection, although the mechanisms underlying this increased susceptibility to infection are unknown. In this study, using an experimental mouse model of ischemic stroke, we showed that older (12-15 mo of age) mice had elevated lung bacterial infection and inflammatory damage after stroke when compared with young (8-10 wk of age) counterparts, despite undergoing the same degree of brain injury. Intravital microscopy of the lung microvasculature revealed that in younger mice, stroke promoted neutrophil arrest in pulmonary microvessels, but this response was not seen in older poststroke mice. In addition, bacterial phagocytosis by neutrophils in the lung microvasculature was reduced by both aging and stroke, such that neutrophils in aged poststroke mice showed the greatest impairment in this function. Analysis of neutrophil migration in vitro and in the cremaster muscle demonstrated that stroke alone did not negatively impact neutrophil migration, but that the combination of increased age and stroke led to reduced effectiveness of neutrophil chemotaxis. Transcriptomic analysis of pulmonary neutrophils using RNA sequencing identified 79 genes that were selectively altered in the context of combined aging and stroke, and they were associated with pathways that control neutrophil chemotaxis. Taken together, the findings of this study show that stroke in older animals results in worsening of neutrophil antibacterial responses and changes in neutrophil gene expression that have the potential to underpin elevated risk of stroke-associated infection in the context of increased age.

Blood-Derived Extracellular Vesicle-Associated miR-3182 Detects Non-Small Cell Lung Cancer Patients.

With five-year survival rates as low as 3%, lung cancer is the most common cause of cancer-related mortality worldwide. The severity of the disease at presentation is accredited to the lack of early detection capacities, resulting in the reliance on low-throughput diagnostic measures, such as tissue biopsy and imaging. Interest in the development and use of liquid biopsies has risen, due to non-invasive sample collection, and the depth of information it can provide on a disease. Small extracellular vesicles (sEVs) as viable liquid biopsies are of particular interest due to their potential as cancer biomarkers. To validate the use of sEVs as cancer biomarkers, we characterised cancer sEVs using miRNA sequencing analysis. We found that miRNA-3182 was highly enriched in sEVs derived from the blood of patients with invasive breast carcinoma and NSCLC. The enrichment of sEV miR-3182 was confirmed in oncogenic, transformed lung cells in comparison to isogenic, untransformed lung cells. Most importantly, miR-3182 can successfully distinguish early-stage NSCLC patients from those with benign lung conditions. Therefore, miR-3182 provides potential to be used for the detection of NSCLC in blood samples, which could result in earlier therapy and thus improved outcomes and survival for patients.

Integrating deep mutational scanning and low-throughput mutagenesis data to predict the impact of amino acid variants.

BACKGROUND: Evaluating the impact of amino acid variants has been a critical challenge for studying protein function and interpreting genomic data. High-throughput experimental methods like deep mutational scanning (DMS) can measure the effect of large numbers of variants in a target protein, but because DMS studies have not been performed on all proteins, researchers also model DMS data computationally to estimate variant impacts by predictors. RESULTS: In this study, we extended a linear regression-based predictor to explore whether incorporating data from alanine scanning (AS), a widely used low-throughput mutagenesis method, would improve prediction results. To evaluate our model, we collected 146 AS datasets, mapping to 54 DMS datasets across 22 distinct proteins. CONCLUSIONS: We show that improved model performance depends on the compatibility of the DMS and AS assays, and the scale of improvement is closely related to the correlation between DMS and AS results.

svaRetro and svaNUMT: modular packages for annotating retrotransposed transcripts and nuclear integration of mitochondrial DNA in genome sequencing data.

Nuclear integration of mitochondrial genomes and retrocopied transcript insertion are biologically important but often-overlooked aspects of structural variant (SV) annotation. While tools for their detection exist, these typically rely on reanalysis of primary data using specialised detectors rather than leveraging calls from general purpose structural variant callers. Such reanalysis potentially leads to additional computational expense and does not take advantage of advances in general purpose structural variant calling. Here, we present svaRetro and svaNUMT; R packages that provide functions for annotating novel genomic events, such as nonreference retrocopied transcripts and nuclear integration of mitochondrial DNA. The packages were developed to work within the Bioconductor framework. We evaluate the performance of these packages to detect events using simulations and public benchmarking datasets, and annotate processed transcripts in a public structural variant database. svaRetro and svaNUMT provide modular, SV-caller agnostic tools for downstream annotation of structural variant calls.

Ductal variant prostate carcinoma is associated with a significantly shorter metastasis-free survival.

BACKGROUND: Ductal adenocarcinoma is an uncommon prostate cancer variant. Previous studies suggest that ductal variant histology may be associated with worse clinical outcomes, but these are difficult to interpret. To address this, we performed an international, multi-institutional study to describe the characteristics of ductal adenocarcinoma, particularly focussing on the effect of presence of ductal variant cancer on metastasis-free survival. METHODS: Patients with ductal variant histology from two institutional databases who underwent radical prostatectomies were identified and compared with an independent acinar adenocarcinoma cohort. After propensity score matching, the effect of the presence of ductal adenocarcinoma on time to biochemical recurrence, initiation of salvage therapy and the development of metastatic disease was determined. Deep whole-exome sequencing was performed for selected cases (n = 8). RESULTS: A total of 202 ductal adenocarcinoma and 2037 acinar adenocarcinoma cases were analysed. Survival analysis after matching demonstrated that patients with ductal variant histology had shorter salvage-free survival (8.1 versus 22.0 months, p = 0.03) and metastasis-free survival (6.7 versus 78.6 months, p < 0.0001). Ductal variant histology was consistently associated with RB1 loss, as well as copy number gains in TAP1, SLC4A2 and EHHADH. CONCLUSIONS: The presence of any ductal variant adenocarcinoma at the time of prostatectomy portends a worse clinical outcome than pure acinar cancers, with significantly shorter times to initiation of salvage therapies and the onset of metastatic disease. These features appear to be driven by uncoupling of chromosomal duplication from cell division, resulting in widespread copy number aberration with specific gain of genes implicated in treatment resistance.

Unifying package managers, workflow engines, and containers: Computational reproducibility with BioNix.

MOTIVATION: A challenge for computational biologists is to make our analyses reproducible-i.e. to rerun, combine, and share, with the assurance that equivalent runs will generate identical results. Current best practice aims at this using a combination of package managers, workflow engines, and containers. RESULTS: We present BioNix, a lightweight library built on the Nix deployment system. BioNix manages software dependencies, computational environments, and workflow stages together using a single abstraction: pure functions. This lets users specify workflows in a clean, uniform way, with strong reproducibility guarantees. AVAILABILITY AND IMPLEMENTATION: BioNix is implemented in the Nix expression language and is released on GitHub under the 3-clause BSD license: https://github.com/PapenfussLab/bionix (biotools:BioNix) (BioNix, RRID:SCR_017662).

Improved HLA-based prediction of coeliac disease identifies two novel genetic interactions.

Human Leucocyte Antigen (HLA) testing is useful in the clinical work-up of coeliac disease (CD) with high negative but low positive predictive value. We construct a genomic risk score (GRS) using HLA risk genotypes to improve CD prediction and guide exclusion criteria. Imputed HLA genotypes for five European CD case-control GWAS (n > 15,000) were used to construct and validate an interpretable HLA-based risk model (HDQ15), which shows statistically significant improvements in predictive performance upon all previous HLA-based risk models. Conditioning on this model, we find two novel associations, HLA-DQ6.2 and HLA-DQ7.3, that interact significantly with HLA-DQ2.5 (p = 2.51 × 10-9, 1.99 × 10-7, respectively). Integrating these novel alleles into a new risk model (HDQ17) leads to predictive performance equivalent or better than the strongest reported GRS (GRS228) using 228 single nucleotide polymorphisms (SNPs). We also demonstrate that our proposed HLA-based models can be implemented using only six HLA tagging SNPs with statistically equivalent predictive performance. Using insights from our model to guide exclusionary criteria, we find the positive predictive value of CD testing in high-risk populations can be increased by 55%, from 17.5 to 27.1%, while maintaining a negative predictive value above 99%. Our results suggest that HLA typing is currently undervalued in CD assessment.

MTOR signaling orchestrates stress-induced mutagenesis, facilitating adaptive evolution in cancer.

In microorganisms, evolutionarily conserved mechanisms facilitate adaptation to harsh conditions through stress-induced mutagenesis (SIM). Analogous processes may underpin progression and therapeutic failure in human cancer. We describe SIM in multiple in vitro and in vivo models of human cancers under nongenotoxic drug selection, paradoxically enhancing adaptation at a competing intrinsic fitness cost. A genome-wide approach identified the mechanistic target of rapamycin (MTOR) as a stress-sensing rheostat mediating SIM across multiple cancer types and conditions. These observations are consistent with a two-phase model for drug resistance, in which an initially rapid expansion of genetic diversity is counterbalanced by an intrinsic fitness penalty, subsequently normalizing to complete adaptation under the new conditions. This model suggests synthetic lethal strategies to minimize resistance to anticancer therapy.

BioShake: a Haskell EDSL for bioinformatics workflows.

Typical bioinformatics analyses comprise of long running computational workflows. An important part of reproducible research is the management and execution of these workflows to allow robust execution and to minimise errors. BioShake is an embedded domain specific language in Haskell for specifying and executing computational workflows for bioinformatics that significantly reduces the possibility of errors occurring. Unlike other workflow frameworks, BioShake raises many properties to the type level allowing the correctness of a workflow to be statically checked during compilation, catching errors before any lengthy execution process. BioShake builds on the Shake build tool to provide robust dependency tracking, parallel execution, reporting, and resumption capabilities. Finally, BioShake abstracts execution so that jobs can either be executed directly or submitted to a cluster. BioShake is available at http://github.com/PapenfussLab/bioshake.

Spleen Volume Variation in Patients with Locally Advanced Non-Small Cell Lung Cancer Receiving Platinum-Based Chemo-Radiotherapy.

There is renewed interest in the immune regulatory role of the spleen in oncology. To date, very few studies have examined macroscopic variations of splenic volume in the setting of cancer, prior to or during therapy, especially in humans. Changes in splenic volume may be associated with changes in splenic function. The purpose of this study was to investigate variations in spleen volume in NSCLC patients during chemo-radiotherapy. Sixty patients with stage I-IIIB NSCLC underwent radiotherapy (60 Gy/30 fractions) for six weeks with concomitant carboplatin/paclitaxel (Ca/P; n = 32) or cisplatin/etoposide (Ci/E; n = 28). A baseline PET/CT scan was performed within 2 weeks prior to treatment and during Weeks 2 and 4 of chemo-radiotherapy. Spleen volume was measured by contouring all CT slices. Significant macroscopic changes in splenic volume occurred early after the commencement of treatment. A significant decrease in spleen volume was observed for 66% of Ca/P and 79% of Ci/E patients between baseline and Week 2. Spleen volume was decreased by 14.2% for Ca/P (p<0.001) and 19.3% for Ci/E (p<0.001) patients. By Week 4, spleen volume was still significantly decreased for Ca/P patients compared to baseline, while for Ci/E patients, spleen volume returned to above baseline levels. This is the first report demonstrating macroscopic changes in the spleen in NSCLC patients undergoing radical chemo-radiotherapy that can be visualized by non-invasive imaging.

RAR/RXR binding dynamics distinguish pluripotency from differentiation associated cis-regulatory elements.

In mouse embryonic cells, ligand-activated retinoic acid receptors (RARs) play a key role in inhibiting pluripotency-maintaining genes and activating some major actors of cell differentiation. To investigate the mechanism underlying this dual regulation, we performed joint RAR/RXR ChIP-seq and mRNA-seq time series during the first 48 h of the RA-induced Primitive Endoderm (PrE) differentiation process in F9 embryonal carcinoma (EC) cells. We show here that this dual regulation is associated with RAR/RXR genomic redistribution during the differentiation process. In-depth analysis of RAR/RXR binding sites occupancy dynamics and composition show that in undifferentiated cells, RAR/RXR interact with genomic regions characterized by binding of pluripotency-associated factors and high prevalence of the non-canonical DR0-containing RA response element. By contrast, in differentiated cells, RAR/RXR bound regions are enriched in functional Sox17 binding sites and are characterized with a higher frequency of the canonical DR5 motif. Our data offer an unprecedentedly detailed view on the action of RA in triggering pluripotent cell differentiation and demonstrate that RAR/RXR action is mediated via two different sets of regulatory regions tightly associated with cell differentiation status.

Development and validation of a gene expression tumour classifier for cancer of unknown primary.

Accurate identification of the primary tumour in cancer of unknown primary (CUP) is required for effective treatment selection and improved patient outcomes. The aim of this study was to develop and validate a gene expression tumour classifier and integrate it with histopathology to identify the likely site of origin in CUP.RNA was extracted from 450 formalin fixed, paraffin embedded samples of known origin comprising 18 tumour groups. Whole genome expression analysis was performed using a bead-based array. Classification of the tumours made use of a binary support vector machine, together with recursive feature elimination. A hierarchical tumour classifier was developed and incorporated with conventional histopathology to identify the origins of metastatic tumours.The classifier demonstrated an accuracy of 88% for correctly predicting the tumour type on a validation set of known tumours (n = 94). For CUP samples (n = 49) having a final clinical diagnosis, the classifier improved the accuracy of histology alone for both single and multiple predictions. Furthermore, where histology alone could not suggest any specific diagnosis, the classifier was able to correctly predict the primary site of origin.We demonstrate the integration of gene expression profiling with conventional histopathology to aid the investigation of CUP.

Stability of bivariate GWAS biomarker detection.

Given the difficulty and effort required to confirm candidate causal SNPs detected in genome-wide association studies (GWAS), there is no practical way to definitively filter false positives. Recent advances in algorithmics and statistics have enabled repeated exhaustive search for bivariate features in a practical amount of time using standard computational resources, allowing us to use cross-validation to evaluate the stability. We performed 10 trials of 2-fold cross-validation of exhaustive bivariate analysis on seven Wellcome-Trust Case-Control Consortium GWAS datasets, comparing the traditional [Formula: see text] test for association, the high-performance GBOOST method and the recently proposed GSS statistic (Available at http://bioinformatics.research.nicta.com.au/software/gwis/). We use Spearman's correlation to measure the similarity between the folds of cross validation. To compare incomplete lists of ranks we propose an extension to Spearman's correlation. The extension allows us to consider a natural threshold for feature selection where the correlation is zero. This is the first reported cross-validation study of exhaustive bivariate GWAS feature selection. We found that stability between ranked lists from different cross-validation folds was higher for GSS in the majority of diseases. A thorough analysis of the correlation between SNP-frequency and univariate [Formula: see text] score demonstrated that the [Formula: see text] test for association is highly confounded by main effects: SNPs with high univariate significance replicably dominate the ranked results. We show that removal of the univariately significant SNPs improves [Formula: see text] replicability but risks filtering pairs involving SNPs with univariate effects. We empirically confirm that the stability of GSS and GBOOST were not affected by removal of univariately significant SNPs. These results suggest that the GSS and GBOOST tests are successfully targeting bivariate association with phenotype and that GSS is able to reliably detect a larger set of SNP-pairs than GBOOST in the majority of the data we analysed. However, the [Formula: see text] test for association was confounded by main effects.