Divergent molecular networks program functionally distinct CD8+ skin-resident memory T cells

Skin-resident CD8+ T cells include distinct interferon-γ-producing [tissue-resident memory T type 1 (TRM1)] and interleukin-17 (IL-17)-producing (TRM17) subsets that differentially contribute to immune responses. However, whether these populations use common mechanisms to establish tissue residence is unknown. In this work, we show that TRM1 and TRM17 cells navigate divergent trajectories to acquire tissue residency in the skin. TRM1 cells depend on a T-bet-Hobit-IL-15 axis, whereas TRM17 cells develop independently of these factors. Instead, c-Maf commands a tissue-resident program in TRM17 cells parallel to that induced by Hobit in TRM1 cells, with an ICOS-c-Maf-IL-7 axis pivotal to TRM17 cell commitment. Accordingly, by targeting this pathway, skin TRM17 cells can be ablated without compromising their TRM1 counterparts. Thus, skin-resident T cells rely on distinct molecular circuitries, which can be exploited to strategically modulate local immunity.

Establishing the Link between X-Chromosome Aberrations and TP53 Status, with Breast Cancer Patient Outcomes

Ubiquitous to normal female human somatic cells, X-chromosome inactivation (XCI) tightly regulates the transcriptional silencing of a single X chromosome from each pair. Some genes escape XCI, including crucial tumour suppressors. Cancer susceptibility can be influenced by the variability in the genes that escape XCI. The mechanisms of XCI dysregulation remain poorly understood in complex diseases, including cancer. Using publicly available breast cancer next-generation sequencing data, we show that the status of the major tumour suppressor TP53 from Chromosome 17 is highly associated with the genomic integrity of the inactive X (Xi) and the active X (Xa) chromosomes. Our quantification of XCI and XCI escape demonstrates that aberrant XCI is linked to poor survival. We derived prognostic gene expression signatures associated with either large deletions of Xi; large amplifications of Xa; or abnormal X-methylation. Our findings expose a novel insight into female cancer risks, beyond those associated with the standard molecular subtypes.

sccomp: Robust differential composition and variability analysis for single-cell data

Cellular omics such as single-cell genomics, proteomics, and microbiomics allow the characterization of tissue and microbial community composition, which can be compared between conditions to identify biological drivers. This strategy has been critical to revealing markers of disease progression, such as cancer and pathogen infection. A dedicated statistical method for differential variability analysis is lacking for cellular omics data, and existing methods for differential composition analysis do not model some compositional data properties, suggesting there is room to improve model performance. Here, we introduce sccomp, a method for differential composition and variability analyses that jointly models data count distribution, compositionality, group-specific variability, and proportion mean-variability association, being aware of outliers. sccomp provides a comprehensive analysis framework that offers realistic data simulation and cross-study knowledge transfer. Here, we demonstrate that mean-variability association is ubiquitous across technologies, highlighting the inadequacy of the very popular Dirichlet-multinomial distribution. We show that sccomp accurately fits experimental data, significantly improving performance over state-of-the-art algorithms. Using sccomp, we identified differential constraints and composition in the microenvironment of primary breast cancer.

An open protocol for modeling T Cell Clonotype repertoires using TCRβ CDR3 sequences

T cell receptor repertoires can be profiled using next generation sequencing (NGS) to measure and monitor adaptive dynamical changes in response to disease and other perturbations. Genomic DNA-based bulk sequencing is cost-effective but necessitates multiplex target amplification using multiple primer pairs with highly variable amplification efficiencies. Here, we utilize an equimolar primer mixture and propose a single statistical normalization step that efficiently corrects for amplification bias post sequencing. Using samples analyzed by both our open protocol and a commercial solution, we show high concordance between bulk clonality metrics. This approach is an inexpensive and open-source alternative to commercial solutions.

Spatial analysis with SPIAT and spaSim to characterize and simulate tissue microenvironments

Spatial proteomics technologies have revealed an underappreciated link between the location of cells in tissue microenvironments and the underlying biology and clinical features, but there is significant lag in the development of downstream analysis methods and benchmarking tools. Here we present SPIAT (spatial image analysis of tissues), a spatial-platform agnostic toolkit with a suite of spatial analysis algorithms, and spaSim (spatial simulator), a simulator of tissue spatial data. SPIAT includes multiple colocalization, neighborhood and spatial heterogeneity metrics to characterize the spatial patterns of cells. Ten spatial metrics of SPIAT are benchmarked using simulated data generated with spaSim. We show how SPIAT can uncover cancer immune subtypes correlated with prognosis in cancer and characterize cell dysfunction in diabetes. Our results suggest SPIAT and spaSim as useful tools for quantifying spatial patterns, identifying and validating correlates of clinical outcomes and supporting method development.

Early immune pressure initiated by tissue-resident memory T cells sculpts tumor evolution in non-small cell lung cancer

Tissue-resident memory T (T_RM) cells provide immune defense against local infection and can inhibit cancer progression. However, it is unclear to what extent chronic inflammation impacts T_RM activation and whether T_RM cells existing in tissues before tumor onset influence cancer evolution in humans. We performed deep profiling of healthy lungs and lung cancers in never-smokers (NSs) and ever-smokers (ESs), finding evidence of enhanced immunosurveillance by cells with a T_RM-like phenotype in ES lungs. In preclinical models, tumor-specific or bystander T_RM-like cells present prior to tumor onset boosted immune cell recruitment, causing tumor immune evasion through loss of MHC class I protein expression and resistance to immune checkpoint inhibitors. In humans, only tumors arising in ES patients underwent clonal immune evasion, unrelated to tobacco-associated mutagenic signatures or oncogenic drivers. These data demonstrate that enhanced T_RM-like activity prior to tumor development shapes the evolution of tumor immunogenicity and can impact immunotherapy outcomes.

Genomic characterisation of hormone receptor-positive breast cancer arising in very young women

BACKGROUND

Very young premenopausal women diagnosed with hormone receptor-positive, human epidermal growth factor receptor 2-negative (HR+HER2-) early breast cancer (EBC) have higher rates of recurrence and death for reasons that remain largely unexplained.

PATIENTS AND METHODS

Genomic sequencing was applied to HR+HER2- tumours from patients enrolled in the Suppression of Ovarian Function Trial (SOFT) to determine genomic drivers that are enriched in young premenopausal women. Genomic alterations were characterised using next-generation sequencing from a subset of 1276 patients (deep targeted sequencing, n = 1258; whole-exome sequencing in a young-age, case-control subsample, n = 82). We defined copy number (CN) subgroups and assessed for features suggestive of homologous recombination deficiency (HRD). Genomic alteration frequencies were compared between young premenopausal women (<40 years) and older premenopausal women (≥40 years), and assessed for associations with distant recurrence-free interval (DRFI) and overall survival (OS).

RESULTS

