Benchmark of cellular deconvolution methods using a multi-assay reference dataset from postmortem human prefrontal cortex

Background

Cellular deconvolution of bulk RNA-sequencing (RNA-seq) data using single cell or nuclei RNA-seq (sc/snRNA-seq) reference data is an important strategy for estimating cell type composition in heterogeneous tissues, such as human brain. Computational methods for deconvolution have been developed and benchmarked against simulated data, pseudobulked sc/snRNA-seq data, or immunohistochemistry reference data. A major limitation in developing improved deconvolution algorithms has been the lack of integrated datasets with orthogonal measurements of gene expression and estimates of cell type proportions on the same tissue sample. Deconvolution algorithm performance has not yet been evaluated across different RNA extraction methods (cytosolic, nuclear, or whole cell RNA), different library preparation types (mRNA enrichment vs. ribosomal RNA depletion), or with matched single cell reference datasets.

Results

A rich multi-assay dataset was generated in postmortem human dorsolateral prefrontal cortex (DLPFC) from 22 tissue blocks. Assays included spatially-resolved transcriptomics, snRNA-seq, bulk RNA-seq (across six library/extraction RNA-seq combinations), and RNAScope/Immunofluorescence (RNAScope/IF) for six broad cell types. The Mean Ratio method, implemented in the DeconvoBuddies R package, was developed for selecting cell type marker genes. Six computational deconvolution algorithms were evaluated in DLPFC and predicted cell type proportions were compared to orthogonal RNAScope/IF measurements.

Conclusions

Bisque and hspe were the most accurate methods, were robust to differences in RNA library types and extractions. This multi-assay dataset showed that cell size differences, marker genes differentially quantified across RNA libraries, and cell composition variability in reference snRNA-seq impact the accuracy of current deconvolution methods.

lute: estimating the cell composition of heterogeneous tissue with varying cell sizes using gene expression

Relative cell type fraction estimates in bulk RNA-sequencing data are important to control for cell composition differences across heterogenous tissue samples. Current computational tools estimate relative RNA abundances rather than cell type proportions in tissues with varying cell sizes, leading to biased estimates. We present lute, a computational tool to accurately deconvolute cell types with varying sizes. Our software wraps existing deconvolution algorithms in a standardized framework. Using simulated and real datasets, we demonstrate how lute adjusts for differences in cell sizes to improve the accuracy of cell composition. Software is available from https://bioconductor.org/packages/lute.

Challenges and opportunities to computationally deconvolve heterogeneous tissue with varying cell sizes using single-cell RNA-sequencing datasets

Deconvolution of cell mixtures in "bulk" transcriptomic samples from homogenate human tissue is important for understanding disease pathologies. However, several experimental and computational challenges impede transcriptomics-based deconvolution approaches using single-cell/nucleus RNA-seq reference atlases. Cells from the brain and blood have substantially different sizes, total mRNA, and transcriptional activities, and existing approaches may quantify total mRNA instead of cell type proportions. Further, standards are lacking for the use of cell reference atlases and integrative analyses of single-cell and spatial transcriptomics data. We discuss how to approach these key challenges with orthogonal "gold standard" datasets for evaluating deconvolution methods.

Challenges and opportunities to computationally deconvolve heterogeneous tissue with varying cell sizes using single cell RNA-sequencing datasets

Deconvolution of cell mixtures in "bulk" transcriptomic samples from homogenate human tissue is important for understanding the pathologies of diseases. However, several experimental and computational challenges remain in developing and implementing transcriptomics-based deconvolution approaches, especially those using a single cell/nuclei RNA-seq reference atlas, which are becoming rapidly available across many tissues. Notably, deconvolution algorithms are frequently developed using samples from tissues with similar cell sizes. However, brain tissue or immune cell populations have cell types with substantially different cell sizes, total mRNA expression, and transcriptional activity. When existing deconvolution approaches are applied to these tissues, these systematic differences in cell sizes and transcriptomic activity confound accurate cell proportion estimates and instead may quantify total mRNA content. Furthermore, there is a lack of standard reference atlases and computational approaches to facilitate integrative analyses, including not only bulk and single cell/nuclei RNA-seq data, but also new data modalities from spatial -omic or imaging approaches. New multi-assay datasets need to be collected with orthogonal data types generated from the same tissue block and the same individual, to serve as a "gold standard" for evaluating new and existing deconvolution methods. Below, we discuss these key challenges and how they can be addressed with the acquisition of new datasets and approaches to analysis.

recountmethylation enables flexible analysis of public blood DNA methylation array data

Summary

