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.

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.

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.

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-30). 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.

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.

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.

Subtypes of Barrett's oesophagus and oesophageal adenocarcinoma based on genome-wide methylation analysis.

OBJECTIVE: To identify and characterise DNA methylation subtypes in oesophageal adenocarcinoma (EAC) and its precursor Barrett's oesophagus (BE). DESIGN: We performed genome-wide DNA methylation profiling on samples of non-dysplastic BE from cancer-free patients (n=59), EAC (n=23), normal squamous oesophagus (n=33) and normal fundus (n=9), and identified methylation subtypes using a recursively partitioned mixture model. We assessed genomic alterations for 9 BE and 22 EAC samples with massively parallel sequencing of 243 EAC-associated genes, and we conducted integrative analyses with transcriptome data to identify epigenetically repressed genes. We also carried out in vitro experiments treating EAC cell lines with 5-Aza-2'-Deoxycytidine (5-Aza-dC), short hairpin RNA knockdown and anticancer therapies. RESULTS: We identified and validated four methylation subtypes of EAC and BE. The high methylator subtype (HM) of EAC had the greatest number of activating events in ERBB2 (p<0.05, Student's t-test) and the highest global mutation load (p<0.05, Fisher's exact test). PTPN13 was silenced by aberrant methylation in the HM subtype preferentially and in 57% of EACs overall. In EAC cell lines, 5-Aza-dC treatment restored PTPN13 expression and significantly decreased its promoter methylation in HM cell lines (p<0.05, Welch's t-test). Inhibition of PTPN13 expression in the SK-GT-4 EAC cell line promoted proliferation, colony formation and migration, and increased phosphorylation in ERBB2/EGFR/Src kinase pathways. Finally, EAC cell lines showed subtype-specific responses to topotecan, SN-38 and palbociclib treatment. CONCLUSIONS: We identified and characterised methylator subtypes in BE and EAC. We further demonstrated the biological and clinical relevance of EAC methylator subtypes, which may ultimately help guide clinical management of patients with EAC.

Discovery of methylated circulating DNA biomarkers for comprehensive non-invasive monitoring of treatment response in metastatic colorectal cancer.

OBJECTIVE: Mutations in cell-free circulating DNA (cfDNA) have been studied for tracking disease relapse in colorectal cancer (CRC). This approach requires personalised assay design due to the lack of universally mutated genes. In contrast, early methylation alterations are restricted to defined genomic loci allowing comprehensive assay design for population studies. Our objective was to identify cancer-specific methylated biomarkers which could be measured longitudinally in cfDNA (liquid biopsy) to monitor therapeutic outcome in patients with metastatic CRC (mCRC). DESIGN: Genome-wide methylation microarrays of CRC cell lines (n=149) identified five cancer-specific methylated loci (EYA4, GRIA4, ITGA4, MAP3K14-AS1, MSC). Digital PCR assays were employed to measure methylation of these genes in tumour tissue DNA (n=82) and cfDNA from patients with mCRC (n=182). Plasma longitudinal assessment was performed in a patient subset treated with chemotherapy or targeted therapy. RESULTS: Methylation in at least one marker was detected in all tumour tissue samples and in 156 mCRC patient cfDNA samples (85.7%). Plasma marker prevalence was 71.4% for EYA4, 68.5% for GRIA4, 69.7% for ITGA4, 69.1% for MAP3K14-AS1% and 65.1% for MSC. Dynamics of methylation markers was not affected by treatment type and correlated with objective tumour response and progression-free survival. CONCLUSION: This five-gene methylation panel can be used to circumvent the absence of patient-specific mutations for monitoring tumour burden dynamics in liquid biopsy under different therapeutic regimens. This method might be proposed for assessing pharmacodynamics in clinical trials or when conventional imaging has limitations.

Identification of a key role of widespread epigenetic drift in Barrett's esophagus and esophageal adenocarcinoma.

BACKGROUND: Recent studies have identified age-related changes in DNA methylation patterns in normal and cancer tissues in a process that is called epigenetic drift. However, the evolving patterns, functional consequences, and dynamics of epigenetic drift during carcinogenesis remain largely unexplored. Here we analyze the evolution of epigenetic drift patterns during progression from normal squamous esophagus tissue to Barrett's esophagus (BE) to esophageal adenocarcinoma (EAC) using 173 tissue samples from 100 (nonfamilial) BE patients, along with publically available datasets including The Cancer Genome Atlas (TCGA). RESULTS: Our analysis reveals extensive methylomic drift between normal squamous esophagus and BE tissues in nonprogressed BE patients, with differential drift affecting 4024 (24%) of 16,984 normally hypomethylated cytosine-guanine dinucleotides (CpGs) occurring in CpG islands. The majority (63%) of islands that include drift CpGs are associated with gene promoter regions. Island CpGs that drift have stronger pairwise correlations than static islands, reflecting collective drift consistent with processive DNA methylation maintenance. Individual BE tissues are extremely heterogeneous in their distribution of methylomic drift and encompass unimodal low-drift to bimodal high-drift patterns, reflective of differences in BE tissue age. Further analysis of longitudinally collected biopsy samples from 20 BE patients confirm the time-dependent evolution of these drift patterns. Drift patterns in EAC are similar to those in BE, but frequently exhibit enhanced bimodality and advanced mode drift. To better understand the observed drift patterns, we developed a multicellular stochastic model at the CpG island level. Importantly, we find that nonlinear feedback in the model between mean island methylation and CpG methylation rates is able to explain the widely heterogeneous collective drift patterns. Using matched gene expression and DNA methylation data in EAC from TCGA and other publically available data, we also find that advanced methylomic drift is correlated with significant transcriptional repression of ~ 200 genes in important regulatory and developmental pathways, including several checkpoint and tumor suppressor-like genes. CONCLUSIONS: Taken together, our findings suggest that epigenetic drift evolution acts to significantly reduce the expression of developmental genes that may alter tissue characteristics and improve functional adaptation during BE to EAC progression.

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.