Younger women (<40 years, n = 359) compared with older women (≥40 years, n = 917) had significantly higher frequencies of mutations in GATA3 (19% versus 16%) and CN amplifications (CNAs) (47% versus 26%), but significantly lower frequencies of mutations in PIK3CA (32% versus 47%), CDH1 (3% versus 9%), and MAP3K1 (7% versus 12%). Additionally, they had significantly higher frequencies of features suggestive of HRD (27% versus 21%) and a higher proportion of PIK3CA mutations with concurrent CNAs (23% versus 11%). Genomic features suggestive of HRD, PIK3CA mutations with CNAs, and CNAs were associated with significantly worse DRFI and OS compared with those without these features. These poor prognostic features were enriched in younger patients: present in 72% of patients aged <35 years, 54% aged 35-39 years, and 40% aged ≥40 years. Poor prognostic features [n = 584 (46%)] versus none [n = 692 (54%)] had an 8-year DRFI of 84% versus 94% and OS of 88% versus 96%. Younger women (<40 years) had the poorest outcomes: 8-year DRFI 74% versus 85% and OS 80% versus 93%, respectively.

CONCLUSION

These results provide insights into genomic alterations that are enriched in young women with HR+HER2- EBC, provide rationale for genomic subgrouping, and highlight priority molecular targets for future clinical trials.

Intratumoral CD8+ T cells with a tissue-resident memory phenotype mediate local immunity and immune checkpoint responses in breast cancer

CD8⁺ tumor-infiltrating lymphocytes with a tissue-resident memory T (T_RM) cell phenotype are associated with favorable prognosis in patients with triple-negative breast cancer (TNBC). However, the relative contribution of CD8⁺ T_RM cells to anti-tumor immunity and immune checkpoint blockade efficacy in breast cancer remains unknown. Here, we show that intratumoral CD8⁺ T cells in murine mammary tumors transcriptionally resemble those from TNBC patients. Phenotypic and transcriptional studies established two intratumoral sub-populations: one more enriched in markers of terminal exhaustion (T_EX-like) and the other with a bona fide resident phenotype (T_RM-like). Treatment with anti-PD-1 and anti-CTLA-4 therapy resulted in expansion of these intratumoral populations, with the T_RM-like subset displaying significantly enhanced cytotoxic capacity. T_RM-like CD8⁺ T cells could also provide local immune protection against tumor rechallenge and a T_RM gene signature extracted from tumor-free tissue was significantly associated with improved clinical outcomes in TNBC patients treated with checkpoint inhibitors.

An open protocol for modeling T Cell Clonotype repertoires using TCRβ CDR3 sequences

T cell receptor repertoires can be profiled using next generation sequencing (NGS) to measure and monitor adaptive dynamical changes in response to disease and other perturbations. Genomic DNA-based bulk sequencing is cost-effective but necessitates multiplex target amplification using multiple primer pairs with highly variable amplification efficiencies. Here, we utilize an equimolar primer mixture and propose a single statistical normalization step that efficiently corrects for amplification bias post sequencing. Using samples analyzed by both our open protocol and a commercial solution, we show high concordance between bulk clonality metrics. This approach is an inexpensive and open-source alternative to commercial solutions.

Removing unwanted variation from large-scale RNA sequencing data with PRPS

Accurate identification and effective removal of unwanted variation is essential to derive meaningful biological results from RNA sequencing (RNA-seq) data, especially when the data come from large and complex studies. Using RNA-seq data from The Cancer Genome Atlas (TCGA), we examined several sources of unwanted variation and demonstrate here how these can significantly compromise various downstream analyses, including cancer subtype identification, association between gene expression and survival outcomes and gene co-expression analysis. We propose a strategy, called pseudo-replicates of pseudo-samples (PRPS), for deploying our recently developed normalization method, called removing unwanted variation III (RUV-III), to remove the variation caused by library size, tumor purity and batch effects in TCGA RNA-seq data. We illustrate the value of our approach by comparing it to the standard TCGA normalizations on several TCGA RNA-seq datasets. RUV-III with PRPS can be used to integrate and normalize other large transcriptomic datasets coming from multiple laboratories or platforms.

Runx3 drives a CD8+ T cell tissue residency program that is absent in CD4+ T cells

Tissue-resident memory T cells (TRM cells) provide rapid and superior control of localized infections. While the transcription factor Runx3 is a critical regulator of CD8+ T cell tissue residency, its expression is repressed in CD4+ T cells. Here, we show that, as a direct consequence of this Runx3-deficiency, CD4+ TRM cells lacked the transforming growth factor (TGF)-β-responsive transcriptional network that underpins the tissue residency of epithelial CD8+ TRM cells. While CD4+ TRM cell formation required Runx1, this, along with the modest expression of Runx3 in CD4+ TRM cells, was insufficient to engage the TGF-β-driven residency program. Ectopic expression of Runx3 in CD4+ T cells incited this TGF-β-transcriptional network to promote prolonged survival, decreased tissue egress, a microanatomical redistribution towards epithelial layers and enhanced effector functionality. Thus, our results reveal distinct programming of tissue residency in CD8+ and CD4+ TRM cell subsets that is attributable to divergent Runx3 activity.

RUV-III-NB: normalization of single cell RNA-seq data

Normalization of single cell RNA-seq data remains a challenging task. The performance of different methods can vary greatly between datasets when unwanted factors and biology are associated. Most normalization methods also only remove the effects of unwanted variation for the cell embedding but not from gene-level data typically used for differential expression (DE) analysis to identify marker genes. We propose RUV-III-NB, a method that can be used to remove unwanted variation from both the cell embedding and gene-level counts. Using pseudo-replicates, RUV-III-NB explicitly takes into account potential association with biology when removing unwanted variation. The method can be used for both UMI or read counts and returns adjusted counts that can be used for downstream analyses such as clustering, DE and pseudotime analyses. Using published datasets with different technological platforms, kinds of biology and levels of association between biology and unwanted variation, we show that RUV-III-NB manages to remove library size and batch effects, strengthen biological signals, improve DE analyses, and lead to results exhibiting greater concordance with independent datasets of the same kind. The performance of RUV-III-NB is consistent and is not sensitive to the number of factors assumed to contribute to the unwanted variation.

A mouse-specific retrotransposon drives a conserved Cdk2ap1 isoform essential for development

Retrotransposons mediate gene regulation in important developmental and pathological processes. Here, we characterized the transient retrotransposon induction during preimplantation development of eight mammals. Induced retrotransposons exhibit similar preimplantation profiles across species, conferring gene regulatory activities, particularly through long terminal repeat (LTR) retrotransposon promoters. A mouse-specific MT2B2 retrotransposon promoter generates an N-terminally truncated Cdk2ap1ΔN that peaks in preimplantation embryos and promotes proliferation. In contrast, the canonical Cdk2ap1 peaks in mid-gestation and represses cell proliferation. This MT2B2 promoter, whose deletion abolishes Cdk2ap1ΔN production, reduces cell proliferation and impairs embryo implantation, is developmentally essential. Intriguingly, Cdk2ap1ΔN is evolutionarily conserved in sequence and function yet is driven by different promoters across mammals. The distinct preimplantation Cdk2ap1ΔN expression in each mammalian species correlates with the duration of its preimplantation development. Hence, species-specific transposon promoters can yield evolutionarily conserved, alternative protein isoforms, bestowing them with new functions and species-specific expression to govern essential biological divergence.

