Search

Significant variations have been observed in viral copies generated during SARS-CoV-2 infections. However, the factors that impact viral copies and infection dynamics are not fully understood, and may be inherently dependent upon different viral and host factors. Here, we conducted virus whole genome sequencing and measured viral copies using RT-qPCR from 9,902 SARS-CoV-2 infections over a 2-year period to examine the impact of virus genetic variation on changes in viral copies adjusted for host age and vaccination status. Using a genome-wide association study (GWAS) approach, we identified multiple single-nucleotide polymorphisms (SNPs) corresponding to amino acid changes in the SARS-CoV-2 genome associated with variations in viral copies. We further applied a marginal epistasis test to detect interactions among SNPs and identified multiple pairs of substitutions located in the spike gene that have non-linear effects on viral copies. We also analyzed the temporal patterns and found that SNPs associated with increased viral copies were predominantly observed in Delta and Omicron BA.2/BA.4/BA.5/XBB infections, whereas those associated with decreased viral copies were only observed in infections with Omicron BA.1 variants. Our work showcases how GWAS can be a useful tool for probing phenotypes related to SNPs in viral genomes that are worth further exploration. We argue that this approach can be used more broadly across pathogens to characterize emerging variants and monitor therapeutic interventions. Author summary: Our study explores why viral load (copies measured by RT-qPCR) varies during SARS-CoV-2 infections by analyzing viral mutations and measuring viral copies in 9,902 individuals over two years. We aimed to understand how genetic differences in SARS-CoV-2 influence viral copies, considering host age and vaccination status. Using a genome-wide association study (GWAS), we identified several single-nucleotide polymorphisms (SNPs) in the virus linked to variations in viral levels. Notably, interactions between certain SNPs in the spike gene had non-linear effects on viral copies. Our analysis revealed that SNPs associated with higher viral copies were common in Delta and Omicron BA.2/BA.4/BA.5/XBB variants, while those linked to lower levels were mainly found in Omicron BA.1. This research highlights GWAS as a powerful tool for exploring virus genetics and suggests it can be broadly applied to monitor new variants of COVID-19 and other infectious diseases. [ABSTRACT FROM AUTHOR]

In addition to the classical HLA genes, the major histocompatibility complex (MHC) harbors a high number of other polymorphic genes with less established roles in disease associations and transplantation matching. To facilitate studies of the non-classical and non-HLA genes in large patient and biobank cohorts, we trained imputation models for MICA, MICB, HLA-E, HLA-F and HLA-G alleles on genome SNP array data. We show, using both population-specific and multi-population 1000 Genomes references, that the alleles of these genes can be accurately imputed for screening and research purposes. The best imputation model for MICA, MICB, HLA-E, -F and -G achieved a mean accuracy of 99.3% (min, max: 98.6, 99.9). Furthermore, validation of the 1000 Genomes exome short-read sequencing-based allele calling against a clinical-grade reference data showed an average accuracy of 99.8%, testifying for the quality of the 1000 Genomes data as an imputation reference. We also fitted the models for Infinium Global Screening Array (GSA, Illumina, Inc.) and Axiom Precision Medicine Research Array (PMRA, Thermo Fisher Scientific Inc.) SNP content, with mean accuracies of 99.1% (97.2, 100) and 98.9% (97.4, 100), respectively. Author summary: The major histocompatibility complex (MHC) region on chromosome 6 significantly influences disease risk, particularly in autoimmune conditions. To improve fine-mapping of potentially causal genetic variants within this region, we developed accurate imputation methods for inferring functional allelic variation from MHC SNP data. While existing tools primarily focus on classical HLA genes, our study extends to non-classical HLA genes (HLA-E, HLA-F, and HLA-G) and MHC Class I Chain-Related MIC genes (MICA and MICB) which have specific functions in innate and adaptive immunity. Leveraging population-specific Finnish and multi-population 1000 Genomes references, our imputation models demonstrate high accuracy. Moreover, we tailored models for two widely used genome SNP arrays: the Infinium Global Screening Array (Illumina, Inc.) and the Axiom Precision Medicine Research Array (Thermo Fisher Scientific Inc.). These freely available, multi-population models empower researchers to explore genetic MHC associations in more detail and contribute to our understanding of immune-related disease mechanisms. [ABSTRACT FROM AUTHOR]

Human induced pluripotent stem cell-derived cardiomyocytes (iPSC-CMs) have gained traction as a powerful model in cardiac disease and therapeutics research, since iPSCs are self-renewing and can be derived from healthy and diseased patients without invasive surgery. However, current iPSC-CM differentiation methods produce cardiomyocytes with immature, fetal-like electrophysiological phenotypes, and the variety of maturation protocols in the literature results in phenotypic differences between labs. Heterogeneity of iPSC donor genetic backgrounds contributes to additional phenotypic variability. Several mathematical models of iPSC-CM electrophysiology have been developed to help to predict cell responses, but these models individually do not capture the phenotypic variability observed in iPSC-CMs. Here, we tackle these limitations by developing a computational pipeline to calibrate cell preparation-specific iPSC-CM electrophysiological parameters. We used the genetic algorithm (GA), a heuristic parameter calibration method, to tune ion channel parameters in a mathematical model of iPSC-CM physiology. To systematically optimize an experimental protocol that generates sufficient data for parameter calibration, we created in silico datasets by simulating various protocols applied to a population of models with known conductance variations, and then fitted parameters to those datasets. We found that calibrating to voltage and calcium transient data under 3 varied experimental conditions, including electrical pacing combined with ion channel blockade and changing buffer ion concentrations, improved model parameter estimates and model predictions of unseen channel block responses. This observation also held when the fitted data were normalized, suggesting that normalized fluorescence recordings, which are more accessible and higher throughput than patch clamp recordings, could sufficiently inform conductance parameters. Therefore, this computational pipeline can be applied to different iPSC-CM preparations to determine cell line-specific ion channel properties and understand the mechanisms behind variability in perturbation responses. Author summary: Many drug treatments or environmental factors can trigger cardiac arrhythmias, which are dangerous and often unpredictable. Human cardiomyocytes derived from donor stem cells have proven to be a promising model for studying these events, but variability in donor genetic background and cell maturation methods, as well as overall immaturity of stem cell-derived cardiomyocytes relative to the adult heart, have hindered reproducibility and reliability of these studies. Mathematical models of these cells can aid in understanding the underlying electrophysiological contributors to this variability, but determining these models' parameters for multiple cell preparations is challenging. In this study, we tackle these limitations by developing a computational method to simultaneously estimate multiple model parameters using data from imaging-based experiments, which can be easily scaled to rapidly characterize multiple cell lines. This method can generate many personalized models of individual cell preparations, improving drug response predictions and revealing specific differences in electrophysiological properties that contribute to variability in cardiac maturity and arrhythmia susceptibility. [ABSTRACT FROM AUTHOR]

Alzheimer's disease (AD) involves aggregation of amyloid β and tau, neuron loss, cognitive decline, and neuroinflammatory responses. Both resident microglia and peripheral immune cells have been associated with the immune component of AD. However, the relative contribution of resident and peripheral immune cell types to AD predisposition has not been thoroughly explored due to their similarity in gene expression and function. To study the effects of AD associated variants on cis-regulatory elements, we train convolutional neural network (CNN) regression models that link genome sequence to cell type-specific levels of open chromatin, a proxy for regulatory element activity. We then use in silico mutagenesis of regulatory sequences to predict the relative impact of candidate variants across these cell types. We develop and apply criteria for evaluating our models and refine our models using massively parallel reporter assay (MPRA) data. Our models identify multiple AD-associated variants with a greater predicted impact in peripheral cells relative to microglia or neurons. Our results support their use as models to study the effects of AD-associated variants and even suggest that peripheral immune cells themselves may mediate a component of AD predisposition. We make our library of CNN models and predictions available as a resource for the community to study immune and neurological disorders. Author summary: Human genetic association studies continue to identify a large number of genetic variants associated with common diseases. However, genetic linkage renders it difficult to identify truly causal variants from these associations. The problem is further complicated by the fact that most disease-associated variants lie in non-coding regions of the genome which are involved in regulating gene expression. Gene regulatory elements are highly tissue- and cell-type-specific which makes it important to study disease variant effects in the right tissue or cell type. In this study, we propose a computational framework to identify functional variants that influence disease predisposition by disruption of gene regulatory element activity. We first identify potential cell types that may influence disease predisposition through cell-type-specific gene regulatory element activity. Then, we use machine learning (ML) models that relate human genetic sequence to cell-type-specific gene regulatory activity to characterize potential effects of genetic variants. Finally, using in silico mutagenesis, we identify putative functional variants and predict their ability to disrupt regulatory activity in various cell types. Further, we define best practices to train, benchmark, improve and apply these models to interpret the cell-type-specific determinants of disease. In a specific use case, we apply this framework to interpret the cell type-specific contributions of regulatory genetic variants in Alzheimer's Disease (AD). For AD, we show that we can separate out peripheral immune and brain resident immune microglial signals to understand their specific contributions to AD predisposition. We find that gene regulatory elements active in peripheral immune cells strongly colocalize with AD genetic variants compared to regulatory elements active in neural cell types. In addition, our models predict that many functional AD variants have variant effects that are specific to peripheral immune cells and not brain resident immune cells (microglia). Our study opens up opportunities to experimentally study Alzheimer's Disease variants in peripheral immune cell models that are often easier to manipulate than many neural cell types and yet not usually considered in functional studies on AD variants. [ABSTRACT FROM AUTHOR]

Cancer is a significant global health issue, with treatment challenges arising from intratumor heterogeneity. This heterogeneity stems mainly from somatic evolution, causing genetic diversity within the tumor, and phenotypic plasticity of tumor cells leading to reversible phenotypic changes. However, the interplay of both factors has not been rigorously investigated. Here, we examine the complex relationship between somatic evolution and phenotypic plasticity, explicitly focusing on the interplay between cell migration and proliferation. This type of phenotypic plasticity is essential in glioblastoma, the most aggressive form of brain tumor. We propose that somatic evolution alters the regulation of phenotypic plasticity in tumor cells, specifically the reaction to changes in the microenvironment. We study this hypothesis using a novel, spatially explicit model that tracks individual cells' phenotypic and genetic states. We assume cells change between migratory and proliferative states controlled by inherited and mutation-driven genotypes and the cells' microenvironment. We observe that cells at the tumor edge evolve to favor migration over proliferation and vice versa in the tumor bulk. Notably, different genetic configurations can result in this pattern of phenotypic heterogeneity. We analytically predict the outcome of the evolutionary process, showing that it depends on the tumor microenvironment. Synthetic tumors display varying levels of genetic and phenotypic heterogeneity, which we show are predictors of tumor recurrence time after treatment. Interestingly, higher phenotypic heterogeneity predicts poor treatment outcomes, unlike genetic heterogeneity. Our research offers a novel explanation for heterogeneous patterns of tumor recurrence in glioblastoma patients. Author summary: Intratumor heterogeneity presents a significant barrier to effective cancer therapy. This heterogeneity stems from the evolution of cancer cells and their capability for phenotypic plasticity. However, the interplay between these two factors still needs to be fully understood. This study examines the interaction between cancer cell evolution and phenotypic plasticity, focusing on the phenotypic switch between migration and proliferation. Such plasticity is particularly relevant to glioblastoma, the most aggressive form of brain tumor. By employing a novel model, we explore how tumor cell evolution, influenced by both genotype and microenvironment, contributes to tumor heterogeneity. We observe that cells at the tumor periphery tend to migrate, while those within the tumor are more inclined to proliferate. Interestingly, our analysis reveals that distinct genetic configurations of the tumor can lead to this observed pattern. Further, we delve into the implications for cancer treatment and discover that it is phenotypic, rather than genetic, heterogeneity that more accurately predicts tumor recurrence following therapy. Our findings offer insights into the significant variability observed in glioblastoma recurrence times post-treatment. [ABSTRACT FROM AUTHOR]

For decades, the 16S rRNA gene has been used to taxonomically classify prokaryotic species and to taxonomically profile microbial communities. However, the 16S rRNA gene has been criticized for being too conserved to differentiate between distinct species. We argue that the inability to differentiate between species is not a unique feature of the 16S rRNA gene. Rather, we observe the gradual loss of species-level resolution for other nearly-universal prokaryotic marker genes as the number of gene sequences increases in reference databases. This trend was strongly correlated with how represented a taxonomic group was in the database and indicates that, at the gene-level, the boundaries between many species might be fuzzy. Through our study, we argue that any approach that relies on a single marker to distinguish bacterial taxa is fraught even if some markers appear to be discriminative in current databases. Author summary: The use of reference databases for assigning taxonomic labels to genomic and metagenomic sequences is a fundamental bioinformatic task in the characterization of microbial communities. The increasing accessibility of high throughput sequencing has led to a rapid increase in the size and number of sequences in databases. This has been beneficial for improving our understanding of the global microbial genetic diversity. However, there is evidence that as the microbial diversity is more densely sampled, increasingly longer genomic segments are needed to differentiate between distinct bacterial species. The scientific community needs to be aware of this issue and needs to develop methods that better account for it when assigning taxonomic labels to metagenomic sequences from microbial communities. [ABSTRACT FROM AUTHOR]