Thousands of DNA methylation (DNAm) array samples from human blood are publicly available on the Gene Expression Omnibus (GEO), but they remain underutilized for experiment planning, replication and cross-study and cross-platform analyses. To facilitate these tasks, we augmented our recountmethylation R/Bioconductor package with 12 537 uniformly processed EPIC and HM450K blood samples on GEO as well as several new features. We subsequently used our updated package in several illustrative analyses, finding (i) study ID bias adjustment increased variation explained by biological and demographic variables, (ii) most variation in autosomal DNAm was explained by genetic ancestry and CD4+ T-cell fractions and (iii) the dependence of power to detect differential methylation on sample size was similar for each of peripheral blood mononuclear cells (PBMC), whole blood and umbilical cord blood. Finally, we used PBMC and whole blood to perform independent validations, and we recovered 38-46% of differentially methylated probes between sexes from two previously published epigenome-wide association studies.

Availability and implementation

Source code to reproduce the main results are available on GitHub (repo: recountmethylation_flexible-blood-analysis_manuscript; url: https://github.com/metamaden/recountmethylation_flexible-blood-analysis_manuscript). All data was publicly available and downloaded from the Gene Expression Omnibus (https://www.ncbi.nlm.nih.gov/geo/). Compilations of the analyzed public data can be accessed from the website recount.bio/data (preprocessed HM450K array data: https://recount.bio/data/remethdb_h5se-gm_epic_0-0-2_1589820348/; preprocessed EPIC array data: https://recount.bio/data/remethdb_h5se-gm_epic_0-0-2_1589820348/).

Supplementary information

Supplementary data are available at Bioinformatics Advances online.

Retained introns in long RNA-seq reads are not reliably detected in sample-matched short reads

BACKGROUND

There is growing interest in retained introns in a variety of disease contexts including cancer and aging. Many software tools have been developed to detect retained introns from short RNA-seq reads, but reliable detection is complicated by overlapping genes and transcripts as well as the presence of unprocessed or partially processed RNAs.

RESULTS

We compared introns detected by 8 tools using short RNA-seq reads with introns observed in long RNA-seq reads from the same biological specimens. We found significant disagreement among tools (Fleiss' [Formula: see text]) such that 47.7% of all detected intron retentions were not called by more than one tool. We also observed poor performance of all tools, with none achieving an F1-score greater than 0.26, and qualitatively different behaviors between general-purpose alternative splicing detection tools and tools confined to retained intron detection.

CONCLUSIONS

Short-read tools detect intron retention with poor recall and precision, calling into question the completeness and validity of a large percentage of putatively retained introns called by commonly used methods.

Human methylome variation across Infinium 450K data on the Gene Expression Omnibus

While DNA methylation (DNAm) is the most-studied epigenetic mark, few recent studies probe the breadth of publicly available DNAm array samples. We collectively analyzed 35 360 Illumina Infinium HumanMethylation450K DNAm array samples published on the Gene Expression Omnibus. We learned a controlled vocabulary of sample labels by applying regular expressions to metadata and used existing models to predict various sample properties including epigenetic age. We found approximately two-thirds of samples were from blood, one-quarter were from brain and one-third were from cancer patients. About 19% of samples failed at least one of Illumina's 17 prescribed quality assessments; signal distributions across samples suggest modifying manufacturer-recommended thresholds for failure would make these assessments more informative. We further analyzed DNAm variances in seven tissues (adipose, nasal, blood, brain, buccal, sperm and liver) and characterized specific probes distinguishing them. Finally, we compiled DNAm array data and metadata, including our learned and predicted sample labels, into database files accessible via the recountmethylation R/Bioconductor companion package. Its vignettes walk the user through some analyses contained in this paper.

Human Leukocyte Antigen Susceptibility Map for Severe Acute Respiratory Syndrome Coronavirus 2

Genetic variability across the three major histocompatibility complex (MHC) class I genes (human leukocyte antigen A [HLA-A], -B, and -C genes) may affect susceptibility to and severity of the disease caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), the virus responsible for coronavirus disease 2019 (COVID-19). We performed a comprehensive in silico analysis of viral peptide-MHC class I binding affinity across 145 HLA-A, -B, and -C genotypes for all SARS-CoV-2 peptides. We further explored the potential for cross-protective immunity conferred by prior exposure to four common human coronaviruses. The SARS-CoV-2 proteome was successfully sampled and was represented by a diversity of HLA alleles. However, we found that HLA-B*46:01 had the fewest predicted binding peptides for SARS-CoV-2, suggesting that individuals with this allele may be particularly vulnerable to COVID-19, as they were previously shown to be for SARS (M. Lin, H.-T. Tseng, J. A. Trejaut, H.-L. Lee, et al., BMC Med Genet 4:9, 2003, https://bmcmedgenet.biomedcentral.com/articles/10.1186/1471-2350-4-9). Conversely, we found that HLA-B*15:03 showed the greatest capacity to present highly conserved SARS-CoV-2 peptides that are shared among common human coronaviruses, suggesting that it could enable cross-protective T-cell-based immunity. Finally, we reported global distributions of HLA types with potential epidemiological ramifications in the setting of the current pandemic.IMPORTANCE Individual genetic variation may help to explain different immune responses to a virus across a population. In particular, understanding how variation in HLA may affect the course of COVID-19 could help identify individuals at higher risk from the disease. HLA typing can be fast and inexpensive. Pairing HLA typing with COVID-19 testing where feasible could improve assessment of severity of viral disease in the population. Following the development of a vaccine against SARS-CoV-2, the virus that causes COVID-19, individuals with high-risk HLA types could be prioritized for vaccination.

Putatively cancer-specific exon-exon junctions are shared across patients and present in developmental and other non-cancer cells

This study probes the distribution of putatively cancer-specific junctions across a broad set of publicly available non-cancer human RNA sequencing (RNA-seq) datasets. We compared cancer and non-cancer RNA-seq data from The Cancer Genome Atlas (TCGA), the Genotype-Tissue Expression (GTEx) Project and the Sequence Read Archive. We found that (i) averaging across cancer types, 80.6% of exon–exon junctions thought to be cancer-specific based on comparison with tissue-matched samples (σ = 13.0%) are in fact present in other adult non-cancer tissues throughout the body; (ii) 30.8% of junctions not present in any GTEx or TCGA normal tissues are shared by multiple samples within at least one cancer type cohort, and 87.4% of these distinguish between different cancer types; and (iii) many of these junctions not found in GTEx or TCGA normal tissues (15.4% on average, σ = 2.4%) are also found in embryological and other developmentally associated cells. These findings refine the meaning of RNA splicing event novelty, particularly with respect to the human neoepitope repertoire. Ultimately, cancer-specific exon–exon junctions may have a substantial causal relationship with the biology of disease.

Dysfunctional epigenetic aging of the normal colon and colorectal cancer risk

Background

Chronological age is a prominent risk factor for many types of cancers including colorectal cancer (CRC). Yet, the risk of CRC varies substantially between individuals, even within the same age group, which may reflect heterogeneity in biological tissue aging between people. Epigenetic clocks based on DNA methylation are a useful measure of the biological aging process with the potential to serve as a biomarker of an individual’s susceptibility to age-related diseases such as CRC.

Methods

We conducted a genome-wide DNA methylation study on samples of normal colon mucosa (N = 334). Subjects were assigned to three cancer risk groups (low, medium, and high) based on their personal adenoma or cancer history. Using previously established epigenetic clocks (Hannum, Horvath, PhenoAge, and EpiTOC), we estimated the biological age of each sample and assessed for epigenetic age acceleration in the samples by regressing the estimated biological age on the individual’s chronological age. We compared the epigenetic age acceleration between different risk groups using a multivariate linear regression model with the adjustment for gender and cell-type fractions for each epigenetic clock. An epigenome-wide association study (EWAS) was performed to identify differential methylation changes associated with CRC risk.

Results

Each epigenetic clock was significantly correlated with the chronological age of the subjects, and the Horvath clock exhibited the strongest correlation in all risk groups (r > 0.8, p < 1 × 10⁻³⁰). The PhenoAge clock (p = 0.0012) revealed epigenetic age deceleration in the high-risk group compared to the low-risk group.

Conclusions

Among the four DNA methylation-based measures of biological age, the Horvath clock is the most accurate for estimating the chronological age of individuals. Individuals with a high risk for CRC have epigenetic age deceleration in their normal colons measured by the PhenoAge clock, which may reflect a dysfunctional epigenetic aging process.

Senescence-associated tissue microenvironment promotes colon cancer formation through the secretory factor GDF15

The risk of colorectal cancer (CRC) varies between people, and the cellular mechanisms mediating the differences in risk are largely unknown. Senescence has been implicated as a causative cellular mechanism for many diseases, including cancer, and may affect the risk for CRC. Senescent fibroblasts that accumulate in tissues secondary to aging and oxidative stress have been shown to promote cancer formation via a senescence‐associated secretory phenotype (SASP). In this study, we assessed the role of senescence and the SASP in CRC formation. Using primary human colon tissue, we found an accumulation of senescent fibroblasts in normal tissues from individuals with advanced adenomas or carcinomas in comparison with individuals with no polyps or CRC. In in vitro and ex vivo model systems, we induced senescence using oxidative stress in colon fibroblasts and demonstrated that the senescent fibroblasts secrete GDF15 as an essential SASP factor that promotes cell proliferation, migration, and invasion in colon adenoma and CRC cell lines as well as primary colon organoids via the MAPK and PI3K signaling pathways. In addition, we observed increased mRNA expression of GDF15 in primary normal colon tissue from people at increased risk for CRC in comparison with average risk individuals. These findings implicate the importance of a senescence‐associated tissue microenvironment and the secretory factor GDF15 in promoting CRC formation.

Implications of Epigenetic Drift in Colorectal Neoplasia

Many normal tissues undergo age-related drift in DNA methylation, providing a quantitative measure of tissue age. Here, we identify and validate 781 CpG islands (CGI) that undergo significant methylomic drift in 232 normal colorectal tissues and show that these CGI continue to drift in neoplasia while retaining significant correlations across samples. However, compared with normal colon, this drift advanced (∼3-4-fold) faster in neoplasia, consistent with increased cell proliferation during neoplastic progression. The observed drift patterns were broadly consistent with modeled adenoma-to-carcinoma sojourn time distributions from colorectal cancer incidence data. These results support the hypothesis that, beginning with the founder premalignant cell, cancer precursors frequently sojourn for decades before turning into cancer, implying that the founder cell typically arises early in life. At least 77% to 89% of the observed drift variance in distal and rectal tumors was explained by stochastic variability associated with neoplastic progression, whereas only 55% of the variance was explained for proximal tumors. However, gene-CGI pairs in the proximal colon that underwent drift were significantly and primarily negatively correlated with cancer gene expression, suggesting that methylomic drift participates in the clonal evolution of colorectal cancer. Methylomic drift advanced in colorectal neoplasia, consistent with extended sojourn time distributions, which accounts for a significant fraction of epigenetic heterogeneity in colorectal cancer. Importantly, these estimated long-duration premalignant sojourn times suggest that early dietary and lifestyle interventions may be more effective than later changes in reducing colorectal cancer incidence. SIGNIFICANCE: These findings present age-related methylomic drift in colorectal neoplasia as evidence that premalignant cells can persist for decades before becoming cancerous.See related commentary by Sapienza, p. 437.

Dietary fatty acid intake is associated with paraoxonase 1 activity in a cohort-based analysis of 1,548 subjects

BACKGROUND

Paraoxonase 1 (PON1) is a cardioprotective, HDL-associated glycoprotein enzyme with broad substrate specificity. Our previous work found associations between dietary cholesterol and vitamin C with PON1 activity. The goal of this study was to determine the effect of specific dietary fatty acid (DFA) intake on PON1 activity.

METHODS

1,548 participants with paraoxonase activity measures completed the Harvard Standardized Food Frequency Questionnaire to determine their daily nutrient intake over the past year. Eight saturated, 3 monounsaturated, and 6 polyunsaturated DFAs were measured by the questionnaire. To reduce the number of observations tested, only specific fatty acids that were not highly correlated (r < 0.8) with other DFAs or that were representative of other DFAs through high correlation within each respective group (saturated, monounsaturated, or polyunsaturated) were retained for analysis. Six specific DFA intakes - myristic acid (14 carbon atoms, no double bonds - 14:0), oleic acid (18:1), gadoleic acid (20:1), α-linolenic acid (18:3), arachidonic acid (20:4), and eicosapentaenoic acid (20:5) - were carried forward to stepwise linear regression, which evaluated the effect of each specific DFA on covariate-adjusted PON1 enzyme activity.

RESULTS

Four of the 6 tested DFA intakes - myristic acid (p = 0.038), gadoleic acid (p = 6.68 × 10(-7)), arachidonic acid (p = 0.0007), and eicosapentaenoic acid (p = 0.013) - were independently associated with covariate-adjusted PON1 enzyme activity. Myristic acid, a saturated fat, and gadoleic acid, a monounsaturated fat, were both positively associated with PON1 activity. Both of the tested polyunsaturated fats, arachidonic acid and eicosapentaenoic acid, were negatively associated with PON1 activity.

CONCLUSIONS

This study presents the largest cohort-based analysis of the relationship between dietary lipids and PON1 enzyme activity. Further research is necessary to elucidate and understand the specific biological mechanisms, whether direct or regulatory, through which DFAs affect PON1 activity.