The Ratio of Exhausted to Resident Infiltrating Lymphocytes Is Prognostic for Colorectal Cancer Patient Outcome

Immunotherapy success in colorectal cancer is mainly limited to patients whose tumors exhibit high microsatellite instability (MSI). However, there is variability in treatment outcomes within this group, which is in part driven by the frequency and characteristics of tumor-infiltrating immune cells. Indeed, the presence of specific infiltrating immune-cell subsets has been shown to correlate with immunotherapy response and is in many cases prognostic of treatment outcome. Tumor-infiltrating lymphocytes (TIL) can undergo distinct differentiation programs, acquiring features of tissue-residency or exhaustion, a process during which T cells upregulate inhibitory receptors, such as PD-1, and lose functionality. Although residency and exhaustion programs of CD8⁺ T cells are relatively well studied, these programs have only recently been appreciated in CD4⁺ T cells and remain largely unknown in tumor-infiltrating natural killer (NK) cells. In this study, we used single-cell RNA sequencing (RNA-seq) data to identify signatures of residency and exhaustion in colorectal cancer-infiltrating lymphocytes, including CD8⁺, CD4⁺, and NK cells. We then tested these signatures in independent single-cell data from tumor and normal tissue-infiltrating immune cells. Furthermore, we used versions of these signatures designed for bulk RNA-seq data to explore tumor-intrinsic mutations associated with residency and exhaustion from TCGA data. Finally, using two independent transcriptomic datasets from patients with colon adenocarcinoma, we showed that combinations of these signatures, in particular combinations of NK-cell activity signatures, together with tumor-associated signatures, such as TGFβ signaling, were associated with distinct survival outcomes in patients with colon adenocarcinoma.

Finding a suitable library size to call variants in RNA-Seq

BACKGROUND

RNA sequencing allows the study of both gene expression changes and transcribed mutations, providing a highly effective way to gain insight into cancer biology. When planning the sequencing of a large cohort of samples, library size is a fundamental factor affecting both the overall cost and the quality of the results. Here we specifically address how overall library size influences the detection of somatic mutations in RNA-seq data in two acute myeloid leukaemia datasets. RESULTS : We simulated shallower sequencing depths by downsampling 45 acute myeloid leukaemia samples (100 bp PE) that are part of the Leucegene project, which were originally sequenced at high depth. We compared the sensitivity of six methods of recovering validated mutations on the same samples. The methods compared are a combination of three popular callers (MuTect, VarScan, and VarDict) and two filtering strategies. We observed an incremental loss in sensitivity when simulating libraries of 80M, 50M, 40M, 30M and 20M fragments, with the largest loss detected with less than 30M fragments (below 90%, average loss of 7%). The sensitivity in recovering insertions and deletions varied markedly between callers, with VarDict showing the highest sensitivity (60%). Single nucleotide variant sensitivity is relatively consistent across methods, apart from MuTect, whose default filters need adjustment when using RNA-Seq. We also analysed 136 RNA-Seq samples from the TCGA-LAML cohort (50 bp PE) and assessed the change in sensitivity between the initial libraries (average 59M fragments) and after downsampling to 40M fragments. When considering single nucleotide variants in recurrently mutated myeloid genes we found a comparable performance, with a 6% average loss in sensitivity using 40M fragments.

CONCLUSIONS

Between 30M and 40M 100 bp PE reads are needed to recover 90-95% of the initial variants on recurrently mutated myeloid genes. To extend this result to another cancer type, an exploration of the characteristics of its mutations and gene expression patterns is suggested.

Correction to: Serum microRNA is a biomarker for post-operative monitoring in glioma

For the reference citation '[57]' in the second paragraph of the Results section of the original article there was no corresponding entry in the References section. It should have referred to the below mentioned article by Ebrahimkhani et al. (2018).

Removing unwanted variation with CytofRUV to integrate multiple CyTOF datasets

Mass cytometry (CyTOF) is a technology that has revolutionised single-cell biology. By detecting over 40 proteins on millions of single cells, CyTOF allows the characterisation of cell subpopulations in unprecedented detail. However, most CyTOF studies require the integration of data from multiple CyTOF batches usually acquired on different days and possibly at different sites. To date, the integration of CyTOF datasets remains a challenge due to technical differences arising in multiple batches. To overcome this limitation, we developed an approach called CytofRUV for analysing multiple CyTOF batches, which includes an R-Shiny application with diagnostic plots. CytofRUV can correct for batch effects and integrate data from large numbers of patients and conditions across batches, to confidently compare cellular changes and correlate these with clinically relevant outcomes.

Strategies to enable large-scale proteomics for reproducible research

Reproducible research is the bedrock of experimental science. To enable the deployment of large-scale proteomics, we assess the reproducibility of mass spectrometry (MS) over time and across instruments and develop computational methods for improving quantitative accuracy. We perform 1560 data independent acquisition (DIA)-MS runs of eight samples containing known proportions of ovarian and prostate cancer tissue and yeast, or control HEK293T cells. Replicates are run on six mass spectrometers operating continuously with varying maintenance schedules over four months, interspersed with ~5000 other runs. We utilise negative controls and replicates to remove unwanted variation and enhance biological signal, outperforming existing methods. We also design a method for reducing missing values. Integrating these computational modules into a pipeline (ProNorM), we mitigate variation among instruments over time and accurately predict tissue proportions. We demonstrate how to improve the quantitative analysis of large-scale DIA-MS data, providing a pathway toward clinical proteomics.

Clinical proteomics critically depends on the ability to acquire highly reproducible data over an extended period of time. Here, the authors assess reproducibility over four months across different mass spectrometers and develop a computational approach to mitigate variation among instruments over time.

Changes in ferrous iron and glutathione promote ferroptosis and frailty in aging Caenorhabditis elegans

All eukaryotes require iron. Replication, detoxification, and a cancer-protective form of regulated cell death termed ferroptosis, all depend on iron metabolism. Ferrous iron accumulates over adult lifetime in Caenorhabditis elegans. Here, we show that glutathione depletion is coupled to ferrous iron elevation in these animals, and that both occur in late life to prime cells for ferroptosis. We demonstrate that blocking ferroptosis, either by inhibition of lipid peroxidation or by limiting iron retention, mitigates age-related cell death and markedly increases lifespan and healthspan. Temporal scaling of lifespan is not evident when ferroptosis is inhibited, consistent with this cell death process acting at specific life phases to induce organismal frailty, rather than contributing to a constant aging rate. Because excess age-related iron elevation in somatic tissue, particularly in brain, is thought to contribute to degenerative disease, post-developmental interventions to limit ferroptosis may promote healthy aging.