A core problem in genetics is molecular quantitative trait locus (QTL) mapping, in which genetic variants associated with changes in the molecular phenotypes are identified. One of the most-studied molecular QTL mapping problems is expression QTL (eQTL) mapping, in which the molecular phenotype is gene expression. It is common in eQTL mapping to compute gene expression by aggregating the expression levels of individual isoforms from the same gene and then performing linear regression between SNPs and this aggregated gene expression level. However, SNPs may regulate isoforms from the same gene in different directions due to alternative splicing, or only regulate the expression level of one isoform, causing this approach to lose power. Here, we examine a broader question: which genes have at least one isoform whose expression level is regulated by genetic variants? In this study, we propose and evaluate several approaches to answering this question, demonstrating that "isoform-aware" methods—those that account for the expression levels of individual isoforms—have substantially greater power to answer this question than standard "gene-level" eQTL mapping methods. We identify settings in which different approaches yield an inflated number of false discoveries or lose power. In particular, we show that calling an eGene if there is a significant association between a SNP and any isoform fails to control False Discovery Rate, even when applying standard False Discovery Rate correction. We show that similar trends are observed in real data from the GEUVADIS and GTEx studies, suggesting the possibility that similar effects are present in these consortia. [ABSTRACT FROM AUTHOR]

Genome-wide association studies (GWAS) have successfully identified over two hundred thousand genotype-trait associations. Yet some challenges remain. First, complex traits are often associated with many single nucleotide polymorphisms (SNPs), most with small or moderate effect sizes, making them difficult to detect. Second, many complex traits share a common genetic basis due to 'pleiotropy' and and though few methods consider it, leveraging pleiotropy can improve statistical power to detect genotype-trait associations with weaker effect sizes. Third, currently available statistical methods are limited in explaining the functional mechanisms through which genetic variants are associated with specific or multiple traits. We propose multi-GPA-Tree to address these challenges. The multi-GPA-Tree approach can identify risk SNPs associated with single as well as multiple traits while also identifying the combinations of functional annotations that can explain the mechanisms through which risk-associated SNPs are linked with the traits. First, we implemented simulation studies to evaluate the proposed multi-GPA-Tree method and compared its performance with existing statistical approaches. The results indicate that multi-GPA-Tree outperforms existing statistical approaches in detecting risk-associated SNPs for multiple traits. Second, we applied multi-GPA-Tree to a systemic lupus erythematosus (SLE) and rheumatoid arthritis (RA), and to a Crohn's disease (CD) and ulcertive colitis (UC) GWAS, and functional annotation data including GenoSkyline and GenoSkylinePlus. Our results demonstrate that multi-GPA-Tree can be a powerful tool that improves association mapping while facilitating understanding of the underlying genetic architecture of complex traits and potential mechanisms linking risk-associated SNPs with complex traits. Author summary: In spite of continued success in developing statistical methodologies that integrate GWAS summary statistics and functional annotation data, existing methods are unable to pinpoint the interactions between functional annotations that influence one or more traits. Hence, the underlying interactions between biological mechanisms linking risk-associated SNPs to traits remain unknown. We propose multi-GPA-Tree to identify risk-associated SNPs and the combinations of functional annotations related to one or more trait risk-associated SNPs. Notably, multi-GPA-Tree requires only GWAS p-value summary statistics, instead of individual level genotype-phenotype data, making it more viable to implement. Compared to the existing state-of-the-art methods, multi-GPA-Tree showed improved performance in simulation studies and validated results for several auto-immune diseases in real data application. These combined results suggest that multi-GPA-Tree is an effective tool for integrative analysis and can potentially be valuable to clinical genomic researchers for hypothesis generation and validation. [ABSTRACT FROM AUTHOR]

Phenotype prediction is at the center of many questions in biology. Prediction is often achieved by determining statistical associations between genetic and phenotypic variation, ignoring the exact processes that cause the phenotype. Here, we present a framework based on genome-scale metabolic reconstructions to reveal the mechanisms behind the associations. We calculated a polygenic score (PGS) that identifies a set of enzymes as predictors of growth, the phenotype. This set arises from the synergy of the functional mode of metabolism in a particular setting and its evolutionary history, and is suitable to infer the phenotype across a variety of conditions. We also find that there is optimal genetic variation for predictability and demonstrate how the linear PGS can still explain phenotypes generated by the underlying nonlinear biochemistry. Therefore, the explicit model interprets the black box statistical associations of the genotype-to-phenotype map and helps to discover what limits the prediction in metabolism. Author summary: Predicting phenotypes from genotypes is fundamental in biology. The prevalent method, quantitative genetics, leverages genetic and phenotypic variation data to create predictive statistical tools like polygenic scores -weighted sums of phenotype-linked alleles. However, this approach overlooks the nonlinear mechanisms that significantly contribute to the genesis of the phenotype. In this work, we used comprehensive genome-wide metabolic models as a direct genotype-to-phenotype map that enables us to scrutinize the mechanistic basis behind statistically forecasting growth (biomass). We introduce a procedure within this framework for generating the necessary genetic and phenotypic diversity essential for calculating polygenic scores and assess the impact of the extent of genetic variation and environmental factors that curtail predictability. Our discoveries underscore that the structure of the polygenic score relies on the interplay among three pivotal components: the functional mode governing phenotype generation, the imprint of evolutionary history on essential parameters, and the genetic architecture at play. Our study exemplifies how integrating prediction and interpretation can yield a broader understanding of how genetic information and phenotypic outcomes interact. This approach transcends linear assumptions, enriching our insight into the complex nonlinear nature of genotype-phenotype relationships. [ABSTRACT FROM AUTHOR]

Lung adenocarcinoma (LUAD) is a deadly tumor with dynamic evolutionary process. Although much endeavors have been made in identifying the temporal patterns of cancer progression, it remains challenging to infer and interpret the molecular alterations associated with cancer development and progression. To this end, we developed a computational approach to infer the progression trajectory based on cross-sectional transcriptomic data. Analysis of the LUAD data using our approach revealed a linear trajectory with three different branches for malignant progression, and the results showed consistency in three independent cohorts. We used the progression model to elucidate the potential molecular events in LUAD progression. Further analysis showed that overexpression of BUB1B, BUB1 and BUB3 promoted tumor cell proliferation and metastases by disturbing the spindle assembly checkpoint (SAC) in the mitosis. Aberrant mitotic spindle checkpoint signaling appeared to be one of the key factors promoting LUAD progression. We found the inferred cancer trajectory allows to identify LUAD susceptibility genetic variations using genome-wide association analysis. This result shows the opportunity for combining analysis of candidate genetic factors with disease progression. Furthermore, the trajectory showed clear evident mutation accumulation and clonal expansion along with the LUAD progression. Understanding how tumors evolve and identifying mutated genes will help guide cancer management. We investigated the clonal architectures and identified distinct clones and subclones in different LUAD branches. Validation of the model in multiple independent data sets and correlation analysis with clinical results demonstrate that our method is effective and unbiased. Author summary: Lung adenocarcinoma (LUAD) is a deadly tumor that remains approximately 15% survival rate in 5 years. The diverse reasons that made cancer progression and metastasis include genetic mutation, gene expression alteration, and so on. Understanding this dynamic process and identifying pivotal molecular events driving tumor progression is essential for improving LUAD diagnosis and treatment. Studying time-series data can simulate cancer evolution and determine the temporal patterns of molecular alterations. Unfortunately, it is difficult to collect complete time-series data from individual patients due to various reasons. In this study, we developed a trajectory model based on gene expression of LUAD patients and identified distinct progression branches. Moreover, we found a key gene BUB1B that could lead to aberrant spindle assembly checkpoint (SAC) signaling, and this appeared to be one of the key factors promoting LUAD progression. The trajectory shows clear evident mutation accumulation and clonal expansion along with the LUAD progression. We also investigated the clonal architectures and identified distinct clones and subclones in different LUAD branches. This clones and subclones can be used to identify combination therapy for LUAD. Our results demonstrate the biological utility and clinical application prospects of this progression model. [ABSTRACT FROM AUTHOR]

Thermodynamic stability is a crucial fitness constraint in protein evolution and is a central factor in shaping the sequence landscapes of proteins. The correlation between stability and molecular fitness depends on the mechanism that relates the biophysical property with biological function. In the simplest case, stability and fitness are related by the amount of folded protein. However, when proteins are toxic in the unfolded state, the fitness function shifts, resulting in higher stability under mutation-selection balance. Likewise, a higher population size results in a similar change in protein stability, as it magnifies the effect of the selection pressure in evolutionary dynamics. This study investigates how such factors affect the evolution of protein stability, site-specific mutation rates, and residue-residue covariation. To simulate evolutionary trajectories with realistic modeling of protein energetics, we develop an all-atom simulator of protein evolution, RosettaEvolve. By evolving proteins under different fitness functions, we can study how the fitness function affects the distribution of proposed and accepted mutations, site-specific rates, and the prevalence of correlated amino acid substitutions. We demonstrate that fitness pressure affects the proposal distribution of mutational effects, that changes in stability can largely explain variations in site-specific substitution rates in evolutionary trajectories, and that increased fitness pressure results in a stronger covariation signal. Our results give mechanistic insight into the evolutionary consequences of variation in protein stability and provide a basis to rationalize the strong covariation signal observed in natural sequence alignments. Author summary: Modern-day proteins are the result of the process of evolution. The fate of random substitutions at the nucleotide level is dependent on the fitness of the new gene variant. One of the strongest fitness pressures shaping the sequences of protein is thermodynamic stability; proteins must typically be stable to carry out its function and misfolded proteins can be toxic. To understand the importance of thermodynamic stability in protein evolution and to what extent it can explain natural sequence variation we have developed a method for simulating protein evolution using a three-dimensional structure and structure-based stability calculations. In the simulations, the strength of selection can be varied, and complete phylogenetic trees of a protein family can be generated. Using these simulations, we demonstrate how mutation rates at individual sites in a protein are coupled to the overall stability of the protein, and how the spectrum of accepted mutations is shaped by stability, and how strong interactions between residues in a protein can result in sequence covariation. [ABSTRACT FROM AUTHOR]

The majority of gene expression studies focus on the search for genes whose mean expression is different between two or more populations of samples in the so-called "differential expression analysis" approach. However, a difference in variance in gene expression may also be biologically and physiologically relevant. In the classical statistical model used to analyze RNA-sequencing (RNA-seq) data, the dispersion, which defines the variance, is only considered as a parameter to be estimated prior to identifying a difference in mean expression between conditions of interest. Here, we propose to evaluate four recently published methods, which detect differences in both the mean and dispersion in RNA-seq data. We thoroughly investigated the performance of these methods on simulated datasets and characterized parameter settings to reliably detect genes with a differential expression dispersion. We applied these methods to The Cancer Genome Atlas datasets. Interestingly, among the genes with an increased expression dispersion in tumors and without a change in mean expression, we identified some key cellular functions, most of which were related to catabolism and were overrepresented in most of the analyzed cancers. In particular, our results highlight autophagy, whose role in cancerogenesis is context-dependent, illustrating the potential of the differential dispersion approach to gain new insights into biological processes and to discover new biomarkers. Author summary: Gene expression is the process by which genetic information is translated into functional molecules. Transcription is the first step of this process, consisting of synthesizing messenger RNAs. During recent decades, genome-wide transcriptional profiling technologies have made it possible to assess the expression levels of thousands of genes in parallel in a variety of biological contexts. In statistical analyses, the expression of a gene is estimated by counting sequencing reads over a set of samples and is defined by two dimensions: mean and variance. The overwhelming majority of gene expression studies focus on identifying genes whose mean expression significantly changes when comparing samples of different conditions of interest to gain knowledge of biological processes. In this classical approach, the variance is usually considered only as a noise parameter to be estimated before assessing the mean expression. However, finely estimating the variance of expression may be biologically relevant since a modification of this parameter may reflect a change in gene expression regulation. Here, we propose to evaluate the performance of statistical methods that identify such differentially variant genes. We highlighted the potential of this approach by analyzing cancer datasets, thus identifying key cellular functions in tumor progression. [ABSTRACT FROM AUTHOR]

Transcription factors read the genome, fundamentally connecting DNA sequence to gene expression across diverse cell types. Determining how, where, and when TFs bind chromatin will advance our understanding of gene regulatory networks and cellular behavior. The 2017 ENCODE-DREAM in vivo Transcription-Factor Binding Site (TFBS) Prediction Challenge highlighted the value of chromatin accessibility data to TFBS prediction, establishing state-of-the-art methods for TFBS prediction from DNase-seq. However, the more recent Assay-for-Transposase-Accessible-Chromatin (ATAC)-seq has surpassed DNase-seq as the most widely-used chromatin accessibility profiling method. Furthermore, ATAC-seq is the only such technique available at single-cell resolution from standard commercial platforms. While ATAC-seq datasets grow exponentially, suboptimal motif scanning is unfortunately the most common method for TFBS prediction from ATAC-seq. To enable community access to state-of-the-art TFBS prediction from ATAC-seq, we (1) curated an extensive benchmark dataset (127 TFs) for ATAC-seq model training and (2) built "maxATAC", a suite of user-friendly, deep neural network models for genome-wide TFBS prediction from ATAC-seq in any cell type. With models available for 127 human TFs, maxATAC is the largest collection of high-performance TFBS prediction models for ATAC-seq. maxATAC performance extends to primary cells and single-cell ATAC-seq, enabling improved TFBS prediction in vivo. We demonstrate maxATAC's capabilities by identifying TFBS associated with allele-dependent chromatin accessibility at atopic dermatitis genetic risk loci. Author summary: Proteins called transcription factors interpret the genome, reading both DNA sequence and chromatin state, to orchestrate gene expression across the diversity of human cell types. In any given cell type, most chromatin is "inaccessible", and only those parts of the genetic code needed or likely to be needed soon are "accessible" for transcription factor binding to affect gene expression and cellular behavior. Hundreds of transcription factors are expressed in a given cell type and context (e.g., age, disease), and knowledge of their context-specific DNA binding sites is key to uncovering how transcription factors regulate cellular behaviors in health or disease. However, experimentally profiling the >1,600 human transcription factors across all cell types and contexts is infeasible. We built a suite of computational models "maxATAC" to predict transcription factor binding from a measurement of accessible chromatin, ATAC-seq. Importantly, ATAC-seq is feasible even at single-cell resolution. Thus, this data type, in combination with maxATAC, can be used to infer transcription factor binding sites in directly-relevant cell types isolated from physiological and disease settings, enabling insights into disease mechanisms, including how genetic variants and cellular context impact transcription factor binding, gene expression patterns and disease risk. [ABSTRACT FROM AUTHOR]