DECENT: differential expression with capture efficiency adjustmeNT for single-cell RNA-seq data

Motivation

Dropout is a common phenomenon in single-cell RNA-seq (scRNA-seq) data, and when left unaddressed it affects the validity of the statistical analyses. Despite this, few current methods for differential expression (DE) analysis of scRNA-seq data explicitly model the process that gives rise to the dropout events. We develop DECENT, a method for DE analysis of scRNA-seq data that explicitly and accurately models the molecule capture process in scRNA-seq experiments.

Results

We show that DECENT demonstrates improved DE performance over existing DE methods that do not explicitly model dropout. This improvement is consistently observed across several public scRNA-seq datasets generated using different technological platforms. The gain in improvement is especially large when the capture process is overdispersed. DECENT maintains type I error well while achieving better sensitivity. Its performance without spike-ins is almost as good as when spike-ins are used to calibrate the capture model.

Availability and implementation

The method is implemented as a publicly available R package available from https://github.com/cz-ye/DECENT.

Supplementary information

Supplementary data are available at Bioinformatics online.

Multiple sclerosis risk variants regulate gene expression in innate and adaptive immune cells

At least 200 single-nucleotide polymorphisms (SNPs) are associated with multiple sclerosis (MS) risk. A key function that could mediate SNP-encoded MS risk is their regulatory effects on gene expression. We performed microarrays using RNA extracted from purified immune cell types from 73 untreated MS cases and 97 healthy controls and then performed Cis expression quantitative trait loci mapping studies using additive linear models. We describe MS risk expression quantitative trait loci associations for 129 distinct genes. By extending these models to include an interaction term between genotype and phenotype, we identify MS risk SNPs with opposing effects on gene expression in cases compared with controls, namely, rs2256814 MYT1 in CD4 cells (q = 0.05) and rs12087340 RF00136 in monocyte cells (q = 0.04). The rs703842 SNP was also associated with a differential effect size on the expression of the METTL21B gene in CD8 cells of MS cases relative to controls (q = 0.03). Our study provides a detailed map of MS risk loci that function by regulating gene expression in cell types relevant to MS.

scClassify: sample size estimation and multiscale classification of cells using single and multiple reference

Automated cell type identification is a key computational challenge in single-cell RNA-sequencing (scRNA-seq) data. To capitalise on the large collection of well-annotated scRNA-seq datasets, we developed scClassify, a multiscale classification framework based on ensemble learning and cell type hierarchies constructed from single or multiple annotated datasets as references. scClassify enables the estimation of sample size required for accurate classification of cell types in a cell type hierarchy and allows joint classification of cells when multiple references are available. We show that scClassify consistently performs better than other supervised cell type classification methods across 114 pairs of reference and testing data, representing a diverse combination of sizes, technologies and levels of complexity, and further demonstrate the unique components of scClassify through simulations and compendia of experimental datasets. Finally, we demonstrate the scalability of scClassify on large single-cell atlases and highlight a novel application of identifying subpopulations of cells from the Tabula Muris data that were unidentified in the original publication. Together, scClassify represents state-of-the-art methodology in automated cell type identification from scRNA-seq data.

A new normalization for Nanostring nCounter gene expression data

The Nanostring nCounter gene expression assay uses molecular barcodes and single molecule imaging to detect and count hundreds of unique transcripts in a single reaction. These counts need to be normalized to adjust for the amount of sample, variations in assay efficiency and other factors. Most users adopt the normalization approach described in the nSolver analysis software, which involves background correction based on the observed values of negative control probes, a within-sample normalization using the observed values of positive control probes and normalization across samples using reference (housekeeping) genes. Here we present a new normalization method, Removing Unwanted Variation-III (RUV-III), which makes vital use of technical replicates and suitable control genes. We also propose an approach using pseudo-replicates when technical replicates are not available. The effectiveness of RUV-III is illustrated on four different datasets. We also offer suggestions on the design and analysis of studies involving this technology.

Identification of cancer sex-disparity in the functional integrity of p53 and its X chromosome network

The disproportionately high prevalence of male cancer is poorly understood. We tested for sex-disparity in the functional integrity of the major tumor suppressor p53 in sporadic cancers. Our bioinformatics analyses expose three novel levels of p53 impact on sex-disparity in 12 non-reproductive cancer types. First, TP53 mutation is more frequent in these cancers among US males than females, with poorest survival correlating with its mutation. Second, numerous X-linked genes are associated with p53, including vital genomic regulators. Males are at unique risk from alterations of their single copies of these genes. High expression of X-linked negative regulators of p53 in wild-type TP53 cancers corresponds with reduced survival. Third, females exhibit an exceptional incidence of non-expressed mutations among p53-associated X-linked genes. Our data indicate that poor survival in males is contributed by high frequencies of TP53 mutations and an inability to shield against deregulated X-linked genes that engage in p53 networks.

Using long-read sequencing to detect imprinted DNA methylation

Systematic variation in the methylation of cytosines at CpG sites plays a critical role in early development of humans and other mammals. Of particular interest are regions of differential methylation between parental alleles, as these often dictate monoallelic gene expression, resulting in parent of origin specific control of the embryonic transcriptome and subsequent development, in a phenomenon known as genomic imprinting. Using long-read nanopore sequencing we show that, with an average genomic coverage of ∼10, it is possible to determine both the level of methylation of CpG sites and the haplotype from which each read arises. The long-read property is exploited to characterize, using novel methods, both methylation and haplotype for reads that have reduced basecalling precision compared to Sanger sequencing. We validate the analysis both through comparison of nanopore-derived methylation patterns with those from Reduced Representation Bisulfite Sequencing data and through comparison with previously reported data. Our analysis successfully identifies known imprinting control regions (ICRs) as well as some novel differentially methylated regions which, due to their proximity to hitherto unknown monoallelically expressed genes, may represent new ICRs.

Evaluation of cross-platform and interlaboratory concordance via consensus modelling of genomic measurements

Motivation

A synoptic view of the human genome benefits chiefly from the application of nucleic acid sequencing and microarray technologies. These platforms allow interrogation of patterns such as gene expression and DNA methylation at the vast majority of canonical loci, allowing granular insights and opportunities for validation of original findings. However, problems arise when validating against a “gold standard” measurement, since this immediately biases all subsequent measurements towards that particular technology or protocol. Since all genomic measurements are estimates, in the absence of a ”gold standard” we instead empirically assess the measurement precision and sensitivity of a large suite of genomic technologies via a consensus modelling method called the row-linear model. This method is an application of the American Society for Testing and Materials Standard E691 for assessing interlaboratory precision and sources of variability across multiple testing sites. Both cross-platform and cross-locus comparisons can be made across all common loci, allowing identification of technology- and locus-specific tendencies.

Results

We assess technologies including the Infinium MethylationEPIC BeadChip, whole genome bisulfite sequencing (WGBS), two different RNA-Seq protocols (PolyA+ and Ribo-Zero) and five different gene expression array platforms. Each technology thus is characterised herein, relative to the consensus. We showcase a number of applications of the row-linear model, including correlation with known interfering traits. We demonstrate a clear effect of cross-hybridisation on the sensitivity of Infinium methylation arrays. Additionally, we perform a true interlaboratory test on a set of samples interrogated on the same platform across twenty-one separate testing laboratories.

Availability and implementation

A full implementation of the row-linear model, plus extra functions for visualisation, are found in the R package consensus at https://github.com/timpeters82/consensus.

Supplementary information

Supplementary data are available at Bioinformatics online.

Dysregulated IL-1β-GM-CSF Axis in Acute Rheumatic Fever That Is Limited by Hydroxychloroquine

BACKGROUND

Acute rheumatic fever (ARF) and rheumatic heart disease are autoimmune consequences of group A streptococcus infection and remain major causes of cardiovascular morbidity and mortality around the world. Improved treatment has been stymied by gaps in understanding key steps in the immunopathogenesis of ARF and rheumatic heart disease. This study aimed to identify (1) effector T cell cytokine(s) that might be dysregulated in the autoimmune response of patients with ARF by group A streptococcus, and (2) an immunomodulatory agent that suppresses this response and could be clinically translatable to high-risk patients with ARF.

METHODS

The immune response to group A streptococcus was analyzed in peripheral blood mononuclear cells from an Australian Aboriginal ARF cohort by a combination of multiplex cytokine array, flow cytometric analysis, and global gene expression analysis by RNA sequencing. The immunomodulatory drug hydroxychloroquine was tested for effects on this response.

RESULTS

We found a dysregulated interleukin-1β-granulocyte-macrophage colony-stimulating factor (GM-CSF) cytokine axis in ARF peripheral blood mononuclear cells exposed to group A streptococcus in vitro, whereby persistent interleukin-1β production is coupled to overproduction of GM-CSF and selective expansion of CXCR3⁺CCR4^-CCR6^- CD4 T cells. CXCR3⁺CCR4^-CCR6^- CD4 T cells are the major source of GM-CSF in human CD4 T cells and CXCL10, a CXCR3 ligand and potent T helper 1 chemoattractant, was elevated in sera from patients with ARF. GM-CSF has recently emerged as a key T cell-derived effector cytokine in numerous autoimmune diseases, including myocarditis, and the production of CXCL10 may explain selective trafficking of these cells to the heart. We provide evidence that interleukin-1β amplifies the expansion of GM-CSF-expressing CD4 T cells, which is effectively suppressed by hydroxychloroquine. RNA sequencing showed shifts in gene expression profiles and differentially expressed genes in peripheral blood mononuclear cells derived from patients at different clinical stages of ARF.

CONCLUSIONS

Given the safety profile of hydroxychloroquine and its clinical pedigree in treating autoimmune diseases such as rheumatoid arthritis, where GM-CSF plays a pivotal role, we propose that hydroxychloroquine could be repurposed to reduce the risk of rheumatic heart disease after ARF.

Abstract PD5-03: Characterization of high TIL breast cancers reveals a prognostic and functionally distinct tissue-resident memory subpopulation

Background: Tumor infiltrating lymphocytes (TILs) assessed via light microscopy are prognostic and predictive in the early stage and advanced triple negative and HER2-amplified breast cancer (BC). Higher TILs can also identify patients more likely to benefit from anti-PD-1 therapy. In this study we interrogated T cell subsets that comprise high TILs to determine if distinct subpopulations are key mediators of anti-tumor immunity.

Methods: We characterised TILs with a focus on CD3+ T cells in 129 primary and metastatic BC samples using flow cytometry, bulk RNASeq on flow sorted T cell populations, multiplex immunohistochemistry and microdroplet based single cell 3' mRNA sequencing on the 10X Genomics Chromium platform. Cell type specific gene expression signatures were determined from differential expression between putative T cell subpopulations. These signatures were investigated in clinical cohorts, including trial cohorts treated with pembrolizumab.

Results: High TIL Infiltrates consisted primarily of CD3+ T cells, with both CD8 and CD4 populations. Unsupervised clustering of single cell sequencing identified 9 CD8 and CD4 subpopulations with distinct gene expression profiles. In addition to Tregs and CD8 effector memory (TEM) T cells, we found a CD8+ tissue resident memory (TRM) population expressing greater levels of T-cell checkpoints and cytotoxic markers compared to effector memory cells. In 2 primary tumours and 1 liver metastasis, bulk RNASeq of flow sorted TEM and TRM corroborated the single cell mRNASeq results. T cell receptor profiling (TCR) in the 3 samples found non-overlapping repertoires in the 2 primary tumours, but overlap in one metastatic lesion, suggesting divergent developmental origins in the breast, but the potential for nascent TRM differentiation in a metastatic niche. Clustering of these TCRs suggested differing antigen specificities between TRM and non-TRM CD8 T cells. Using Metabric data, the CD8 TRM gene expression signature was prognostic for disease free survival (DFS) in primary TNBCs (n=329, log-rank p=0.003), and was able to further stratify cases with high and low CD8A expression for DFS (log-rank p = 0.03). The CD8 TRM signature was enriched in baseline tumour samples of responders (n = 9) compared with non-responders (n=36) in 45 patients with metastatic melanoma treated with T cell checkpoint blockade (p &lt; 0.0001). Additional single cell sequencing data with TCR sequencing will be combined with these initial results, and an independent data set of single cell mRNASeq and TCR Seq on CD3+ BC TILs will be used to confirm our findings. Cell type specific signatures will be explored in additional clinical cohorts including KEYNOTE-086, and presented at the meeting.

Conclusion: Using single cell profiling of the immune microenvironment in BC we demonstrate that high TIL BCs contain multiple T cell subpopulations with different functional and prognostic significance. Our approach identified a CD8 TRM population with a distinct gene expression profile and strong expression of key immune checkpoints likely representing the presence of true tumor specific immunity. This population may be a key target of immune checkpoint blockade.

Citation Format: Savas P, Virassamy B, Ye C, Salim A, Mintoff CP, Caramia F, Salgado R, Teo ZL, Dushyanthen S, Byrne A, Luen SJ, Fox SB, Speed TP, Mackay LK, Neeson PJ, Loi S. Characterization of high TIL breast cancers reveals a prognostic and functionally distinct tissue-resident memory subpopulation [abstract]. In: Proceedings of the 2018 San Antonio Breast Cancer Symposium; 2018 Dec 4-8; San Antonio, TX. Philadelphia (PA): AACR; Cancer Res 2019;79(4 Suppl):Abstract nr PD5-03.

Maps of variability in cell lineage trees