Background: Infection with human immunodeficiency virus type 1 (HIV) typically results from transmission of a small and genetically uniform viral population. Following transmission, the virus population becomes more diverse because of recombination and acquired mutations through genetic drift and selection. Viral intrahost genetic diversity remains a major obstacle to the cure the HIV; however, the association between intrahost diversity and disease progression markers has not been investigated in large and diverse cohorts for which the majority of the genome has been deep-sequenced. Viral load (VL) is a key progression marker and understanding of its relationship to viral intrahost genetic diversity could help design future strategies for HIV monitoring and treatment. Methods: We analysed deep-sequenced viral genomes from 2,650 treatment-naive HIV-infected persons to measure the intrahost genetic diversity of 2,447 genomic codon positions as calculated by Shannon entropy. We tested for associations between VL and amino acid (AA) entropy accounting for sex, age, race, duration of infection, and HIV population structure. Results: We confirmed that the intrahost genetic diversity is highest in the env gene. Furthermore, we showed that mean Shannon entropy is significantly associated with VL, especially in infections of >24 months duration. We identified 16 significant associations between VL (p-value<2.0x10-5) and Shannon entropy at AA positions which in our association analysis explained 13% of the variance in VL. Finally, equivalent analysis based on variation in HIV consensus sequences explained only 2% of VL variance. Conclusions: Our results elucidate that viral intrahost genetic diversity is associated with VL and could be used as a better disease progression marker than HIV consensus sequence variants, especially in infections of longer duration. We emphasize that viral intrahost diversity should be considered when studying viral genomes and infection outcomes. Trial registration: Samples included in this study were derived from participants who consented in the clinical trial, START (NCT00867048) (23), run by the International Network for Strategic Initiatives in Global HIV Trials (INSIGHT). All the participant sites are listed here: http://www.insight-trials.org/start/my_phpscript/participating.php?by=site Author summary: Viral intrahost genetic diversity complicates the cure of HIV; nevertheless, there is a lack of large and diverse cohort studies based on near-full and deep-sequenced HIV genomes. Here, we analysed deep-sequenced viral genomes from 2,650 demographically diverse and treatment-naive HIV-infected persons to measure the intrahost genetic diversity of 2,447 genomic codon positions as calculated by Shannon entropy. First, we verified that intrahost genetic diversity is highest in the env gene. Then, we identified that mean Shannon entropy positively associates with viral load, mostly in infections of >24 months duration. Finally, we showed that intrahost diversity in 16 HIV genomic positions significantly associated with viral load and explained 13% of the variance in viral load, whereas equivalent analysis based on variation in HIV consensus sequences explained only 2% of VL variance. Overall, our study shows that higher viral intrahost genetic diversity positively associates with viral load and is a better disease progression marker than variation of HIV consensus sequences. Our findings further suggest that better understanding of pathogen genomics is required to better address infectious diseases. [ABSTRACT FROM AUTHOR]

To date, ancient genome analyses have been largely confined to the study of single nucleotide polymorphisms (SNPs). Copy number variants (CNVs) are a major contributor of disease and of evolutionary adaptation, but identifying CNVs in ancient shotgun-sequenced genomes is hampered by typical low genome coverage (<1×) and short fragments (<80 bps), precluding standard CNV detection software to be effectively applied to ancient genomes. Here we present CONGA, tailored for genotyping CNVs at low coverage. Simulations and down-sampling experiments suggest that CONGA can genotype deletions >1 kbps with F-scores >0.75 at ≥1×, and distinguish between heterozygous and homozygous states. We used CONGA to genotype 10,002 outgroup-ascertained deletions across a heterogenous set of 71 ancient human genomes spanning the last 50,000 years, produced using variable experimental protocols. A fraction of these (21/71) display divergent deletion profiles unrelated to their population origin, but attributable to technical factors such as coverage and read length. The majority of the sample (50/71), despite originating from nine different laboratories and having coverages ranging from 0.44×-26× (median 4×) and average read lengths 52-121 bps (median 69), exhibit coherent deletion frequencies. Across these 50 genomes, inter-individual genetic diversity measured using SNPs and CONGA-genotyped deletions are highly correlated. CONGA-genotyped deletions also display purifying selection signatures, as expected. CONGA thus paves the way for systematic CNV analyses in ancient genomes, despite the technical challenges posed by low and variable genome coverage. Author summary: In parallel with developments in genomic technologies over the last decades, ancient genomics opened a new era in understanding the evolutionary history of populations and species. However, the field still needs novel computational methods for accurate and effective use of ancient genome data, which is mostly low-coverage and more challenging to analyse than modern-day genomes. Single nucleotide polymorphisms (SNPs), to date, have yet been the main source of information analysed in ancient genome studies. This is despite copy number variants (CNVs) harboring at least as much information as SNPs, especially with respect to natural selection. Here we developed CONGA, an algorithm for genotyping deletions and duplications in low-coverage genomes. We assessed its accuracy using simulations (with ancient-like data), and also studied its performance among 71 real ancient human genomes from different laboratories. We found that the common practice of authors filtering their ancient genome data before publishing prevents the reliable identification of duplications. Meanwhile, large (>1,000 base-pair) deletions can be detected even at quite low coverage (e.g. 0.5×). Deletions called in ancient genomes reflect population history and also show signs of negative selection. [ABSTRACT FROM AUTHOR]