New approaches to lineage tracking have allowed the study of differentiation in multicellular organisms over many generations of cells. Understanding the phenotypic variability observed in these lineage trees requires new statistical methods. Whereas an invariant cell lineage, such as that for the nematode Caenorhabditis elegans, can be described by a lineage map, defined as the pattern of phenotypes overlaid onto the binary tree, a traditional lineage map is static and does not describe the variability inherent in the cell lineages of higher organisms. Here, we introduce lineage variability maps which describe the pattern of second-order variation in lineage trees. These maps can be undirected graphs of the partial correlations between every lineal position, or directed graphs showing the dynamics of bifurcated patterns in each subtree. We show how to infer these graphical models for lineages of any depth from sample sizes of only a few pedigrees. This required developing the generalized spectral analysis for a binary tree, the natural framework for describing tree-structured variation. When tested on pedigrees from C. elegans expressing a marker for pharyngeal differentiation potential, the variability maps recover essential features of the known lineage map. When applied to highly-variable pedigrees monitoring cell size in T lymphocytes, the maps show that most of the phenotype is set by the founder naive T cell. Lineage variability maps thus elevate the concept of the lineage map to the population level, addressing questions about the potency and dynamics of cell lineages and providing a way to quantify the progressive restriction of cell fate with increasing depth in the tree.

Joint Modeling of Mixed Plasmodium Species Infections Using a Bivariate Poisson Lognormal Model

Infectious diseases often present as coinfections that may affect each other in positive or negative ways. Understanding the relationship between two coinfecting pathogens is thus important to understand the risk of infection and burden of disease caused by each pathogen. Although coinfections with Plasmodium falciparum and Plasmodium vivax are very common outside Africa, it is yet unclear whether infections by the two parasite species are positively associated or if infection by one parasite suppresses the other. In this study, we use bivariate Poisson lognormal models (BPLM) to estimate covariate-adjusted associations between the incidence of infections (as measured by the force of blood-stage infections, _molFOI) and clinical episodes caused by both P. falciparum and P. vivax in a cohort of Papua New Guinean children. A BPLM permits estimation of either positive or negative correlation, unlike most other multivariate Poisson models. Our results demonstrated a moderately positive association between P. falciparum and P. vivax infection rates, arguing against the hypothesis that P. vivax infections protect against P. falciparum infections. Our findings also suggest that the BPLM is only useful for counts with suitably large means and overdispersion.

Publisher Correction: Single-cell profiling of breast cancer T cells reveals a tissue-resident memory subset associated with improved prognosis

In the version of this article originally published, the institution in affiliation 10 was missing. Affiliation 10 was originally listed as Department of Surgery, Royal Melbourne Hospital and Royal Womens' Hospital, Melbourne, Victoria, Australia. It should have been Department of Surgery, Royal Melbourne Hospital and Royal Womens' Hospital, University of Melbourne, Melbourne, Victoria, Australia. The error has been corrected in the HTML and PDF versions of this article.

The healthy ageing gene expression signature for Alzheimer's disease diagnosis: a random sampling perspective

In a recent publication, Sood et al. (Genome Biol 16:185, 2015) presented a set of 150 probe sets that could be used in the diagnosis of Alzheimer’s disease (AD) based on gene expression. We reproduce some of their experiments and show that their signature is indeed able to discriminate between AD and control patients using blood gene expression in two cohorts. We also show that its performance does not stand out compared to randomly sampled sets of 150 probe sets from the same array.

Single-cell profiling of breast cancer T cells reveals a tissue-resident memory subset associated with improved prognosis

The quantity of tumor-infiltrating lymphocytes (TILs) in breast cancer (BC) is a robust prognostic factor for improved patient survival, particularly in triple-negative and HER2-overexpressing BC subtypes¹. Although T cells are the predominant TIL population², the relationship between quantitative and qualitative differences in T cell subpopulations and patient prognosis remains unknown. We performed single-cell RNA sequencing (scRNA-seq) of 6,311 T cells isolated from human BCs and show that significant heterogeneity exists in the infiltrating T cell population. We demonstrate that BCs with a high number of TILs contained CD8⁺ T cells with features of tissue-resident memory T (T_RM) cell differentiation and that these CD8⁺ T_RM cells expressed high levels of immune checkpoint molecules and effector proteins. A CD8⁺ T_RM gene signature developed from the scRNA-seq data was significantly associated with improved patient survival in early-stage triple-negative breast cancer (TNBC) and provided better prognostication than CD8 expression alone. Our data suggest that CD8⁺ T_RM cells contribute to BC immunosurveillance and are the key targets of modulation by immune checkpoint inhibition. Further understanding of the development, maintenance and regulation of T_RM cells will be crucial for successful immunotherapeutic development in BC.

Correction to: A statistical framework for analyzing deep mutational scanning data

CORRECTION

After publication of our article [1] it was brought to our attention that a line of code was missing from our program to combine the within-replicate variance and between-replicate variance. This led to an overestimation of the standard errors calculated using the Enrich2 random-effects model.

Accurate RNA Sequencing From Formalin-Fixed Cancer Tissue To Represent High-Quality Transcriptome From Frozen Tissue

Purpose

Accurate transcriptional sequencing (RNA-seq) from formalin-fixation and paraffin-embedding (FFPE) tumor samples presents an important challenge for translational research and diagnostic development. In addition, there are now several different protocols to prepare a sequencing library from total RNA. We evaluated the accuracy of RNA-seq data generated from FFPE samples in terms of expression profiling.

Methods

We designed a biospecimen study to directly compare gene expression results from different protocols to prepare libraries for RNA-seq from human breast cancer tissues, with randomization to fresh-frozen (FF) or FFPE conditions. The protocols were compared using multiple computational methods to assess alignment of reads to reference genome, and the uniformity and continuity of coverage; as well as the variance and correlation, of overall gene expression and patterns of measuring coding sequence, phenotypic patterns of gene expression, and measurements from representative multigene signatures.

Results

The principal determinant of variance in gene expression was use of exon capture probes, followed by the conditions of preservation (FF versus FFPE), and phenotypic differences between breast cancers. One protocol, with RNase H-based rRNA depletion, exhibited least variability of gene expression measurements, strongest correlation between FF and FFPE samples, and was generally representative of the transcriptome from standard FF RNA-seq protocols.

Conclusion

Method of RNA-seq library preparation from FFPE samples had marked effect on the accuracy of gene expression measurement compared to matched FF samples. Nevertheless, some protocols produced highly concordant expression data from FFPE RNA-seq data, compared to RNA-seq results from matched frozen samples.

A statistical framework for analyzing deep mutational scanning data

Deep mutational scanning is a widely used method for multiplex measurement of functional consequences of protein variants. We developed a new deep mutational scanning statistical model that generates error estimates for each measurement, capturing both sampling error and consistency between replicates. We apply our model to one novel and five published datasets comprising 243,732 variants and demonstrate its superiority in removing noisy variants and conducting hypothesis testing. Simulations show our model applies to scans based on cell growth or binding and handles common experimental errors. We implemented our model in Enrich2, software that can empower researchers analyzing deep mutational scanning data.

Correction: The Subclonal Architecture of Metastatic Breast Cancer: Results from a Prospective Community-Based Rapid Autopsy Program "CASCADE"

[This corrects the article DOI: 10.1371/journal.pmed.1002204.].

Deficiency of microRNA miR-34a expands cell fate potential in pluripotent stem cells

Embryonic stem cells (ESCs) and induced pluripotent stem cells (iPSCs) efficiently generate all embryonic cell lineages but rarely generate extraembryonic cell types. We found that microRNA miR-34a deficiency expands the developmental potential of mouse pluripotent stem cells, yielding both embryonic and extraembryonic lineages and strongly inducing MuERV-L (MERVL) endogenous retroviruses, similar to what is seen with features of totipotent two-cell blastomeres. miR-34a restricts the acquisition of expanded cell fate potential in pluripotent stem cells, and it represses MERVL expression through transcriptional regulation, at least in part by targeting the transcription factor Gata2. Our studies reveal a complex molecular network that defines and restricts pluripotent developmental potential in cultured ESCs and iPSCs.

The Subclonal Architecture of Metastatic Breast Cancer: Results from a Prospective Community-Based Rapid Autopsy Program "CASCADE"

BACKGROUND

Understanding the cancer genome is seen as a key step in improving outcomes for cancer patients. Genomic assays are emerging as a possible avenue to personalised medicine in breast cancer. However, evolution of the cancer genome during the natural history of breast cancer is largely unknown, as is the profile of disease at death. We sought to study in detail these aspects of advanced breast cancers that have resulted in lethal disease.

METHODS AND FINDINGS

Three patients with oestrogen-receptor (ER)-positive, human epidermal growth factor receptor 2 (HER2)-negative breast cancer and one patient with triple negative breast cancer underwent rapid autopsy as part of an institutional prospective community-based rapid autopsy program (CASCADE). Cases represented a range of management problems in breast cancer, including late relapse after early stage disease, de novo metastatic disease, discordant disease response, and disease refractory to treatment. Between 5 and 12 metastatic sites were collected at autopsy together with available primary tumours and longitudinal metastatic biopsies taken during life. Samples underwent paired tumour-normal whole exome sequencing and single nucleotide polymorphism (SNP) arrays. Subclonal architectures were inferred by jointly analysing all samples from each patient. Mutations were validated using high depth amplicon sequencing. Between cases, there were significant differences in mutational burden, driver mutations, mutational processes, and copy number variation. Within each case, we found dramatic heterogeneity in subclonal structure from primary to metastatic disease and between metastatic sites, such that no single lesion captured the breadth of disease. Metastatic cross-seeding was found in each case, and treatment drove subclonal diversification. Subclones displayed parallel evolution of treatment resistance in some cases and apparent augmentation of key oncogenic drivers as an alternative resistance mechanism. We also observed the role of mutational processes in subclonal evolution. Limitations of this study include the potential for bias introduced by joint analysis of formalin-fixed archival specimens with fresh specimens and the difficulties in resolving subclones with whole exome sequencing. Other alterations that could define subclones such as structural variants or epigenetic modifications were not assessed.

CONCLUSIONS

This study highlights various mechanisms that shape the genome of metastatic breast cancer and the value of studying advanced disease in detail. Treatment drives significant genomic heterogeneity in breast cancers which has implications for disease monitoring and treatment selection in the personalised medicine paradigm.

Removal of unwanted variation reveals novel patterns of gene expression linked to sleep homeostasis in murine cortex

BACKGROUND

Why we sleep is still one of the most perplexing mysteries in biology. Strong evidence indicates that sleep is necessary for normal brain function and that sleep need is a tightly regulated process. Surprisingly, molecular mechanisms that determine sleep need are incompletely described. Moreover, very little is known about transcriptional changes that specifically accompany the accumulation and discharge of sleep need. Several studies have characterized differential gene expression changes following sleep deprivation. Much less is known, however, about changes in gene expression during the compensatory response to sleep deprivation (i.e. recovery sleep).

RESULTS

In this study we present a comprehensive analysis of the effects of sleep deprivation and subsequent recovery sleep on gene expression in the mouse cortex. We used a non-traditional analytical method for normalization of genome-wide gene expression data, Removal of Unwanted Variation (RUV). RUV improves detection of differential gene expression following sleep deprivation. We also show that RUV normalization is crucial to the discovery of differentially expressed genes associated with recovery sleep. Our analysis indicates that the majority of transcripts upregulated by sleep deprivation require 6 h of recovery sleep to return to baseline levels, while the majority of downregulated transcripts return to baseline levels within 1-3 h. We also find that transcripts that change rapidly during recovery (i.e. within 3 h) do so on average with a time constant that is similar to the time constant for the discharge of sleep need.

CONCLUSIONS

We demonstrate that proper data normalization is essential to identify changes in gene expression that are specifically linked to sleep deprivation and recovery sleep. Our results provide the first evidence that recovery sleep is comprised of two waves of transcriptional regulation that occur at different times and affect functionally distinct classes of genes.

Contextual fear conditioning induces differential alternative splicing

The process of memory consolidation requires transcription and translation to form long-term memories. Significant effort has been dedicated to understanding changes in hippocampal gene expression after contextual fear conditioning. However, alternative splicing by differential transcript regulation during this time period has received less attention. Here, we use RNA-seq to determine exon-level changes in expression after contextual fear conditioning and retrieval. Our work reveals that a short variant of Homer1, Ania-3, is regulated by contextual fear conditioning. The ribosome biogenesis regulator Las1l, small nucleolar RNA Snord14e, and the RNA-binding protein Rbm3 also change specific transcript usage after fear conditioning. The changes in Ania-3 and Las1l are specific to either the new context or the context-shock association, while the changes in Rbm3 occur after context or shock only. Our analysis revealed novel transcript regulation of previously undetected changes after learning, revealing the importance of high throughput sequencing approaches in the study of gene expression changes after learning.

Signatures of tumour immunity distinguish Asian and non-Asian gastric adenocarcinomas

OBJECTIVE

Differences in gastric cancer (GC) clinical outcomes between patients in Asian and non-Asian countries has been historically attributed to variability in clinical management. However, recent international Phase III trials suggest that even with standardised treatments, GC outcomes differ by geography. Here, we investigated gene expression differences between Asian and non-Asian GCs, and if these molecular differences might influence clinical outcome.

DESIGN

We compared gene expression profiles of 1016 GCs from six Asian and three non-Asian GC cohorts, using a two-stage meta-analysis design and a novel biostatistical method (RUV-4) to adjust for technical variation between cohorts. We further validated our findings by computerised immunohistochemical analysis on two independent tissue microarray (TMA) cohorts from Asian and non-Asian localities (n=665).

RESULTS