Understanding the benefits and costs of recombination under different scenarios of evolutionary adaptation remains an open problem for theoretical and experimental research. In this study, we focus on finite populations evolving on neutral networks comprising viable and unfit genotypes. We provide a comprehensive overview of the effects of recombination by jointly considering different measures of evolvability and mutational robustness over a broad parameter range, such that many evolutionary regimes are covered. We find that several of these measures vary non-monotonically with the rates of mutation and recombination. Moreover, the presence of unfit genotypes that introduce inhomogeneities in the network of viable states qualitatively alters the effects of recombination. We conclude that conflicting trends induced by recombination can be explained by an emerging trade-off between evolvability on the one hand, and mutational robustness on the other. Finally, we discuss how different implementations of the recombination scheme in theoretical models can affect the observed dependence on recombination rate through a coupling between recombination and genetic drift. Author summary: Many genetic mechanisms have been invoked to explain the advantage of sex, but a coherent picture is still to emerge. Here we present a systematic theoretical and computational investigation of the effects of recombination in populations evolving on neutral fitness landscapes with unfit genotypes. We focus on populations that are large enough to be polymorphic, but nevertheless strongly affected by drift, which causes them to diffuse across the neutral network of viable genotypes. We identify a novel trade-off between evolvability, robustness and fitness that can lead to a dramatic reduction of the genetic diversity at large recombination rates. This disproves the common notion (often referred to as Weismann's hypothesis) that recombination generally increases diversity and evolvability, and instead highlights the interplay of recombination and mutational robustness. [ABSTRACT FROM AUTHOR]

Genetic variations affect behavior and cause disease but understanding how these variants drive complex traits is still an open question. A common approach is to link the genetic variants to intermediate molecular phenotypes such as the transcriptome using RNA-sequencing (RNA-seq). Paradoxically, these variants between the samples are usually ignored at the beginning of RNA-seq analyses of many model organisms. This can skew the transcriptome estimates that are used later for downstream analyses, such as expression quantitative trait locus (eQTL) detection. Here, we assessed the impact of reference-based analysis on the transcriptome and eQTLs in a widely-used mouse genetic population: the BXD panel of recombinant inbred lines. We highlight existing reference bias in the transcriptome data analysis and propose practical solutions which combine available genetic variants, genotypes, and genome reference sequence. The use of custom BXD line references improved downstream analysis compared to classical genome reference. These insights would likely benefit genetic studies with a transcriptomic component and demonstrate that genome references need to be reassessed and improved. Author summary: To understand how genetic variations affect behavior and cause disease it is common to quantify expression of transcripts by sequencing. Transcripts are extracted, fragmented, and the sequence of the fragments read. An important step for their quantification is to virtually assign the different fragments to the transcript they originate from using a reference genome. Reference genomes are costly to build, so usually only one high-quality reference per animal model species is available. When comparing genetically different individuals, using a single reference may introduce a bias because it might be more similar to some individuals than to others. Paradoxically, the variations at the core of genetic studies are thus ignored at the start of the analysis. We built customized references with known genetic variants for each of the mouse lines we had and quantified the impact of the reference at different levels of the bioinformatic analysis. We found that using customized references reduced the bias compared to using a single reference. Our study uses publicly available data and tools, so others can easily implement this improvement in their analyses. [ABSTRACT FROM AUTHOR]

Proximal genetic variants are frequently correlated, implying that the corresponding effect sizes detected by genome-wide association studies (GWAS) are also not independent. Methods already exist to account for this when aggregating effects from a single GWAS across genes or pathways. Here we present a rigorous yet fast method for detecting genes with coherent association signals for two traits, facilitating cross-GWAS analyses. To this end, we devised a new significance test for the covariance of datapoints not drawn independently but with a known inter-sample covariance structure. We show that the distribution of its test statistic is a linear combination of χ2 distributions with positive and negative coefficients. The corresponding cumulative distribution function can be efficiently calculated with Davies' algorithm at high precision. We apply this general framework to test for dependence between SNP-wise effect sizes of two GWAS at the gene level. We extend this test to detect also gene-wise causal links. We demonstrate the utility of our method by uncovering potential shared genetic links between the severity of COVID-19 and (1) being prescribed class M05B medication (drugs affecting bone structure and mineralization), (2) rheumatoid arthritis, (3) vitamin D (25OHD), and (4) serum calcium concentrations. Our method detects a potential role played by chemokine receptor genes linked to TH1 versus TH2 immune response, a gene related to integrin beta-1 cell surface expression, and other genes potentially impacting the severity of COVID-19. Our approach will be useful for similar analyses involving datapoints with known auto-correlation structures. Author summary: Genome-wide association studies (GWAS) deliver effect size estimates of a given trait for millions of Single Nucleotide Polymorphisms (SNPs). Powerful tools already exist using these summary statistics to elucidate the global joint genetic contribution to a pair of traits, such as cross-trait LD-score regression, but these methods cannot reveal the joint contributions at the level of genes and pathways. Here we present a novel methodology to co-analyze the association data from a pair of GWAS to identify genes and pathways that may be relevant to both individual traits. Our test novelty is that a gene is considered co-relevant if the SNP-wise effects from both GWAS tend to have the same sign and magnitude in the gene window. This is different from the commonly used approach asking only for the aggregate signals from two GWAS to be jointly significant. Our method is feasible due to novel insight into the product-normal distribution. We apply our new method to test for co-significant genes for severe COVID-19 and conditions leading to the prescription of common medications. Of the 23 medication classes we tested for coherent co-significant genes, only one, M05B (drugs affecting bone structure and mineralization), yielded Bonferroni significant hits. We then searched for available GWAS data for related conditions and found that also rheumatoid arthritis, calcium concentration, and vitamin D are traits pointing to several co-relevant genes in our new coherence analysis. Furthermore, testing for anti-coherence showed that the medication classes H03A (thyroid preparations), R03A, and R03BA (drugs for obstructive airway diseases) feature Bonferroni co-significant genes. Our joint analysis provides new insights into potential COVID-19 disease mechanisms. [ABSTRACT FROM AUTHOR]

Unraveling the complex demographic histories of natural populations is a central problem in population genetics. Understanding past demographic events is of general anthropological interest, but is also an important step in establishing accurate null models when identifying adaptive or disease-associated genetic variation. An important class of tools for inferring past population size changes from genomic sequence data are Coalescent Hidden Markov Models (CHMMs). These models make efficient use of the linkage information in population genomic datasets by using the local genealogies relating sampled individuals as latent states that evolve along the chromosome in an HMM framework. Extending these models to large sample sizes is challenging, since the number of possible latent states increases rapidly. Here, we present our method CHIMP (CHMM History-Inference Maximum-Likelihood Procedure), a novel CHMM method for inferring the size history of a population. It can be applied to large samples (hundreds of haplotypes) and only requires unphased genomes as input. The two implementations of CHIMP that we present here use either the height of the genealogical tree (TMRCA) or the total branch length, respectively, as the latent variable at each position in the genome. The requisite transition and emission probabilities are obtained by numerically solving certain systems of differential equations derived from the ancestral process with recombination. The parameters of the population size history are subsequently inferred using an Expectation-Maximization algorithm. In addition, we implement a composite likelihood scheme to allow the method to scale to large sample sizes. We demonstrate the efficiency and accuracy of our method in a variety of benchmark tests using simulated data and present comparisons to other state-of-the-art methods. Specifically, our implementation using TMRCA as the latent variable shows comparable performance and provides accurate estimates of effective population sizes in intermediate and ancient times. Our method is agnostic to the phasing of the data, which makes it a promising alternative in scenarios where high quality data is not available, and has potential applications for pseudo-haploid data. Author summary: The demograpic history of natural populations shapes their genetic variation. The genomes of contemporary individuals can thus be used to unravel past migration events and population size changes, which is of anthropological interest. Moreover, it is also important to uncover these past events for studies investigating disease related genetic variation, since past demographic events can confound such analyses. Here we present a novel method for inferring the size history of a given population from full-genome sequencing data of contemporary individuals. Our method is based on a Coalescent Hidden Markov model framework, a model frequently applied to this type of inference. A key component of the model is the representation of unobserved local genealogical relationships among the sampled individuals as latent states. This is achieved by numerically solving certain differential equations that describe the distributions of these quantities and ultimately enables inference of past population size changes. Other methods performing similar inference rely on availability of high quality genomic data, whereas we demonstrate that our method can be applied in situations with limited data quality. [ABSTRACT FROM AUTHOR]

Movement of individuals between populations or demes is often restricted, especially between geographically isolated populations. The structured coalescent provides an elegant theoretical framework for describing how movement between populations shapes the genealogical history of sampled individuals and thereby structures genetic variation within and between populations. However, in the presence of recombination an individual may inherit different regions of their genome from different parents, resulting in a mosaic of genealogical histories across the genome, which can be represented by an Ancestral Recombination Graph (ARG). In this case, different genomic regions may have different ancestral histories and so different histories of movement between populations. Recombination therefore poses an additional challenge to phylogeographic methods that aim to reconstruct the movement of individuals from genealogies, although also a potential benefit in that different loci may contain additional information about movement. Here, we introduce the Structured Coalescent with Ancestral Recombination (SCAR) model, which builds on recent approximations to the structured coalescent by incorporating recombination into the ancestry of sampled individuals. The SCAR model allows us to infer how the migration history of sampled individuals varies across the genome from ARGs, and improves estimation of key population genetic parameters such as population sizes, recombination rates and migration rates. Using the SCAR model, we explore the potential and limitations of phylogeographic inference using full ARGs. We then apply the SCAR to lineages of the recombining fungus Aspergillus flavus sampled across the United States to explore patterns of recombination and migration across the genome. Author summary: Phylogeographic methods are widely used to reconstruct the historical movement of individuals between different populations. When applied to infectious pathogens, these methods are often used to reconstruct the origin or source of novel pathogen lineages. Most existing phylogeographic methods reconstruct movement based on a single phylogenetic tree, which is assumed to reflect the genetic ancestry of all sampled individuals. However in populations undergoing recombination, genetic material can be exchanged between lineages such that individuals may inherit different regions of their genome from different ancestors. In this case, phylogenetic relationships among individuals can only be captured by a reticulated network rather than any single tree. Ancestral Recombination Graphs (ARGs) provide one way of capturing these reticulate relationships and we develop new models that allow for demographic inference of historical population sizes, recombination rates and migration rates between subpopulations from ARGs. By accounting for recombination, our models not only allow for accurate demographic inference, but can take full advantage of the additional information contained in ARGs about how ancestry varies across genomes to more precisely reconstruct the movement of genetic material between populations. [ABSTRACT FROM AUTHOR]

When two influenza viruses co-infect the same cell, they can exchange genome segments in a process known as reassortment. Reassortment is an important source of genetic diversity and is known to have been involved in the emergence of most pandemic influenza strains. However, because of the difficulty in identifying reassortments events from viral sequence data, little is known about its role in the evolution of the seasonal influenza viruses. Here we introduce TreeKnit, a method that infers ancestral reassortment graphs (ARG) from two segment trees. It is based on topological differences between trees, and proceeds in a greedy fashion by finding regions that are compatible in the two trees. Using simulated genealogies with reassortments, we show that TreeKnit performs well in a wide range of settings and that it is as accurate as a more principled bayesian method, while being orders of magnitude faster. Finally, we show that it is possible to use the inferred ARG to better resolve segment trees and to construct more informative visualizations of reassortments. Author summary: Influenza viruses evolve quickly and escape immune defenses which requires frequent update of vaccines. Understanding this evolution is key to an effective public health response. The genome of influenza viruses is made up of 8 pieces called segments, each coding for different viral proteins. Within each segment, evolution is an asexual process in which genetic diversity is generated by mutations. But influenza also diversifies through reassortment which can occur when two different viruses infect the same cell: offsprings can then contain a combination of segments from both viruses. Reassortment is akin to sexual reproduction and can generate viruses that combine segments from diverged viral lineages. Reassortment is a crucial component of viral evolution, but it is challenging to reconstruct where reassortments happened and which segments share history. Here, we develop a method called TreeKnit to detect reassortment events. TreeKnit is based on genealogical trees of single segments that can be reconstructed using standard bioinformatics tools. Inconsistencies between these trees are then used as signs of reassortment. We show that TreeKnit is as accurate as other recent methods, but runs much faster. Our method will facilitate the study of reassortment and its consequences for influenza evolution. [ABSTRACT FROM AUTHOR]

Building an accurate disease risk prediction model is an essential step in the modern quest for precision medicine. While high-dimensional genomic data provides valuable data resources for the investigations of disease risk, their huge amount of noise and complex relationships between predictors and outcomes have brought tremendous analytical challenges. Deep learning model is the state-of-the-art methods for many prediction tasks, and it is a promising framework for the analysis of genomic data. However, deep learning models generally suffer from the curse of dimensionality and the lack of biological interpretability, both of which have greatly limited their applications. In this work, we have developed a deep neural network (DNN) based prediction modeling framework. We first proposed a group-wise feature importance score for feature selection, where genes harboring genetic variants with both linear and non-linear effects are efficiently detected. We then designed an explainable transfer-learning based DNN method, which can directly incorporate information from feature selection and accurately capture complex predictive effects. The proposed DNN-framework is biologically interpretable, as it is built based on the selected predictive genes. It is also computationally efficient and can be applied to genome-wide data. Through extensive simulations and real data analyses, we have demonstrated that our proposed method can not only efficiently detect predictive features, but also accurately predict disease risk, as compared to many existing methods. Author summary: Accurate disease risk prediction is an essential step towards precision medicine. Deep learning models have achieved the state-of-the-art performance for many prediction tasks. However, they generally suffer from the curse of dimensionality and lack of biological interpretability, both of which have greatly limited their applications to the prediction analysis of whole-genome sequencing data. We present here an explainable deep transfer learning model for the analysis of high-dimensional genomic data. Our proposed method can detect predictive genes that harbor genetic variants with both linear and non-linear effects via the proposed group-wise feature importance score. It can also efficiently and accurately model disease risk based on the detected predictive genes using the proposed transfer-learning based network architecture. Our proposed method is built at the gene level, and thus is much more biologically interpretable. It is also computationally efficiently and can be applied to whole-exome sequencing data that have millions of potential predictors. Through both simulation studies and the analysis of whole-exome data obtained from the Alzheimer's Disease Neuroimaging Initiative, we have demonstrated that our method can efficiently detect predictive genes and it has better prediction performance than many existing methods. [ABSTRACT FROM AUTHOR]

Gene-based association analysis is an effective gene-mapping tool. Many gene-based methods have been proposed recently. However, their power depends on the underlying genetic architecture, which is rarely known in complex traits, and so it is likely that a combination of such methods could serve as a universal approach. Several frameworks combining different gene-based methods have been developed. However, they all imply a fixed set of methods, weights and functional annotations. Moreover, most of them use individual phenotypes and genotypes as input data. Here, we introduce sumSTAAR, a framework for gene-based association analysis using summary statistics obtained from genome-wide association studies (GWAS). It is an extended and modified version of STAAR framework proposed by Li and colleagues in 2020. The sumSTAAR framework offers a wider range of gene-based methods to combine. It allows the user to arbitrarily define a set of these methods, weighting functions and probabilities of genetic variants being causal. The methods used in the framework were adapted to analyse genes with large number of SNPs to decrease the running time. The framework includes the polygene pruning procedure to guard against the influence of the strong GWAS signals outside the gene. We also present new improved matrices of correlations between the genotypes of variants within genes. These matrices estimated on a sample of 265,000 individuals are a state-of-the-art replacement of widely used matrices based on the 1000 Genomes Project data. Author summary: Gene-based association analysis is an effective gene mapping tool. Quite a few frameworks have been proposed recently for gene-based association analysis using a combination of different methods. However, all of these frameworks have at least one of the disadvantages: they use a fixed set of methods, they cannot use functional annotations, or they use individual phenotypes and genotypes as input data. To overcome these limitations, we propose sumSTAAR, a framework for gene-based association analysis using GWAS summary statistics. Our framework allows the user to arbitrarily define a set of the methods and functional annotations. Moreover, we adopted the methods for the analysis of genes with a large number of SNPs to decrease the running time. The framework includes the polygene pruning procedure to guard against the influence of the strong GWAS signals outside the gene. We also present new improved matrices of correlations between the genotypes of variants within genes, which now allows to include ultra-rare variants (MAF < 10−4) in analysis. [ABSTRACT FROM AUTHOR]

Genomewide association studies (GWAS) have identified a large number of loci associated with neuropsychiatric traits, however, understanding the molecular mechanisms underlying these loci remains difficult. To help prioritize causal variants and interpret their functions, computational methods have been developed to predict regulatory effects of non-coding variants. An emerging approach to variant annotation is deep learning models that predict regulatory functions from DNA sequences alone. While such models have been trained on large publicly available dataset such as ENCODE, neuropsychiatric trait-related cell types are under-represented in these datasets, thus there is an urgent need of better tools and resources to annotate variant functions in such cellular contexts. To fill this gap, we collected a large collection of neurodevelopment-related cell/tissue types, and trained deep Convolutional Neural Networks (ResNet) using such data. Furthermore, our model, called MetaChrom, borrows information from public epigenomic consortium to improve the accuracy via transfer learning. We show that MetaChrom is substantially better in predicting experimentally determined chromatin accessibility variants than popular variant annotation tools such as CADD and delta-SVM. By combining GWAS data with MetaChrom predictions, we prioritized 31 SNPs for Schizophrenia, suggesting potential risk genes and the biological contexts where they act. In summary, MetaChrom provides functional annotations of any DNA variants in the neuro-development context and the general method of MetaChrom can also be extended to other disease-related cell or tissue types. Author summary: A large number of genetic variants have been statistically associated with the risks of common diseases. However, whether such variants are actual risk variants and when and where they function are often unknown. To address this challenge, machine learning methods have been developed to predict functional variants in specific cellular contexts. These methods correlate DNA sequences with their biological functions, e.g. enhancer activities, and can predict effects of single base mutations. Nevertheless, the training data used by existing methods often lack neurodevelopment-related cell types, thus annotating variant effects in neuropsychiatric genetics remains difficult. In this work, we fill this gap by collecting a large set of regulatory genomic datasets from fetal and adult brain, from iPSC-based cellular models and brain organoids. We trained deep learning models on this data, and further improved its performance by borrowing information from large external datasets, a strategy known as transfer learning. Our tool, MetaChrom, is substantially better at predicting experimentally determined regulatory variants than current methods, and helps us identify candidate risk variants of Schizophrenia. We believe MetaChrom provides a valuable tool for the neuropsychiatric genetic community, and the software can be of interest to researchers in other fields as well. [ABSTRACT FROM AUTHOR]

Genome wide association studies (GWASs) for complex traits have implicated thousands of genetic loci. Most GWAS-nominated variants lie in noncoding regions, complicating the systematic translation of these findings into functional understanding. Here, we leverage convolutional neural networks to assist in this challenge. Our computational framework, peaBrain, models the transcriptional machinery of a tissue as a two-stage process: first, predicting the mean tissue specific abundance of all genes and second, incorporating the transcriptomic consequences of genotype variation to predict individual abundance on a subject-by-subject basis. We demonstrate that peaBrain accounts for the majority (>50%) of variance observed in mean transcript abundance across most tissues and outperforms regularized linear models in predicting the consequences of individual genotype variation. We highlight the validity of the peaBrain model by calculating non-coding impact scores that correlate with nucleotide evolutionary constraint that are also predictive of disease-associated variation and allele-specific transcription factor binding. We further show how these tissue-specific peaBrain scores can be leveraged to pinpoint functional tissues underlying complex traits, outperforming methods that depend on colocalization of eQTL and GWAS signals. We subsequently: (a) derive continuous dense embeddings of genes for downstream applications; (b) highlight the utility of the model in predicting transcriptomic impact of small molecules and shRNA (on par with in vitro experimental replication of external test sets); (c) explore how peaBrain can be used to model difficult-to-study processes (such as neural induction); and (d) identify putatively functional eQTLs that are missed by high-throughput experimental approaches. Author summary: High-throughput assays are the cornerstone of modern drug discovery and a useful tool to translating the hundreds of genetic discoveries associated with human traits and disease into functional understanding. All high-throughput assays can be described as empirical assessments of the activity of biological entities (e.g., genetic variation, DNA sequences, small molecules) by a standardized output, usually in the form of optically detectable labels (i.e., reporters), or more rarely, using (scalable) high-dimensional measurements (e.g., L1000, RNA-seq). Here, we introduce a modular and readily-extensible computational framework, called peaBrain, that leverages convolutional neural network architecture to enable in silico recapitulation of certain features of these high-throughput assays. We show that peaBrain can predict the expression of genes in a tissue-specific manner and outperforms regularized linear models in predicting the consequences of individual genotype variation. We further highlight the utility of the framework in predicting transcriptomic impact of small molecules and shRNA (on par with in vitro experimental replication of external test sets), explore how peaBrain can be used to model difficult-to-study processes (such as neural induction), and finally, identify putatively functional eQTLs that are missed by high-throughput experimental approaches. [ABSTRACT FROM AUTHOR]

To date, genome-wide association studies have identified thousands of statistically-significant associations between genetic variants, and phenotypes related to a myriad of traits and diseases. A key goal for human-genetics research is to translate these associations into functional mechanisms. Popular gene-set analysis tools, like MAGMA, map variants to genes they might affect, and then integrate genome-wide association study data (that is, variant-level associations for a phenotype) to score genes for association with a phenotype. Gene scores are subsequently used in competitive gene-set analyses to identify biological processes that are enriched for phenotype association. By default, variants are mapped to genes in their proximity. However, many variants that affect phenotypes are thought to act at regulatory elements, which can be hundreds of kilobases away from their target genes. Thus, we explored the idea of augmenting a proximity-based mapping scheme with publicly-available datasets of regulatory interactions. We used MAGMA to analyze genome-wide association study data for ten different phenotypes, and evaluated the effects of augmentation by comparing numbers, and identities, of genes and gene sets detected as statistically significant between mappings. We detected several pitfalls and confounders of such "augmented analyses", and introduced ways to control for them. Using these controls, we demonstrated that augmentation with datasets of regulatory interactions only occasionally strengthened the enrichment for phenotype association amongst (biologically-relevant) gene sets for different phenotypes. Still, in such cases, genes and regulatory elements responsible for the improvement could be pinpointed. For instance, using brain regulatory-interactions for augmentation, we were able to implicate two acetylcholine receptor subunits involved in post-synaptic chemical transmission, namely CHRNB2 and CHRNE, in schizophrenia. Collectively, our study presents a critical approach for integrating regulatory interactions into gene-set analyses for genome-wide association study data, by introducing various controls to distinguish genuine results from spurious discoveries. Author summary: Genome-wide association studies have identified thousands of statistically-significant associations between genetic variants and diverse phenotypes, but translating these associations into functional mechanisms remains a major challenge. Gene-set analyses for genome-wide association study data, which consider the signal contained in practically all variants tested for association with a phenotype, seek functionally-annotated collections of genes that are enriched for phenotype association. To do so, variants must first be mapped to genes they are thought to affect. Usually, variants are mapped to genes in their proximity, even if it is believed that the effects of many variants are, in fact, mediated via long-range regulatory interactions. Here, we demonstrate that augmenting a proximity-based mapping scheme with datasets of regulatory interactions, can strengthen the enrichment for phenotype association amongst gene sets that are biologically-relevant to phenotypes. Yet, by implementing controls, we reveal that apparent improvements are often trivial, or can be explained by confounding. We discuss some examples of meaningful improvements, and pinpoint genes and regulatory elements responsible for them. Incorporating regulatory interactions into gene-set analyses for genome-wide association study data is expected to become the norm, and our work outlines a critical approach for doing so. [ABSTRACT FROM AUTHOR]

Identifying copy number variants (CNVs) can provide diagnoses to patients and provide important biological insights into human health and disease. Current exome and targeted sequencing approaches cannot detect clinically and biologically-relevant CNVs outside their target area. We present SavvyCNV, a tool which uses off-target read data from exome and targeted sequencing data to call germline CNVs genome-wide. Up to 70% of sequencing reads from exome and targeted sequencing fall outside the targeted regions. We have developed a new tool, SavvyCNV, to exploit this 'free data' to call CNVs across the genome. We benchmarked SavvyCNV against five state-of-the-art CNV callers using truth sets generated from genome sequencing data and Multiplex Ligation-dependent Probe Amplification assays. SavvyCNV called CNVs with high precision and recall, outperforming the five other tools at calling CNVs genome-wide, using off-target or on-target reads from targeted panel and exome sequencing. We then applied SavvyCNV to clinical samples sequenced using a targeted panel and were able to call previously undetected clinically-relevant CNVs, highlighting the utility of this tool within the diagnostic setting. SavvyCNV outperforms existing tools for calling CNVs from off-target reads. It can call CNVs genome-wide from targeted panel and exome data, increasing the utility and diagnostic yield of these tests. SavvyCNV is freely available at https://github.com/rdemolgen/SavvySuite. Author summary: We have created SavvyCNV, a new tool for calling genetic variants. Large regions of the genome can be deleted or duplicated–these variants can have important consequences, for example causing a patient's genetic disease. However, many standard genetic tests only target a small fraction of the genome and will miss variants outside of these regions. Therefore, we developed a tool to exploit sequencing data which falls outside of these regions (due to flaws in the targeting process) to call large deletions and duplications. This allows large deletions and duplications to be detected anywhere in the genome. Researchers and diagnostic laboratories can use this tool to discover more genetic variants by re-analysing their sequencing data. [ABSTRACT FROM AUTHOR]

All proteomes contain both proteins and polypeptide segments that don't form a defined three-dimensional structure yet are biologically active—called intrinsically disordered proteins and regions (IDPs and IDRs). Most of these IDPs/IDRs lack useful functional annotation limiting our understanding of their importance for organism fitness. Here we characterized IDRs using protein sequence annotations of functional sites and regions available in the UniProt knowledgebase ("UniProt features": active site, ligand-binding pocket, regions mediating protein-protein interactions, etc.). By measuring the statistical enrichment of twenty-five UniProt features in 981 IDRs of 561 human proteins, we identified eight features that are commonly located in IDRs. We then collected the genetic variant data from the general population and patient-based databases and evaluated the prevalence of population and pathogenic variations in IDPs/IDRs. We observed that some IDRs tolerate 2 to 12-times more single amino acid-substituting missense mutations than synonymous changes in the general population. However, we also found that 37% of all germline pathogenic mutations are located in disordered regions of 96 proteins. Based on the observed-to-expected frequency of mutations, we categorized 34 IDRs in 20 proteins (DDX3X, KIT, RB1, etc.) as intolerant to mutation. Finally, using statistical analysis and a machine learning approach, we demonstrate that mutation-intolerant IDRs carry a distinct signature of functional features. Our study presents a novel approach to assign functional importance to IDRs by leveraging the wealth of available genetic data, which will aid in a deeper understating of the role of IDRs in biological processes and disease mechanisms. Author summary: Intrinsically disordered regions (IDRs) in proteins are typically not considered to be functionally as important as the structured parts. However, it is becoming evident that both structured and disordered regions are essential for the repertoire of protein functions. Nevertheless, most of these largely flexible and functionally dynamic protein regions remain uncharacterized. Nevertheless, most of these largely flexible and functionally dynamic protein regions remain uncharacterized. Here, informed by human genetic diversity (i.e., genetic variations from the general population and patients), we identified the IDRs that are more frequently mutated in patients than in relatively healthy individuals, and further show that they carry a set of characteristic functional features. This approach provides a different and effective means to identify unannotated disordered protein segments that are biologically important and lead to pathogenesis upon mutation. [ABSTRACT FROM AUTHOR]

A critical step in studying biological features (e.g., genetic variants, gene families, metabolic capabilities, or taxa) is assessing their diversity and distribution among a sample of individuals. Accurate assessments of these patterns are essential for linking features to traits or outcomes of interest and understanding their functional impact. Consequently, it is of crucial importance that the measures employed for quantifying feature diversity can perform robustly under any evolutionary scenario. However, the standard measures used for quantifying and comparing the distribution of features, such as prevalence, phylogenetic diversity, and related approaches, either do not take into consideration evolutionary history, or assume strictly vertical patterns of inheritance. Consequently, these approaches cannot accurately assess diversity for features that have undergone recombination or horizontal transfer. To address this issue, we have devised RecPD, a novel recombination-aware phylogenetic-diversity statistic for measuring the distribution and diversity of features under all evolutionary scenarios. RecPD utilizes ancestral-state reconstruction to map the presence / absence of features onto ancestral nodes in a species tree, and then identifies potential recombination events in the evolutionary history of the feature. We also derive several related measures from RecPD that can be used to assess and quantify evolutionary dynamics and correlation of feature evolutionary histories. We used simulation studies to show that RecPD reliably reconstructs feature evolutionary histories under diverse recombination and loss scenarios. We then applied RecPD in two diverse real-world scenarios including a preliminary study type III effector protein families secreted by the plant pathogenic bacterium Pseudomonas syringae and growth phenotypes of the Pseudomonas genus and demonstrate that prevalence is an inadequate measure that obscures the potential impact of recombination. We believe RecPD will have broad utility for revealing and quantifying complex evolutionary processes for features at any biological level. Author summary: Phylogenetic diversity is an important concept utilized in evolutionary ecology which has extensive applications in population genetics to help us understand how evolutionary processes have distributed genetic variation among individuals of a species, and how this impacts phenotypic diversification over time. However, existing approaches for studying phylogenetic diversity largely assume that the genetic features follow vertical inheritance, which is frequently violated in the case of microbial genomes due to horizontal transfer. To address this shortcoming, we present RecPD, a recombination-aware phylogenetic diversity measure, which incorporates ancestral state reconstruction to quantify the phylogenetic diversity of genetic features mapped onto a species phylogeny. Through simulation experiments we show that RecPD robustly reconstructs the evolutionary histories of features evolving under various scenarios of recombination and loss. When applied to a real-world example of type III secreted effector protein families from the plant pathogenic bacterium Pseudomonas syringae, RecPD reveals that horizontal transfer has played an important role in shaping the phylogenetic distributions of a substantial proportion of families across the P. syringae species complex. Furthermore, we demonstrate that the traditional measures of feature prevalence are unsuitable as a measure for comparing feature diversity. We also provide a R package implementation of RecPD for public use: https://github.com/cedatorma/recpd. [ABSTRACT FROM AUTHOR]

Insertions and deletions (indels) in human genomes are associated with a wide range of phenotypes, including various clinical disorders. High-throughput, next generation sequencing (NGS) technologies enable the detection of short genetic variants, such as single nucleotide variants (SNVs) and indels. However, the variant calling accuracy for indels remains considerably lower than for SNVs. Here we present a comparative study of the performance of variant calling tools for indel calling, evaluated with a wide repertoire of NGS datasets. While there is no single optimal tool to suit all circumstances, our results demonstrate that the choice of variant calling tool greatly impacts the precision and recall of indel calling. Furthermore, to reliably detect indels, it is essential to choose NGS technologies that offer a long read length and high coverage coupled with specific variant calling tools. Author summary: The development of next generation sequencing (NGS) technologies and computational algorithms enabled the large scale, simultaneous detection of a wide range of genetic variants, such as single nucleotide variants as well as insertions and deletions (indels), which may confer potential clinical significance. Recently, many studies have been conducted to evaluate variant calling tools for indel calling. However, the optimal indel size range for different variant calling tools remains unclear. A good benchmarking dataset for indel calling evaluation should contain biologically representative high-confident indels with a wide size range and preferably come from various sequencing settings. In this article, we created a semi-simulated whole genome sequencing dataset where the sequencing data were computationally generated. The indels in the semi-simulated genome were incorporated from a real human sample to represent biologically realistic indels and to avoid the inclusion of variants due to potential technical sequencing errors. Furthermore, we used three real-world NGS datasets generated by whole genome or targeted sequencing to further evaluate our candidate tools. Our results demonstrated that variant calling tools vary greatly in calling different sizes of indels. Deletion calling and insertion calling also showed differences among the tools. The sequencing settings in coverage and read length also had a great impact on indel calling. Our results suggest that the accuracy of indel calling was dependent on the combination of a variant calling tool, indel size range, and sequencing settings. [ABSTRACT FROM AUTHOR]

The massive assessment of immune evasion due to viral mutations that increase COVID-19 susceptibility can be computationally facilitated. The adaptive cytotoxic T response is critical during primary infection and the generation of long-term protection. Here, potential HLA class I epitopes in the SARS-CoV-2 proteome were predicted for 2,915 human alleles of 71 families using the netMHCIpan EL algorithm. Allele families showed extreme epitopic differences, underscoring genetic variability of protective capacity between humans. Up to 1,222 epitopes were associated with any of the twelve supertypes, that is, allele clusters covering 90% population. Next, from all mutations identified in ~118,000 viral NCBI isolates, those causing significant epitope score reduction were considered epitope escape mutations. These mutations mainly involved non-conservative substitutions at the second and C-terminal position of the ligand core, or total ligand removal by large recurrent deletions. Escape mutations affected 47% of supertype epitopes, which in 21% of cases concerned isolates from two or more sub-continental areas. Some of these changes were coupled, but never surpassed 15% of evaded epitopes for the same supertype in the same isolate, except for B27. In contrast to most supertypes, eight allele families mostly contained alleles with few SARS-CoV-2 ligands. Isolates harboring cytotoxic escape mutations for these families co-existed geographically within sub-Saharan and Asian populations enriched in these alleles according to the Allele Frequency Net Database. Collectively, our findings indicate that escape mutation events have already occurred for half of HLA class I supertype epitopes. However, it is presently unlikely that, overall, it poses a threat to the global population. In contrast, single and double mutations for susceptible alleles may be associated with viral selective pressure and alarming local outbreaks. The integration of genomic, geographical and immunoinformatic information eases the surveillance of variants potentially affecting the global population, as well as minority subpopulations. Author summary: The cytotoxic T cell response, a type of immune response dependent upon an individual's genetics that does not require antibodies, is critical for neutralizing SARS-CoV-2 infection. The potential bypass of the cytotoxic T cell response by mutations acquired by the virus after one year of the pandemic is therefore of maximal concern. We have approached the complexity of human variability and more than 100,000 viral genomes in this respect using a computational strategy. We have detected numerous mutations in these genomes that mask some viral regions involved in the cytotoxic response. However, the accumulation of these changes in independent isolates is still too low to threaten the global human population. In contrast, our protocol has identified mutations that may be relevant for specific populations and minorities with cytotoxic genetic backgrounds susceptible to SARS-CoV-2 infection. Some viral variants co-existed in the same country with these human communities, which warrants deeper surveillance in these cases to prevent local outbreaks. Our study supports the integration of massive data of different natures in the surveillance of viral pandemics. [ABSTRACT FROM AUTHOR]

The protein structure field is experiencing a revolution. From the increased throughput of techniques to determine experimental structures, to developments such as cryo-EM that allow us to find the structures of large protein complexes or, more recently, the development of artificial intelligence tools, such as AlphaFold, that can predict with high accuracy the folding of proteins for which the availability of homology templates is limited. Here we quantify the effect of the recently released AlphaFold database of protein structural models in our knowledge on human proteins. Our results indicate that our current baseline for structural coverage of 48%, considering experimentally-derived or template-based homology models, elevates up to 76% when including AlphaFold predictions. At the same time the fraction of dark proteome is reduced from 26% to just 10% when AlphaFold models are considered. Furthermore, although the coverage of disease-associated genes and mutations was near complete before AlphaFold release (69% of Clinvar pathogenic mutations and 88% of oncogenic mutations), AlphaFold models still provide an additional coverage of 3% to 13% of these critically important sets of biomedical genes and mutations. Finally, we show how the contribution of AlphaFold models to the structural coverage of non-human organisms, including important pathogenic bacteria, is significantly larger than that of the human proteome. Overall, our results show that the sequence-structure gap of human proteins has almost disappeared, an outstanding success of direct consequences for the knowledge on the human genome and the derived medical applications. Author summary: Protein structures are key to understand many biological phenomena at the molecular scale: from the effects of genetic variation to how different proteins interact with each other to create molecular pathways that, together, have a biological function. Obtaining experimental structures, however, is extremely consuming in terms of both, time and resources. For this and other reasons, scientists have long worked to develop computational approaches that predict the structure of a protein using only its sequence as input. Recently, a group of scientists at Deepmind have developed AlphaFold2, a computational tool that is extremely accurate at this task. Moreover, they have used this tool to predict the structures of all human proteins. In this manuscript we provide an overview of the structural coverage of the human proteome before AlphaFold models were released and how much we have gained thanks to these models. We also show how the gain affects our understanding of human pathogenic variants, both germline and somatic. Finally, we provide evidence suggesting that the gain in non-human organisms is larger than for the human proteome, particularly in the case of bacteria. [ABSTRACT FROM AUTHOR]

Genome-wide association studies rely on the statistical inference of untyped variants, called imputation, to increase the coverage of genotyping arrays. However, the results are often suboptimal in populations underrepresented in existing reference panels and array designs, since the selected single nucleotide polymorphisms (SNPs) may fail to capture population-specific haplotype structures, hence the full extent of common genetic variation. Here, we propose to sequence the full genomes of a small subset of an underrepresented study cohort to inform the selection of population-specific add-on tag SNPs and to generate an internal population-specific imputation reference panel, such that the remaining array-genotyped cohort could be more accurately imputed. Using a Tanzania-based cohort as a proof-of-concept, we demonstrate the validity of our approach by showing improvements in imputation accuracy after the addition of our designed add-on tags to the base H3Africa array. Author summary: Genome-wide association studies, which study the association between genetic variants and various phenotypes, typically rely on genotyping arrays. Only a small proportion of genetic variants within the genome are typed on genotyping arrays. Untyped variants are statistically inferred through a process known as genotype imputation, where correlations between variants (haplotypes) observed in external reference panels are leveraged to infer untyped variants in the study population. However, for study populations that are underrepresented in existing reference panels, the quality of imputation is often sub-optimal. This is because typed variants incorporated on existing genotyping arrays can be unsuitable for the study population, and haplotype structures can be different between the reference and the study population. Here, we illustrate an approach to select a custom set of population-specific typed variants to improve genotype imputation in such underrepresented populations. [ABSTRACT FROM AUTHOR]

Expression QTL (eQTL) analyses have suggested many genes mediating genome-wide association study (GWAS) signals but most GWAS signals still lack compelling explanatory genes. We have leveraged an adipose-specific gene regulatory network to infer expression regulator activities and phenotypic master regulators (MRs), which were used to detect activity QTLs (aQTLs) at cardiometabolic trait GWAS loci. Regulator activities were inferred with the VIPER algorithm that integrates enrichment of expected expression changes among a regulator's target genes with confidence in their regulator-target network interactions and target overlap between different regulators (i.e., pleiotropy). Phenotypic MRs were identified as those regulators whose activities were most important in predicting their respective phenotypes using random forest modeling. While eQTLs were typically more significant than aQTLs in cis, the opposite was true among candidate MRs in trans. Several GWAS loci colocalized with MR trans-eQTLs/aQTLs in the absence of colocalized cis-QTLs. Intriguingly, at the 1p36.1 BMI GWAS locus the EPHB2 cis-aQTL was stronger than its cis-eQTL and colocalized with the GWAS signal and 35 BMI MR trans-aQTLs, suggesting the GWAS signal may be mediated by effects on EPHB2 activity and its downstream effects on a network of BMI MRs. These MR and aQTL analyses represent systems genetic methods that may be broadly applied to supplement standard eQTL analyses for suggesting molecular effects mediating GWAS signals. Author summary: Most human genetic variants lie outside of genes (the functional units of the genome that directly affect a cell's biology) making it unclear which genes are responsible for influencing their associated traits. The gold-standard for linking genetic variants to genes is expression QTL (or eQTL) analysis, which tests for associations between genetic variants and the expression of genes. However, this approach often fails to identify gene(s) potentially mediating the effects of trait-associated variants. Here we propose the use of a supplementary approach called activity QTL (or aQTL) analysis using existing eQTL data. We first inferred the activities of genes that affect other genes' expression based on a gene regulatory network and then tested associations between genetic variants and these inferred regulator activities. This can be advantageous when a gene's measured expression level is a poor indicator of its downstream activity or when multiple genetic influences are funneled through key regulators in a gene regulatory network to affect the trait of interest. Using this approach, we identified genes expressed in adipose tissue (i.e., fat) potentially mediating genetic effects on BMI, fat distribution, diabetes risk and blood cholesterol levels. More broadly, this work highlights the benefits of leveraging relational (i.e., topological) information in addressing complex biological problems. [ABSTRACT FROM AUTHOR]