Gene signatures differentially expressed between Asians and non-Asian GCs were related to immune function and inflammation. Non-Asian GCs were significantly enriched in signatures related to T-cell biology, including CTLA-4 signalling. Similarly, in the TMA cohorts, non-Asian GCs showed significantly higher expression of T-cell markers (CD3, CD45R0, CD8) and lower expression of the immunosuppressive T-regulatory cell marker FOXP3 compared to Asian GCs (p<0.05). Inflammatory cell markers CD66b and CD68 also exhibited significant cohort differences (p<0.05). Exploratory analyses revealed a significant relationship between tumour immunity factors, geographic locality-specific prognosis, and postchemotherapy outcomes.

CONCLUSIONS

Analyses of >1600 GCs suggest that Asian and non-Asian GCs exhibit distinct tumour immunity signatures related to T-cell function. These differences may influence geographical differences in clinical outcome, and the design of future trials particularly in immuno-oncology.

Assessment of DNA methylation profiling and copy number variation as indications of clonal relationship in ipsilateral and contralateral breast cancers to distinguish recurrent breast cancer from a second primary tumour

BACKGROUND

Patients with breast cancer have an increased risk of developing subsequent breast cancers. It is important to distinguish whether these tumours are de novo or recurrences of the primary tumour in order to guide the appropriate therapy. Our aim was to investigate the use of DNA methylation profiling and array comparative genomic hybridization (aCGH) to determine whether the second tumour is clonally related to the first tumour.

METHODS

Methylation-sensitive high-resolution melting was used to screen promoter methylation in a panel of 13 genes reported as methylated in breast cancer (RASSF1A, TWIST1, APC, WIF1, MGMT, MAL, CDH13, RARβ, BRCA1, CDH1, CDKN2A, TP73, and GSTP1) in 29 tumour pairs (16 ipsilateral and 13 contralateral). Using the methylation profile of these genes, we employed a Bayesian and an empirical statistical approach to estimate clonal relationship. Copy number alterations were analysed using aCGH on the same set of tumour pairs.

RESULTS

There is a higher probability of the second tumour being recurrent in ipsilateral tumours compared with contralateral tumours (38 % versus 8 %; p <0.05) based on the methylation profile. Using previously reported recurrence rates as Bayesian prior probabilities, we classified 69 % of ipsilateral and 15 % of contralateral tumours as recurrent. The inferred clonal relationship results of the tumour pairs were generally concordant between methylation profiling and aCGH.

CONCLUSION

Our results show that DNA methylation profiling as well as aCGH have potential as diagnostic tools in improving the clinical decisions to differentiate recurrences from a second de novo tumour.

Systematic noise degrades gene co-expression signals but can be corrected

Background

In the past decade, the identification of gene co-expression has become a routine part of the analysis of high-dimensional microarray data. Gene co-expression, which is mostly detected via the Pearson correlation coefficient, has played an important role in the discovery of molecular pathways and networks. Unfortunately, the presence of systematic noise in high-dimensional microarray datasets corrupts estimates of gene co-expression. Removing systematic noise from microarray data is therefore crucial. Many cleaning approaches for microarray data exist, however these methods are aimed towards improving differential expression analysis and their performances have been primarily tested for this application. To our knowledge, the performances of these approaches have never been systematically compared in the context of gene co-expression estimation.

Results

Using simulations we demonstrate that standard cleaning procedures, such as background correction and quantile normalization, fail to adequately remove systematic noise that affects gene co-expression and at times further degrade true gene co-expression. Instead we show that a global version of removal of unwanted variation (RUV), a data-driven approach, removes systematic noise but also allows the estimation of the true underlying gene-gene correlations. We compare the performance of all noise removal methods when applied to five large published datasets on gene expression in the human brain. RUV retrieves the highest gene co-expression values for sets of genes known to interact, but also provides the greatest consistency across all five datasets. We apply the method to prioritize epileptic encephalopathy candidate genes.

Conclusions

Our work raises serious concerns about the quality of many published gene co-expression analyses. RUV provides an efficient and flexible way to remove systematic noise from high-dimensional microarray datasets when the objective is gene co-expression analysis. The RUV method as applicable in the context of gene-gene correlation estimation is available as a BioconductoR-package: RUVcorr.

Electronic supplementary material

The online version of this article (doi:10.1186/s12859-015-0745-3) contains supplementary material, which is available to authorized users.

Correcting gene expression data when neither the unwanted variation nor the factor of interest are observed

When dealing with large scale gene expression studies, observations are commonly contaminated by sources of unwanted variation such as platforms or batches. Not taking this unwanted variation into account when analyzing the data can lead to spurious associations and to missing important signals. When the analysis is unsupervised, e.g. when the goal is to cluster the samples or to build a corrected version of the dataset--as opposed to the study of an observed factor of interest--taking unwanted variation into account can become a difficult task. The factors driving unwanted variation may be correlated with the unobserved factor of interest, so that correcting for the former can remove the latter if not done carefully. We show how negative control genes and replicate samples can be used to estimate unwanted variation in gene expression, and discuss how this information can be used to correct the expression data. The proposed methods are then evaluated on synthetic data and three gene expression datasets. They generally manage to remove unwanted variation without losing the signal of interest and compare favorably to state-of-the-art corrections. All proposed methods are implemented in the bioconductor package RUVnormalize.

How data analysis affects power, reproducibility and biological insight of RNA-seq studies in complex datasets

The sequencing of the full transcriptome (RNA-seq) has become the preferred choice for the measurement of genome-wide gene expression. Despite its widespread use, challenges remain in RNA-seq data analysis. One often-overlooked aspect is normalization. Despite the fact that a variety of factors or ‘batch effects’ can contribute unwanted variation to the data, commonly used RNA-seq normalization methods only correct for sequencing depth. The study of gene expression is particularly problematic when it is influenced simultaneously by a variety of biological factors in addition to the one of interest. Using examples from experimental neuroscience, we show that batch effects can dominate the signal of interest; and that the choice of normalization method affects the power and reproducibility of the results. While commonly used global normalization methods are not able to adequately normalize the data, more recently developed RNA-seq normalization can. We focus on one particular method, RUVSeq and show that it is able to increase power and biological insight of the results. Finally, we provide a tutorial outlining the implementation of RUVSeq normalization that is applicable to a broad range of studies as well as meta-analysis of publicly available data.

Removing unwanted variation in a differential methylation analysis of Illumina HumanMethylation450 array data

Due to their relatively low-cost per sample and broad, gene-centric coverage of CpGs across the human genome, Illumina's 450k arrays are widely used in large scale differential methylation studies. However, by their very nature, large studies are particularly susceptible to the effects of unwanted variation. The effects of unwanted variation have been extensively documented in gene expression array studies and numerous methods have been developed to mitigate these effects. However, there has been much less research focused on the appropriate methodology to use for accounting for unwanted variation in methylation array studies. Here we present a novel 2-stage approach using RUV-inverse in a differential methylation analysis of 450k data and show that it outperforms existing methods.