The growing number of next-generation sequencing (NGS) data presents a unique opportunity to study the combined impact of mitochondrial and nuclear-encoded genetic variation in complex disease. Mitochondrial DNA variants and in particular, heteroplasmic variants, are critical for determining human disease severity. While there are approaches for obtaining mitochondrial DNA variants from NGS data, these software do not account for the unique characteristics of mitochondrial genetics and can be inaccurate even for homoplasmic variants. We introduce MitoScape, a novel, big-data, software for extracting mitochondrial DNA sequences from NGS. MitoScape adopts a novel departure from other algorithms by using machine learning to model the unique characteristics of mitochondrial genetics. We also employ a novel approach of using rho-zero (mitochondrial DNA-depleted) data to model nuclear-encoded mitochondrial sequences. We showed that MitoScape produces accurate heteroplasmy estimates using gold-standard mitochondrial DNA data. We provide a comprehensive comparison of the most common tools for obtaining mtDNA variants from NGS and showed that MitoScape had superior performance to compared tools in every statistically category we compared, including false positives and false negatives. By applying MitoScape to common disease examples, we illustrate how MitoScape facilitates important heteroplasmy-disease association discoveries by expanding upon a reported association between hypertrophic cardiomyopathy and mitochondrial haplogroup T in men (adjusted p-value = 0.003). The improved accuracy of mitochondrial DNA variants produced by MitoScape will be instrumental in diagnosing disease in the context of personalized medicine and clinical diagnostics. Author summary: Recent studies have highlighted the importance of mitochondrial DNA variation in both primary mitochondrial disease and complex, human pathology including COVID-19, and space-flight stress. The vast amount of existing, next-generation sequencing (NGS) data can be leveraged to interrogate both nuclear and mitochondrial DNA (mtDNA) sequence simultaneously, allowing for analysis of the interplay between mitochondrial and nuclear encoded genes in mitochondrial function. Identifying mtDNA sequence accurately is complicated by the presence of nuclear encoded mitochondrial sequences (NUMTs), which are homologous to mtDNA. Current software for analyzing mtDNA from NGS do not accurately model the unique characteristics of mitochondrial genetics. We introduce MitoScape, a novel, big-data, software which models mitochondrial genetics through machine learning to accurately identify mtDNA sequence from NGS data. MitoScape takes advantage of rho-zero cell data to model the characteristics of NUMTs. We show that MitoScape produces more accurate heteroplasmy estimates compared to published software. We provide an example of applying MitoScape in replicating an association between hypertrophic cardiomyopathy and mitochondrial haplogroup T in men. MitoScape is an important contribution to mitochondrial genomics allowing for accurate mtDNA variants, and the ability to tailor mtDNA analysis in different population and disease contexts, which is not available in other software. [ABSTRACT FROM AUTHOR]

Transcription factors (TFs) are proteins that promote or reduce the expression of genes by binding short genomic DNA sequences known as transcription factor binding sites (TFBS). While several tools have been developed to scan for potential occurrences of TFBS in linear DNA sequences or reference genomes, no tool exists to find them in pangenome variation graphs (VGs). VGs are sequence-labelled graphs that can efficiently encode collections of genomes and their variants in a single, compact data structure. Because VGs can losslessly compress large pangenomes, TFBS scanning in VGs can efficiently capture how genomic variation affects the potential binding landscape of TFs in a population of individuals. Here we present GRAFIMO (GRAph-based Finding of Individual Motif Occurrences), a command-line tool for the scanning of known TF DNA motifs represented as Position Weight Matrices (PWMs) in VGs. GRAFIMO extends the standard PWM scanning procedure by considering variations and alternative haplotypes encoded in a VG. Using GRAFIMO on a VG based on individuals from the 1000 Genomes project we recover several potential binding sites that are enhanced, weakened or missed when scanning only the reference genome, and which could constitute individual-specific binding events. GRAFIMO is available as an open-source tool, under the MIT license, at https://github.com/pinellolab/GRAFIMO and https://github.com/InfOmics/GRAFIMO. Author summary: Transcription factors (TFs) are key regulatory proteins and mutations occurring in their binding sites can alter the normal transcriptional landscape of a cell and lead to disease states. Pangenome variation graphs (VGs) efficiently encode genomes from a population of individuals and their genetic variations. GRAFIMO is an open-source tool that extends the traditional PWM scanning procedure to VGs. By scanning for potential TBFS in VGs, GRAFIMO can simultaneously search thousands of genomes while accounting for SNPs, indels, and structural variants. GRAFIMO reports motif occurrences, their statistical significance, frequency, and location within the reference or alternative haplotypes in a given VG. GRAFIMO makes it possible to study how genetic variation affects the binding landscape of known TFs within a population of individuals. [ABSTRACT FROM AUTHOR]

Increasing body of experimental evidence suggests that anticancer and antimicrobial therapies may themselves promote the acquisition of drug resistance by increasing mutability. The successful control of evolving populations requires that such biological costs of control are identified, quantified and included to the evolutionarily informed treatment protocol. Here we identify, characterise and exploit a trade-off between decreasing the target population size and generating a surplus of treatment-induced rescue mutations. We show that the probability of cure is maximized at an intermediate dosage, below the drug concentration yielding maximal population decay, suggesting that treatment outcomes may in some cases be substantially improved by less aggressive treatment strategies. We also provide a general analytical relationship that implicitly links growth rate, pharmacodynamics and dose-dependent mutation rate to an optimal control law. Our results highlight the important, but often neglected, role of fundamental eco-evolutionary costs of control. These costs can often lead to situations, where decreasing the cumulative drug dosage may be preferable even when the objective of the treatment is elimination, and not containment. Taken together, our results thus add to the ongoing criticism of the standard practice of administering aggressive, high-dose therapies and motivate further experimental and clinical investigation of the mutagenicity and other hidden collateral costs of therapies. Author summary: Evolution of drug resistance to anticancer and antimicrobial therapies is widespread among cancer and pathogen cell populations. Classical theory posits strictly that genetic and phenotypic variation is generated in evolving populations independently of the selection pressure. However, recent experimental findings among antimicrobial agents, traditional cytotoxic chemotherapies and targeted cancer therapies suggest that treatment not only imposes selection but can also affect the rate of adaptation by increasing mutability. Here we analyse a model with drug-induced increase in mutation rate and explore its consequences for treatment optimisation. We argue that the true biological cost of treatment is not limited to the harmful side-effects, but instead realises even more profoundly by fundamentally changing the underlying eco-evolutionary dynamics within the microenvironment. Using the concept of evolutionary rescue, we formulate the treatment as an optimal control problem and solve the optimal elimination strategy, which minimises the probability of evolutionary rescue. We show that aggressive elimination strategies, which aim at eradication as fast as possible and which represent the current standard of care, can be detrimental even with modest drug-induced increases (fold change ≤10) to the baseline mutation rate. [ABSTRACT FROM AUTHOR]

The repurposing of biomedical data is inhibited by its fragmented and multi-formatted nature that requires redundant investment of time and resources by data scientists. This is particularly true for Type 1 Diabetes (T1D), one of the most intensely studied common childhood diseases. Intense investigation of the contribution of pancreatic β-islet and T-lymphocytes in T1D has been made. However, genetic contributions from B-lymphocytes, which are known to play a role in a subset of T1D patients, remain relatively understudied. We have addressed this issue through the creation of Biomedical Data Commons (BMDC), a knowledge graph that integrates data from multiple sources into a single queryable format. This increases the speed of analysis by multiple orders of magnitude. We develop a pipeline using B-lymphocyte multi-dimensional epigenome and connectome data and deploy BMDC to assess genetic variants in the context of Type 1 Diabetes (T1D). Pipeline-identified variants are primarily common, non-coding, poorly conserved, and are of unknown clinical significance. While variants and their chromatin connectivity are cell-type specific, they are associated with well-studied disease genes in T-lymphocytes. Candidates include established variants in the HLA-DQB1 and HLA-DRB1 and IL2RA loci that have previously been demonstrated to protect against T1D in humans and mice providing validation for this method. Others are included in the well-established T1D GRS2 genetic risk scoring method. More intriguingly, other prioritized variants are completely novel and form the basis for future mechanistic and clinical validation studies The BMDC community-based platform can be expanded and repurposed to increase the accessibility, reproducibility, and productivity of biomedical information for diverse applications including the prioritization of cell type-specific disease alleles from complex phenotypes. Author summary: The fragmentation of datasets prevents repurposing due to time-intensive data cleaning and joining. This is especially true for Type 1 Diabetes for which the genetic contributions from B-lymphocytes, a specific type of white blood cells, remain understudied. Here, we create Biomedical Data Commons (BMDC), a knowledge graph, which maps datasets to common entities making them easy to search using queries. We also built a genetic variant prioritization pipeline that uses multi-dimensional 'omics data including three-dimensional connectome data. Using B-lymphocyte cell-type specific data as input, we prioritized variants associated with Type 1 Diabetes. The candidate variants identified are primarily of unknown clinical significance and in the non-coding genome. They are also connected with genes previously implicated in Type 1 Diabetes, suggesting that they affect cell type-specific gene regulation. Some variants in the HLA and IL2RA locus, which are important genomic regions for regulation of immune function, have previously been validated in humans and mice. Other variants have been included in a well-established Type 1 Diabetes genetic risk scoring method. This validates our approach and highlights the novel variants identified that should be prioritized for future clinical and experimental validation. BMDC is a community-based platform that increases the accessibility, reproducibility, and productivity of biomedical information for diverse applications, and our approach is widely applicable for prioritizing variants from other complex diseases. [ABSTRACT FROM AUTHOR]

Despite the growing constellation of genetic loci linked to common traits, these loci have yet to account for most heritable variation, and most act through poorly understood mechanisms. Recent machine learning (ML) systems have used hierarchical biological knowledge to associate genetic mutations with phenotypic outcomes, yielding substantial predictive power and mechanistic insight. Here, we use an ontology-guided ML system to map single nucleotide variants (SNVs) focusing on 6 classic phenotypic traits in natural yeast populations. The 29 identified loci are largely novel and account for ~17% of the phenotypic variance, versus <3% for standard genetic analysis. Representative results show that sensitivity to hydroxyurea is linked to SNVs in two alternative purine biosynthesis pathways, and that sensitivity to copper arises through failure to detoxify reactive oxygen species in fatty acid metabolism. This work demonstrates a knowledge-based approach to amplifying and interpreting signals in population genetic studies. Author summary: Genome-wide association studies (GWAS) have identified many important loci for common diseases and other traits. However, the loci identified by these studies are almost always many steps away from an understanding of underlying biological mechanisms. Here we develop an approach using hierarchical biological knowledge to identify genes and pathways responsible for phenotypic traits. Variants identified by the new method could explain a substantially greater fraction of heritability than previously reported. Moreover, we identified mechanistic pathways by which each causal variant affects cellular function. For example, we find that sensitivity to hydroxyurea is tied to genetic variants in two alternative purine biosynthesis pathways, and that sensitivity to copper arises through failure to detoxify reactive oxygen species in fatty acid metabolism. The new approach is a potentially transformative concept for understanding the genetic drivers of phenotypic variance, with potential applications in understanding traits in biomedicine and agriculture. [ABSTRACT FROM AUTHOR]

Ankyrin protein repeats bind to a wide range of substrates and are one of the most common protein motifs in nature. Here, we collate a high-quality alignment of 7,407 ankyrin repeats and examine for the first time, the distribution of human population variants from large-scale sequencing of healthy individuals across this family. Population variants are not randomly distributed across the genome but are constrained by gene essentiality and function. Accordingly, we interpret the population variants in context with evolutionary constraint and structural features including secondary structure, accessibility and protein-protein interactions across 383 three-dimensional structures of ankyrin repeats. We find five positions that are highly conserved across homologues and also depleted in missense variants within the human population. These positions are significantly enriched in intra-domain contacts and so likely to be key for repeat packing. In contrast, a group of evolutionarily divergent positions are found to be depleted in missense variants in human but significantly enriched in protein-protein interactions. Our analysis also suggests the domain has three, not two surfaces, each with different patterns of enrichment in protein-substrate interactions and missense variants. Our findings will be of interest to those studying or engineering ankyrin-repeat containing proteins as well as those interpreting the significance of disease variants. Author summary: Comparison of variation at each position of the amino acid sequence for a protein across different species is a powerful way to identify parts of the protein that are important for its structure and function. Large-scale DNA sequencing of healthy people has recently made it possible to study normal genetic variation within just one species. Our work combines information on genetic differences between over 100,000 people with in-depth analysis of all available three-dimensional structures for Ankyrin repeats, which are a widespread family of binding proteins formed by units with similar amino acid sequence that are found in tandem. Our combined analysis identifies sites critical for ankyrin stability as well as the positions most important for substrate interactions and hence function. Although focused only on the Ankyrins, the principles developed in our work are general and can be applied to any protein family. [ABSTRACT FROM AUTHOR]

There is an abundance of malaria genetic data being collected from the field, yet using these data to understand the drivers of regional epidemiology remains a challenge. A key issue is the lack of models that relate parasite genetic diversity to epidemiological parameters. Classical models in population genetics characterize changes in genetic diversity in relation to demographic parameters, but fail to account for the unique features of the malaria life cycle. In contrast, epidemiological models, such as the Ross-Macdonald model, capture malaria transmission dynamics but do not consider genetics. Here, we have developed an integrated model encompassing both parasite evolution and regional epidemiology. We achieve this by combining the Ross-Macdonald model with an intra-host continuous-time Moran model, thus explicitly representing the evolution of individual parasite genomes in a traditional epidemiological framework. Implemented as a stochastic simulation, we use the model to explore relationships between measures of parasite genetic diversity and parasite prevalence, a widely-used metric of transmission intensity. First, we explore how varying parasite prevalence influences genetic diversity at equilibrium. We find that multiple genetic diversity statistics are correlated with prevalence, but the strength of the relationships depends on whether variation in prevalence is driven by host- or vector-related factors. Next, we assess the responsiveness of a variety of statistics to malaria control interventions, finding that those related to mixed infections respond quickly (∼months) whereas other statistics, such as nucleotide diversity, may take decades to respond. These findings provide insights into the opportunities and challenges associated with using genetic data to monitor malaria epidemiology. Author Summary: Knowledge of how the prevalence of P.falciparum malaria varies, either between regions or through time, is critical to the operation of malaria control programs. Yet obtaining this information through traditional methods presents many challenges. Parasite genetic data is increasingly accessible, and may provide an alternative means to estimate P.falciparum prevalence in the field. However, our understanding of how the genetic diversity of parasite populations relates to prevalence is limited, and suitable models to guide our understanding are largely lacking. Here, we merge two classical models—the Ross-Macondald and the Moran—to produce a framework in which the relationships between parasite genetic diversity and prevalence can be explored. We find that several genetic diversity statistics are correlated with prevalence, although to differing degrees, and over different time scales. Overall, statistics related to mixed infection are robustly and rapidly responsive to changes in prevalence, suggesting they may be a useful focal point for the development of malaria surveillance methods that harness genetic data. [ABSTRACT FROM AUTHOR]

Non-linear MR (Testing for non-linear causal effects using a binary genotype in a Mendelian randomization study: application to alcohol and cardiovascular traits. These different methods include, but are not limited to, random effects MR, the MR modal-based estimator [[23]], weighted median MR [[24]], MR-Egger regression [[25]], MR-robust adjusted profile score (RAPS) [[26]], and simulation-based heterogeneity and outlier tests (e.g., MR-pleiotropy residual sum and outlier (PRESSO)) [[27]]. Mendelian randomization (MR) is an epidemiological technique for estimating causal relationships using observational data, which has become very popular in recent years following publication of a seminal article by Smith and Ebrahim in 2003 [[1]]. [Extracted from the article]

Ethanol is one of the most widely used recreational substances in the world and due to its ubiquitous use, ethanol abuse has been the cause of over 3.3 million deaths each year. In addition to its effects, ethanol's primary metabolite, acetaldehyde, is a carcinogen that can cause symptoms of facial flushing, headaches, and nausea. How strongly ethanol or acetaldehyde affects an individual depends highly on the genetic polymorphisms of certain genes. In particular, the genetic polymorphisms of mitochondrial aldehyde dehydrogenase, ALDH2, play a large role in the metabolism of acetaldehyde. Thus, it is important to characterize how genetic variations can lead to different exposures and responses to ethanol and acetaldehyde. While the pharmacokinetics of ethanol metabolism through alcohol dehydrogenase have been thoroughly explored in previous studies, in this paper, we combined a base physiologically-based pharmacokinetic (PBPK) model with a whole-body genome-scale model (WBM) to gain further insight into the effect of other less explored processes and genetic variations on ethanol metabolism. This combined model was fit to clinical data and used to show the effect of alcohol concentrations, organ damage, ALDH2 enzyme polymorphisms, and ALDH2-inhibiting drug disulfiram on ethanol and acetaldehyde exposure. Through estimating the reaction rates of auxiliary processes with dynamic Flux Balance Analysis, The PBPK-WBM was able to navigate around a lack of kinetic constants traditionally associated with PK modelling and demonstrate the compensatory effects of the body in response to decreased liver enzyme expression. Additionally, the model demonstrated that acetaldehyde exposure increased with higher dosages of disulfiram and decreased ALDH2 efficiency, and that moderate consumption rates of ethanol could lead to unexpected accumulations in acetaldehyde. This modelling framework combines the comprehensive steady-state analyses from genome-scale models with the dynamics of traditional PK models to create a highly personalized form of PBPK modelling that can push the boundaries of precision medicine. Author summary: Alcohol is a widely used recreational drug in many parts of the world and it is often abused or misused, leading to the deaths of millions of people each year from driving under the influence and overdose. Additionally, the body breaks down alcohol into acetaldehyde, a carcinogen that has its own effects ranging from headaches and nausea to liver damage. The effects of ethanol and acetaldehyde vary due to genetic variations that create different forms of the enzymes responsible for breaking them down. Due to these differences, it is important to characterize how these changes affect the metabolism of alcohol and acetaldehyde. To capture these differences, we have created a new model that integrates the traditional pharmacokinetic model with a whole-body genome-scale model that can characterize different genetic variations. In addition, traditional models often require experimentally measured data, yet with this new framework we avoid this tedious process by mathematically solving the genome-scale model with the dynamic Flux Balance Analysis technique, allowing for gap filling. Through this model, we show that the whole-body genome-scale model demonstrates flexibility and robustness that has not been seen before in pharmacokinetic models. Our model combines advantages from both pharmacokinetic and genome-scale modelling and can be personalized to characterize individual reactions to other drugs and further precision medicine. [ABSTRACT FROM AUTHOR]

Driven by the necessity to survive environmental pathogens, the human immune system has evolved exceptional diversity and plasticity, to which several factors contribute including inheritable structural polymorphism of the underlying genes. Characterizing this variation is challenging due to the complexity of these loci, which contain extensive regions of paralogy, segmental duplication and high copy-number repeats, but recent progress in long-read sequencing and optical mapping techniques suggests this problem may now be tractable. Here we assess this by using long-read sequencing platforms from PacBio and Oxford Nanopore, supplemented with short-read sequencing and Bionano optical mapping, to sequence DNA extracted from CD14+ monocytes and peripheral blood mononuclear cells from a single European individual identified as HV31. We use this data to build a de novo assembly of eight genomic regions encoding four key components of the immune system, namely the human leukocyte antigen, immunoglobulins, T cell receptors, and killer-cell immunoglobulin-like receptors. Validation of our assembly using k-mer based and alignment approaches suggests that it has high accuracy, with estimated base-level error rates below 1 in 10 kb, although we identify a small number of remaining structural errors. We use the assembly to identify heterozygous and homozygous structural variation in comparison to GRCh38. Despite analyzing only a single individual, we find multiple large structural variants affecting core genes at all three immunoglobulin regions and at two of the three T cell receptor regions. Several of these variants are not accurately callable using current algorithms, implying that further methodological improvements are needed. Our results demonstrate that assessing haplotype variation in these regions is possible given sufficiently accurate long-read and associated data. Continued reductions in the cost of these technologies will enable application of these methods to larger samples and provide a broader catalogue of germline structural variation at these loci, an important step toward making these regions accessible to large-scale genetic association studies. Author summary: The human immune system is incredibly versatile underlying its capacity to defend the body against thousands of pathogens. At a molecular level, it recognizes pathogens using large libraries of antibodies and related protein receptors. These molecules are encoded by gene families that are particularly difficult to analyze due to their unusually complex patterns of similarities and differences between genes and individuals. To overcome this, we applied several sequencing methods to DNA from a single individual and developed methods to reconstruct the underlying sequence at eight of the immune-associated regions. Importantly, we used DNA extracted from monocytes to avoid capturing the further rearrangements that occur in active immune cells. We generated accurate assemblies by integrating multiple complementary data types, although we noted a small subset of locations that remain challenging. Moreover, we found that this individual contains multiple structural differences between the two inherited chromosomes and compared to previously analyzed genomes, affecting the copy number of immune system genes. Application of these methods in larger numbers of individuals will clearly uncover much more variation than is currently known, and might lead to new understanding of the effect of genetic variation on the broad range of human diseases determined by the immune response. [ABSTRACT FROM AUTHOR]

It is of great interest and potential to discover causal relationships between pairs of exposures and outcomes using genetic variants as instrumental variables (IVs) to deal with hidden confounding in observational studies. Two most popular approaches are Mendelian randomization (MR), which usually use independent genetic variants/SNPs across the genome, and transcriptome-wide association studies (TWAS) (or their generalizations) using cis-SNPs local to a gene (or some genome-wide and likely dependent SNPs), as IVs. In spite of their many promising applications, both approaches face a major challenge: the validity of their causal conclusions depends on three critical assumptions on valid IVs, and more generally on other modeling assumptions, which however may not hold in practice. The most likely as well as challenging situation is due to the wide-spread horizontal pleiotropy, leading to two of the three IV assumptions being violated and thus to biased statistical inference. More generally, we'd like to conduct a goodness-of-fit (GOF) test to check the model being used. Although some methods have been proposed as being robust to various degrees to the violation of some modeling assumptions, they often give different and even conflicting results due to their own modeling assumptions and possibly lower statistical efficiency, imposing difficulties to the practitioner in choosing and interpreting varying results across different methods. Hence, it would help to directly test whether any assumption is violated or not. In particular, there is a lack of such tests for TWAS. We propose a new and general GOF test, called TEDE (TEsting Direct Effects), applicable to both correlated and independent SNPs/IVs (as commonly used in TWAS and MR respectively). Through simulation studies and real data examples, we demonstrate high statistical power and advantages of our new method, while confirming the frequent violation of modeling (including valid IV) assumptions in practice and thus the importance of model checking by applying such a test in MR/TWAS analysis. Author summary: With the increasing availability of large-scale GWAS summary data of various complex traits/diseases and software packages, it has become convenient and popular to apply Mendelian randomization (MR) and transcriptome-wide association studies (TWAS), using genetic variants as instrumental variables (IVs), to address fundamental and significant questions by unraveling causal relationships between complex or molecular traits such as gene expression and other complex traits. However, the validity of such causal conclusions critically depends on the validity of the model being used, including three key IV assumptions. In particular, with the wide-spread horizontal pleiotropy of genetic variants, two of the three IV assumptions may be violated, leading to biased inference from MR and TWAS. This issue may become more severe as more trait-associated genetic variants are used as IVs to increase the power of MR and TWAS. Although there are some methods to check the modeling assumptions for MR with independent genetic variants as IVs, there is barely any powerful one for TWAS (or more generally for MR and similar methods) with correlated SNPs as IVs. We propose such a powerful method applicable to both MR and TWAS with local or genome-wide, possibly correlated/dependent, SNPs as IVs, demonstrating its higher statistical power than several commonly used methods, while confirming the frequent violation of modeling/IV assumptions in TWAS with our example GWAS data of schizophrenia, Alzheimer's disease and blood lipids. An important conclusion is that in practice it is necessary to conduct model checking in MR and TWAS, and our proposed method is expected to be useful for such a task. [ABSTRACT FROM AUTHOR]

The relationship between different levels of integration is a key feature for understanding the genotype-phenotype map. Here, we describe a novel method of integrated data analysis that incorporates protein abundance data into constraint-based modeling to elucidate the biological mechanisms underlying phenotypic variation. Specifically, we studied yeast genetic diversity at three levels of phenotypic complexity in a population of yeast obtained by pairwise crosses of eleven strains belonging to two species, Saccharomyces cerevisiae and S. uvarum. The data included protein abundances, integrated traits (life-history/fermentation) and computational estimates of metabolic fluxes. Results highlighted that the negative correlation between production traits such as population carrying capacity (K) and traits associated with growth and fermentation rates (Jmax) is explained by a differential usage of energy production pathways: a high K was associated with high TCA fluxes, while a high Jmax was associated with high glycolytic fluxes. Enrichment analysis of protein sets confirmed our results. This powerful approach allowed us to identify the molecular and metabolic bases of integrated trait variation, and therefore has a broad applicability domain. Author summary: The integration of data at different levels of cellular organization is an important goal in computational biology for understanding the way the genotypic variation translates into phenotypic variation. Novel profiling technologies and accurate high-throughput phenotyping now allows genomic, transcriptomic, metabolic and proteomic characterization of a large number of individuals under various environmental conditions. However, the metabolic fluxes remain difficult to measure. In this work, we take advantage of recent advances in genome-scale functional annotation and constraint-based metabolic modeling to provide a mathematical framework that allows to estimate internal cellular fluxes from protein abundances and elucidate the biological mechanisms underlying phenotypic variation. Applied to yeast as a model system, this approach highlights that the negative correlation between production traits such as maximum population size and growth and fermentation traits is explained by a differential usage of energy production pathways. The ability to identify molecular and metabolic bases of the variation of integrated traits through population studies has a broad applicability domain. [ABSTRACT FROM AUTHOR]

While advancements in genome sequencing have identified millions of somatic mutations in cancer, their functional impact is poorly understood. We previously developed the expression-based variant impact phenotyping (eVIP) method to use gene expression data to characterize the function of gene variants. The eVIP method uses a decision tree-based algorithm to predict the functional impact of somatic variants by comparing gene expression signatures induced by introduction of wild-type (WT) versus mutant cDNAs in cell lines. The method distinguishes between variants that are gain-of-function, loss-of-function, change-of-function, or neutral. We present eVIP2, software that allows for pathway analysis (eVIP Pathways) and usage with RNA-seq data. To demonstrate the eVIP2 software and approach, we characterized two recurrent frameshift variants in RNF43, a negative regulator of Wnt signaling, frequently mutated in colorectal, gastric, and endometrial cancer. RNF43 WT, RNF43 R117fs, RNF43 G659fs, or GFP control cDNA were overexpressed in HEK293T cells. Analysis with eVIP2 predicted that the frameshift at position 117 was a loss-of-function mutation, as expected. The second frameshift at position 659 has been previously described as a passenger mutation that maintains the RNF43 WT function as a negative regulator of Wnt. Surprisingly, eVIP2 predicted G659fs to be a change-of-function mutation. Additional eVIP Pathways analysis of RNF43 G659fs predicted 10 pathways to be significantly altered, including TNF-α via NFκB signaling, KRAS signaling, and hypoxia, highlighting the benefit of a more comprehensive approach when determining the impact of gene variant function. To validate these predictions, we performed reporter assays and found that each pathway activated by expression of RNF43 G659fs, but not expression of RNF43 WT, was identified as impacted by eVIP2, supporting that RNF43 G659fs is a change-of-function mutation and its effect on the identified pathways. Pathway activation was further validated by Western blot analysis. Lastly, we show primary colon adenocarcinoma patient samples with R117fs and G659fs variants have transcriptional profiles similar to BRAF missense mutations with activated RAS/MAPK signaling, consistent with KRAS signaling pathways being GOF in both variants. The eVIP2 method is an important step towards overcoming the current challenge of variant interpretation in the implementation of precision medicine. eVIP2 is available at https://github.com/BrooksLabUCSC/eVIP2. Author summary: Cancer is often caused by a DNA change, called a mutation, that creates an atypical protein variant. Many cancer-associated mutations remain uninvestigated because experimental validation of their effect on gene function is costly, time-consuming, and requires prior knowledge of a gene's function. We present our computational tool, eVIP2, which advances our previous approach called expression-based variant impact phenotyping (eVIP) and can be used with RNA sequencing data. This method uses gene expression data to characterize a gene variant's function, requiring no prior knowledge of the wild-type gene's function. With eVIP2, we can predict if a mutation causes a gain, loss, or change in function, or if it is neutral. Here, we found two recurrent frameshift mutations in RNF43 that have different effects on gene function where one mutation causes multiple cancer pathways to be activated. This study highlights the usefulness of using the eVIP2 tool for profiling a mutation's impact at the pathway level. [ABSTRACT FROM AUTHOR]

Plasmodium species, the causative agent of malaria, have a complex life cycle involving two hosts. The sporozoite life stage is characterized by an extended phase in the mosquito salivary glands followed by free movement and rapid invasion of hepatocytes in the human host. This transmission stage has been the subject of many transcriptomics and proteomics studies and is also targeted by the most advanced malaria vaccine. We applied Bayesian data integration to determine which proteins are not only present in sporozoites but are also specific to that stage. Transcriptomic and proteomic Plasmodium data sets from 26 studies were weighted for how representative they are for sporozoites, based on a carefully assembled gold standard for Plasmodium falciparum (Pf) proteins known to be present or absent during the sporozoite life stage. Of 5418 Pf genes for which expression data were available at the RNA level or at the protein level, 975 were identified as enriched in sporozoites and 90 specific to them. We show that Pf sporozoites are enriched for proteins involved in type II fatty acid synthesis in the apicoplast and GPI anchor synthesis, but otherwise appear metabolically relatively inactive in the salivary glands of mosquitos. Newly annotated hypothetical sporozoite-specific and sporozoite-enriched proteins highlight sporozoite-specific functions. They include PF3D7_0104100 that we identified to be homologous to the prominin family, which in human has been related to a quiescent state of cancer cells. We document high levels of genetic variability for sporozoite proteins, specifically for sporozoite-specific proteins that elicit antibodies in the human host. Nevertheless, we can identify nine relatively well-conserved sporozoite proteins that elicit antibodies and that together can serve as markers for previous exposure. Our understanding of sporozoite biology benefits from identifying key pathways that are enriched during this life stage. This work can guide studies of molecular mechanisms underlying sporozoite biology and potential well-conserved targets for marker and drug development. Author summary: When a person is bitten by an infectious malaria mosquito, sporozoites are injected into the skin with mosquito saliva. These sporozoites then travel to the liver, invade hepatocytes and multiply before the onset of the symptom-causing blood stage of malaria. By integrating published data, we contrast sporozoite protein expression with other life stages to filter out the unique features of sporozoites that help us understand this stage. We used a "guideline" that we derived from the literature on individual proteins so that we knew which proteins should be present or absent at the sporozoite stage, allowing us to weigh 26 data sets for their relevance to sporozoites. Among the newly discovered sporozoite-specific genes are candidates for fatty acid synthesis while others might play a role keeping the sporozoites in an inactive state in the mosquito salivary glands. Furthermore, we show that most sporozoite-specific proteins are genetically more variable than non-sporozoite proteins. We identify a set of conserved sporozoite proteins against which antibodies can serve as markers of recent exposure to sporozoites or that can serve as vaccine candidates. Our predictions of sporozoite-specific proteins and the assignment of previously unknown functions give new insights into the biology of this life stage. [ABSTRACT FROM AUTHOR]

The development of high-throughput next-generation sequencing technologies and large-scale genetic association studies produced numerous advances in the biostatistics field. Various aggregation tests, i.e. statistical methods that analyze associations of a trait with multiple markers within a genomic region, have produced a variety of novel discoveries. Notwithstanding their usefulness, there is no single test that fits all needs, each suffering from specific drawbacks. Selecting the right aggregation test, while considering an unknown underlying genetic model of the disease, remains an important challenge. Here we propose a new ensemble method, called Excalibur, based on an optimal combination of 36 aggregation tests created after an in-depth study of the limitations of each test and their impact on the quality of result. Our findings demonstrate the ability of our method to control type I error and illustrate that it offers the best average power across all scenarios. The proposed method allows for novel advances in Whole Exome/Genome sequencing association studies, able to handle a wide range of association models, providing researchers with an optimal aggregation analysis for the genetic regions of interest., Competing Interests: The authors have declared that no competing interests exist., (Copyright: © 2023 Boutry et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.)