Seed design is important for sequence similarity search applications such as read mapping and average nucleotide identity (ANI) estimation. Although k-mers and spaced k-mers are likely the most well-known and used seeds, sensitivity suffers at high error rates, particularly when indels are present. Recently, we developed a pseudorandom seeding construct, strobemers, which was empirically shown to have high sensitivity also at high indel rates. However, the study lacked a deeper understanding of why. In this study, we propose a model to estimate the entropy of a seed and find that seeds with high entropy, according to our model, in most cases have high match sensitivity. Our discovered seed randomness-sensitivity relationship explains why some seeds perform better than others, and the relationship provides a framework for designing even more sensitive seeds. We also present three new strobemer seed constructs: mixedstrobes, altstrobes, and multistrobes. We use both simulated and biological data to show that our new seed constructs improve sequence-matching sensitivity to other strobemers. We show that the three new seed constructs are useful for read mapping and ANI estimation. For read mapping, we implement strobemers into minimap2 and observe 30% faster alignment time and 0.2% higher accuracy than using k-mers when mapping reads at high error rates. As for ANI estimation, we find that higher entropy seeds have a higher rank correlation between estimated and true ANI.

The pairwise sequentially Markovian coalescent (PSMC) algorithm and its extensions infer the coalescence time of two homologous chromosomes at each genomic position. This inference is used in reconstructing demographic histories, detecting selection signatures, studying genome-wide associations, constructing ancestral recombination graphs, and more. Inference of coalescence times between each pair of haplotypes in a large data set is of great interest, as they may provide rich information about the population structure and history of the sample. Here, we introduce a new method, Gamma-SMC, which is more than 10 times faster than current methods. To obtain this speed-up, we represent the posterior coalescence time distributions succinctly as a gamma distribution with just two parameters; in contrast, PSMC and its extensions hold these in a vector over discrete intervals of time. Thus, Gamma-SMC has constant time-complexity per site, without dependence on the number of discrete time states. Additionally, because of this continuous representation, our method is able to infer times spanning many orders of magnitude and, as such, is robust to parameter misspecification. We describe how this approach works, show its performance on simulated and real data, and illustrate its use in studying recent positive selection in the 1000 Genomes Project data set.

Tools that classify sequencing reads against a database of reference sequences require efficient index data-structures. The r-index is a compressed full-text index that answers substring presence/absence, count, and locate queries in space proportional to the amount of distinct sequence in the database: [Formula: see text] space, where r is the number of Burrows-Wheeler runs. To date, the r-index has lacked the ability to quickly classify matches according to which reference sequences (or sequence groupings, i.e., taxa) a match overlaps. We present new algorithms and methods for solving this problem. Specifically, given a collection D of d documents, [Formula: see text] over an alphabet of size σ, we extend the r-index with [Formula: see text] additional words to support document listing queries for a pattern [Formula: see text] that occurs in [Formula: see text] documents in D in [Formula: see text] time and [Formula: see text] space, where w is the machine word size. Applied in a bacterial mock community experiment, our method is up to three times faster than a comparable method that uses the standard r-index locate queries. We show that our method classifies both simulated and real nanopore reads at the strain level with higher accuracy compared with other approaches. Finally, we present strategies for compacting this structure in applications in which read lengths or match lengths can be bounded.

Hummingbirds are very well adapted to sustain efficient and rapid metabolic shifts. They oxidize ingested nectar to directly fuel flight when foraging but have to switch to oxidizing stored lipids derived from ingested sugars during the night or long-distance migratory flights. Understanding how this organism moderates energy turnover is hampered by a lack of information regarding how relevant enzymes differ in sequence, expression, and regulation. To explore these questions, we generated a chromosome-scale genome assembly of the ruby-throated hummingbird (A. colubris) using a combination of long- and short-read sequencing, scaffolding it using existing assemblies. We then used hybrid long- and short-read RNA sequencing of liver and muscle tissue in fasted and fed metabolic states for a comprehensive transcriptome assembly and annotation. Our genomic and transcriptomic data found positive selection of key metabolic genes in nectivorous avian species and deletion of critical genes (SLC2A4, GCK) involved in glucostasis in other vertebrates. We found expression of a fructose-specific version of SLC2A5 putatively in place of insulin-sensitive SLC2A5, with predicted protein models suggesting affinity for both fructose and glucose. Alternative isoforms may even act to sequester fructose to preclude limitations from transport in metabolism. Finally, we identified differentially expressed genes from fasted and fed hummingbirds, suggesting key pathways for the rapid metabolic switch hummingbirds undergo.

Genome sequence data are no longer scarce. The UK Biobank alone comprises 200,000 individual genomes, with more on the way, leading the field of human genetics toward sequencing entire populations. Within the next decades, other model organisms will follow suit, especially domesticated species such as crops and livestock. Having sequences from most individuals in a population will present new challenges for using these data to improve health and agriculture in the pursuit of a sustainable future. Existing population genetic methods are designed to model hundreds of randomly sampled sequences but are not optimized for extracting the information contained in the larger and richer data sets that are beginning to emerge, with thousands of closely related individuals. Here we develop a new method called trio-based inference of dominance and selection (TIDES) that uses data from tens of thousands of family trios to make inferences about natural selection acting in a single generation. TIDES further improves on the state of the art by making no assumptions regarding demography, linkage, or dominance. We discuss how our method paves the way for studying natural selection from new angles.

High-throughput methods such as RNA-seq, ChIP-seq, and ATAC-seq have well-established guidelines, commercial kits, and analysis pipelines that enable consistency and wider adoption for understanding genome function and regulation. STARR-seq, a popular assay for directly quantifying the activities of thousands of enhancer sequences simultaneously, has seen limited standardization across studies. The assay is long, with more than 250 steps, and frequent customization of the protocol and variations in bioinformatics methods raise concerns for reproducibility of STARR-seq studies. Here, we assess each step of the protocol and analysis pipelines from published sources and in-house assays, and identify critical steps and quality control (QC) checkpoints necessary for reproducibility of the assay. We also provide guidelines for experimental design, protocol scaling, customization, and analysis pipelines for better adoption of the assay. These resources will allow better optimization of STARR-seq for specific research needs, enable comparisons and integration across studies, and improve the reproducibility of results.

Bacterial genome data are accumulating at an unprecedented speed due to the routine use of sequencing in clinical diagnoses, public health surveillance, and population genetics studies. Genealogical reconstruction is fundamental to many of these uses; however, inferring genealogy from large-scale genome data sets quickly, accurately, and flexibly is still a challenge. Here, we extend an alignment- and annotation-free method, PopPUNK, to increase its flexibility and interpretability across data sets. Our method, iterative-PopPUNK, rapidly produces multiple consistent cluster assignments across a range of sequence identities. By constructing a partially resolved genealogical tree with respect to these clusters, users can select a resolution most appropriate for their needs. We showed the accuracy of clusters at all levels of similarity and genealogical inference of iterative-PopPUNK based on simulated data and obtained phylogenetically concordant results in real data sets from seven bacterial species. Using two example sets of Escherichia/Shigella and Vibrio parahaemolyticus genomes, we show that iterative-PopPUNK can achieve cluster resolutions ranging from phylogroup down to sequence typing (ST). The iterative-PopPUNK algorithm is implemented in the "PopPUNK_iterate" program, available as part of the PopPUNK package.

Sequence-to-graph alignment is crucial for applications such as variant genotyping, read error correction, and genome assembly. We propose a novel seeding approach that relies on long inexact matches rather than short exact matches, and show that it yields a better time-accuracy trade-off in settings with up to a [Formula: see text] mutation rate. We use sketches of a subset of graph nodes, which are more robust to indels, and store them in a k-nearest neighbor index to avoid the curse of dimensionality. Our approach contrasts with existing methods and highlights the important role that sketching into vector space can play in bioinformatics applications. We show that our method scales to graphs with 1 billion nodes and has quasi-logarithmic query time for queries with an edit distance of [Formula: see text] For such queries, longer sketch-based seeds yield a [Formula: see text] increase in recall compared with exact seeds. Our approach can be incorporated into other aligners, providing a novel direction for sequence-to-graph alignment.

Gene expression data provide molecular insights into the functional impact of genetic variation, for example, through expression quantitative trait loci (eQTLs). With an improving understanding of the association between genotypes and gene expression comes a greater concern that gene expression profiles could be matched to genotype profiles of the same individuals in another data set, known as a linking attack. Prior works show such a risk could analyze only a fraction of eQTLs that is independent owing to restrictive model assumptions, leaving the full extent of this risk incompletely understood. To address this challenge, we introduce the discriminative sequence model (DSM), a novel probabilistic framework for predicting a sequence of genotypes based on gene expression data. By modeling the joint distribution over all known eQTLs in a genomic region, DSM improves the power of linking attacks with necessary calibration for linkage disequilibrium and redundant predictive signals. We show greater linking accuracy of DSM compared with existing approaches across a range of attack scenarios and data sets including up to 22,288 individuals, suggesting that DSM helps uncover a substantial additional risk overlooked by previous studies. Our work provides a unified framework for assessing the privacy risks of sharing diverse omics data sets beyond transcriptomics.

Tandem repeats (TRs) are one of the largest sources of polymorphism, and their length is associated with gene regulation. Although previous studies reported several tandem repeats regulating gene splicing in cis (spl-TRs), no large-scale study has been conducted. In this study, we established a genome-wide catalog of 9537 spl-TRs with a total of 58,290 significant TR-splicing associations across 49 tissues (false discovery rate 5%) by using Genotype-Tissue expression (GTex) Project data. Regression models explaining splicing variation by using spl-TRs and other flanking variants suggest that at least some of the spl-TRs directly modulate splicing. In our catalog, two spl-TRs are known loci for repeat expansion diseases, spinocerebellar ataxia 6 (SCA6) and 12 (SCA12). Splicing alterations by these spl-TRs were compatible with those observed in SCA6 and SCA12. Thus, our comprehensive spl-TR catalog may help elucidate the pathomechanism of genetic diseases.

Recombinant plasmid vectors are versatile tools that have facilitated discoveries in molecular biology, genetics, proteomics, and many other fields. As the enzymatic and bacterial processes used to create recombinant DNA can introduce errors, sequence validation is an essential step in plasmid assembly. Sanger sequencing is the current standard for plasmid validation; however, this method is limited by an inability to sequence through complex secondary structure and lacks scalability when applied to full-plasmid sequencing of multiple plasmids owing to read-length limits. Although high-throughput sequencing does provide full-plasmid sequencing at scale, it is impractical and costly when used outside of library-scale validation. Here, we present Oxford nanopore-based rapid analysis of multiplexed plasmids (OnRamp), an alternative method for routine plasmid validation that combines the advantages of high-throughput sequencing's full-plasmid coverage and scalability with Sanger's affordability and accessibility by leveraging nanopore's long-read sequencing technology. We include customized wet-laboratory protocols for plasmid preparation along with a pipeline designed for analysis of read data obtained using these protocols. This analysis pipeline is deployed on the OnRamp web app, which generates alignments between actual and predicted plasmid sequences, quality scores, and read-level views. OnRamp is designed to be broadly accessible regardless of programming experience to facilitate more widespread adoption of long-read sequencing for routine plasmid validation. Here we describe the OnRamp protocols and pipeline and show our ability to obtain full sequences from pooled plasmids while detecting sequence variation even in regions of high secondary structure at less than half the cost of equivalent Sanger sequencing.

A genomic sketch is a small, probabilistic representation of the set of k-mers in a sequencing data set. Sketches are building blocks for large-scale analyses that consider similarities between many pairs of sequences or sequence collections. Although existing tools can easily compare tens of thousands of genomes, data sets can reach millions of sequences and beyond. Popular tools also fail to consider k-mer multiplicities, making them less applicable in quantitative settings. Here, we describe a method called Dashing 2 that builds on the SetSketch data structure. SetSketch is related to HyperLogLog (HLL) but discards use of leading zero count in favor of a truncated logarithm of adjustable base. Unlike HLL, SetSketch can perform multiplicity-aware sketching when combined with the ProbMinHash method. Dashing 2 integrates locality-sensitive hashing to scale all-pairs comparisons to millions of sequences. It achieves superior similarity estimates for the Jaccard coefficient and average nucleotide identity compared with the original Dashing, but in much less time while using the same-sized sketch. Dashing 2 is a free, open source software.

Tandem duplications are frequent structural variations of the genome and play important roles in genetic disease and cancer. However, interpreting the phenotypic consequences of tandem duplications remains challenging, in part owing to the lack of genetic tools to model such variations. Here, we developed a strategy, tandem duplication via prime editing (TD-PE), to create targeted, programmable, and precise tandem duplication in the mammalian genome. In this strategy, we design a pair of in trans prime editing guide RNAs (pegRNAs) for each targeted tandem duplication, which encode the same edits but prime the single-stranded DNA (ssDNA) extension in opposite directions. The reverse transcriptase (RT) template of each extension is designed homologous to the target region of the other single guide RNA (sgRNA) to promote the reannealing of the edited DNA strands and the duplication of the fragment in between. We showed that TD-PE produced robust and precise in situ tandem duplications of genomic fragments ranging from ∼50 bp to ∼10 kb, with a maximal efficiency up to 28.33%. By fine-tuning the pegRNAs, we achieved simultaneous targeted duplication and fragment insertion. Finally, we successfully produced multiple disease-relevant tandem duplications, showing the general utility of TD-PE in genetic research.

We developed an analysis pipeline that can extract microbial sequences from spatial transcriptomic (ST) data and assign taxonomic labels, generating a spatial microbial abundance matrix in addition to the default host expression matrix, enabling simultaneous analysis of host expression and microbial distribution. We called the pipeline spatial metatranscriptome (SMT) and applied it on both human and murine intestinal sections and validated the spatial microbial abundance information with alternative assays. Biological insights were gained from these novel data that showed host-microbe interaction at various spatial scales. Finally, we tested experimental modification that can increase microbial capture while preserving host spatial expression quality and, by use of positive controls, quantitatively showed the capture efficiency and recall of our methods. This proof-of-concept work shows the feasibility of SMT analysis and paves the way for further experimental optimization and application.

The recently constructed mutant libraries of diploid crops by the CRISPR-Cas9 system have provided abundant resources for functional genomics and crop breeding. However, because of the genome complexity, it is a big challenge to accomplish large-scale targeted mutagenesis in polyploid plants. Here, we demonstrate the feasibility of using a pooled CRISPR library to achieve genome-scale targeted editing in an allotetraploid crop of Brassica napus A total of 18,414 sgRNAs were designed to target 10,480 genes of interest, and afterward, 1104 regenerated transgenic plants harboring 1088 sgRNAs were obtained. Editing interrogation results revealed that 93 of the 178 genes were identified as mutated, thus representing an editing efficiency of 52.2%. Furthermore, we have discovered that Cas9-mediated DNA cleavages tend to occur at all the target sites guided by the same individual sgRNA, a novel finding in polyploid plants. Finally, we show the strong capability of reverse genetic screening for various traits with the postgenotyped plants. Several genes, which might dominate the fatty acid profile and seed oil content and have yet to be reported, were unveiled from the forward genetic studies. Our research provides valuable resources for functional genomics, elite crop breeding, and a good reference for high-throughput targeted mutagenesis in other polyploid plants.

Genome browsers have become an intuitive and critical tool to visualize and analyze genomic features and data. Conventional genome browsers display data/annotations on a single reference genome/assembly; there are also genomic alignment viewer/browsers that help users visualize alignment, mismatch, and rearrangement between syntenic regions. However, there is a growing need for a comparative epigenome browser that can display genomic and epigenomic data sets across different species and enable users to compare them between syntenic regions. Here, we present the WashU Comparative Epigenome Browser. It allows users to load functional genomic data sets/annotations mapped to different genomes and display them over syntenic regions simultaneously. The browser also displays genetic differences between the genomes from single-nucleotide variants (SNVs) to structural variants (SVs) to visualize the association between epigenomic differences and genetic differences. Instead of anchoring all data sets to the reference genome coordinates, it creates independent coordinates of different genome assemblies to faithfully present features and data mapped to different genomes. It uses a simple, intuitive genome-align track to illustrate the syntenic relationship between different species. It extends the widely used WashU Epigenome Browser infrastructure and can be expanded to support multiple species. This new browser function will greatly facilitate comparative genomic/epigenomic research, as well as support the recent growing needs to directly compare and benchmark the T2T CHM13 assembly and other human genome assemblies.

Understanding the impact of DNA variation on human traits is a fundamental question in human genetics. Variable number tandem repeats (VNTRs) make up ∼3% of the human genome but are often excluded from association analysis owing to poor read mappability or divergent repeat content. Although methods exist to estimate VNTR length from short-read data, it is known that VNTRs vary in both length and repeat (motif) composition. Here, we use a repeat-pangenome graph (RPGG) constructed on 35 haplotype-resolved assemblies to detect variation in both VNTR length and repeat composition. We align population-scale data from the Genotype-Tissue Expression (GTEx) Consortium to examine how variations in sequence composition may be linked to expression, including cases independent of overall VNTR length. We find that 9422 out of 39,125 VNTRs are associated with nearby gene expression through motif variations, of which only 23.4% are accessible from length. Fine-mapping identifies 174 genes to be likely driven by variation in certain VNTR motifs and not overall length. We highlight two genes, CACNA1C and RNF213, that have expression associated with motif variation, showing the utility of RPGG analysis as a new approach for trait association in multiallelic and highly variable loci.

The zebrafish telencephalon is composed of highly specialized subregions that regulate complex behaviors such as learning, memory, and social interactions. The transcriptional signatures of the neuronal cell types in the telencephalon and the timeline of their emergence from larva to adult remain largely undescribed. Using an integrated analysis of single-cell transcriptomes of approximately 64,000 cells obtained from 6-day-postfertilization (dpf), 15-dpf, and adult telencephalon, we delineated nine main neuronal cell types in the pallium and eight in the subpallium and nominated novel marker genes. Comparing zebrafish and mouse neuronal cell types revealed both conserved and absent types and marker genes. Mapping of cell types onto a spatial larval reference atlas created a resource for anatomical and functional studies. Using this multiage approach, we discovered that although most neuronal subtypes are established early in the 6-dpf fish, some emerge or expand in number later in development. Analyzing the samples from each age separately revealed further complexity in the data, including several cell types that expand substantially in the adult forebrain and do not form clusters at the larval stages. Together, our work provides a comprehensive transcriptional analysis of the cell types in the zebrafish telencephalon and a resource for dissecting its development and function.

Cancer results from an evolutionary process that typically yields multiple clones with varying sets of mutations within the same tumor. Accurately modeling this process is key to understanding and predicting cancer evolution. Here, we introduce clone to mutation (CloMu), a flexible and low-parameter tree generative model of cancer evolution. CloMu uses a two-layer neural network trained via reinforcement learning to determine the probability of new mutations based on the existing mutations on a clone. CloMu supports several prediction tasks, including the determination of evolutionary trajectories, tree selection, causality and interchangeability between mutations, and mutation fitness. Importantly, previous methods support only some of these tasks, and many suffer from overfitting on data sets with a large number of mutations. Using simulations, we show that CloMu either matches or outperforms current methods on a wide variety of prediction tasks. In particular, for simulated data with interchangeable mutations, current methods are unable to uncover causal relationships as effectively as CloMu. On breast cancer and leukemia cohorts, we show that CloMu determines similarities and causal relationships between mutations as well as the fitness of mutations. We validate CloMu's inferred mutation fitness values for the leukemia cohort by comparing them to clonal proportion data not used during training, showing high concordance. In summary, CloMu's low-parameter model facilitates a wide range of prediction tasks regarding cancer evolution on increasingly available cohort-level data sets.

Multiple sequence alignment (MSA) is a critical step in the study of protein sequence and function. Typically, MSA algorithms progressively align pairs of sequences and combine these alignments with the aid of a guide tree. These alignment algorithms use scoring systems based on substitution matrices to measure amino acid similarities. Although successful, standard methods struggle on sets of proteins with low sequence identity: the so-called twilight zone of protein alignment. For these difficult cases, another source of information is needed. Protein language models are a powerful new approach that leverages massive sequence data sets to produce high-dimensional contextual embeddings for each amino acid in a sequence. These embeddings have been shown to reflect physicochemical and higher-order structural and functional attributes of amino acids within proteins. Here, we present a novel approach to MSA, based on clustering and ordering amino acid contextual embeddings. Our method for aligning semantically consistent groups of proteins circumvents the need for many standard components of MSA algorithms, avoiding initial guide tree construction, intermediate pairwise alignments, gap penalties, and substitution matrices. The added information from contextual embeddings leads to higher accuracy alignments for structurally similar proteins with low amino-acid similarity. We anticipate that protein language models will become a fundamental component of the next generation of algorithms for generating MSAs.

HIV-1 integration introduces ectopic transcription factor binding sites into host chromatin. We postulate that the integrated provirus serves as an ectopic enhancer that recruits additional transcription factors to the integration locus, increases chromatin accessibility, changes 3D chromatin interactions, and enhances both retroviral and host gene expression. We used four well-characterized HIV-1-infected cell line clones having unique integration sites and low to high levels of HIV-1 expression. Using single-cell DOGMA-seq, which captured the heterogeneity of HIV-1 expression and host chromatin accessibility, we found that HIV-1 transcription correlated with HIV-1 accessibility and host chromatin accessibility. HIV-1 integration increased local host chromatin accessibility within an ∼5- to 30-kb distance. CRISPRa- and CRISPRi-mediated HIV-1 promoter activation and inhibition confirmed integration site-dependent HIV-1-driven changes of host chromatin accessibility. HIV-1 did not drive chromatin confirmation changes at the genomic level (by Hi-C) or the enhancer connectome (by H3K27ac HiChIP). Using 4C-seq to interrogate HIV-1-chromatin interactions, we found that HIV-1 interacted with host chromatin ∼100-300 kb from the integration site. By identifying chromatin regions having both increased transcription factor activity (by ATAC-seq) and HIV-1-chromatin interaction (by 4C-seq), we identified enrichment of ETS, RUNT, and ZNF-family transcription factor binding that may mediate HIV-1-host chromatin interactions. Our study has found that HIV-1 promoter activity increases host chromatin accessibility, and HIV-1 interacted with host chromatin within the existing chromatin boundaries in an integration site-dependent manner.

Rare species are vital members of a microbial community, but retrieving their genomes is difficult because of their low abundance. The ReadUntil (RU) approach allows nanopore devices to sequence specific DNA molecules selectively in real time, which provides an opportunity for enriching rare species. Despite the robustness of enriching rare species by reducing the sequencing depth of known host sequences, such as the human genome, there is still a gap in RU-based enriching of rare species in environmental samples whose community composition is unclear, and many rare species have poor or incomplete reference genomes in public databases. Therefore, here we present metaRUpore to overcome this challenge. When we applied metaRUpore to a thermophilic anaerobic digester (TAD) community and human gut microbial community, it reduced coverage of the high-abundance populations and modestly increased (∼2×) the genome coverage of the rare taxa, facilitating successful recovery of near-finished metagenome-assembled genomes (nf-MAGs) of rare species. The simplicity and robustness of the approach make it accessible for laboratories with moderate computational resources, and hold the potential to become the standard practice in future metagenomic sequencing of complicated microbiomes.

Tn5 transposon tagments double-stranded DNA and RNA/DNA hybrids to generate nucleic acids that are ready to be amplified for high-throughput sequencing. The nucleic acid substrates for the Tn5 transposon must be explored to increase the applications of Tn5. Here, we found that the Tn5 transposon can transpose oligos into the 5' end of single-stranded DNA longer than 140 nucleotides. Based on this property of Tn5, we developed a tagmentation-based and ligation-enabled single-stranded DNA sequencing method called TABLE-seq. Through a series of reaction temperature, time, and enzyme concentration tests, we applied TABLE-seq to strand-specific RNA sequencing, starting with as little as 30 pg of total RNA. Moreover, compared with traditional dUTP-based strand-specific RNA sequencing, this method detects more genes, has a higher strand specificity, and shows more evenly distributed reads across genes. Together, our results provide insights into the properties of Tn5 transposons and expand the applications of Tn5 in cutting-edge sequencing techniques.

High-throughput genotyping enables the large-scale analysis of genetic diversity in population genomics and genome-wide association studies that combine the genotypic and phenotypic characterization of large collections of accessions. Sequencing-based approaches for genotyping are progressively replacing traditional genotyping methods because of the lower ascertainment bias. However, genome-wide genotyping based on sequencing becomes expensive in species with large genomes and a high proportion of repetitive DNA. Here we describe the use of CRISPR-Cas9 technology to deplete repetitive elements in the 3.76-Gb genome of lentil (Lens culinaris), 84% consisting of repeats, thus concentrating the sequencing data on coding and regulatory regions (single-copy regions). We designed a custom set of 566,766 gRNAs targeting 2.9 Gbp of repeats and excluding repetitive regions overlapping annotated genes and putative regulatory elements based on ATAC-seq data. The novel depletion method removed ∼40% of reads mapping to repeats, increasing those mapping to single-copy regions by ∼2.6-fold. When analyzing 25 million fragments, this repeat-to-single-copy shift in the sequencing data increased the number of genotyped bases of ∼10-fold compared to nondepleted libraries. In the same condition, we were also able to identify ∼12-fold more genetic variants in the single-copy regions and increased the genotyping accuracy by rescuing thousands of heterozygous variants that otherwise would be missed because of low coverage. The method performed similarly regardless of the multiplexing level, type of library or genotypes, including different cultivars and a closely related species (L. orientalis). Our results showed that CRISPR-Cas9-driven repeat depletion focuses sequencing data on single-copy regions, thus improving high-density and genome-wide genotyping in large and repetitive genomes.

Macaques provide the most widely used nonhuman primate models for studying the immunology and pathogenesis of human diseases. Although the macaque major histocompatibility complex (MHC) region shares most features with the human leukocyte antigen (HLA) region, macaques have an expanded repertoire of MHC class I genes. Although a chimera of two rhesus macaque MHC haplotypes was first published in 2004, the structural diversity of MHC genomic organization in macaques remains poorly understood owing to a lack of adequate genomic reference sequences. We used ultralong Oxford Nanopore and high-accuracy Pacific Biosciences (PacBio) HiFi sequences to fully assemble the ∼5.2-Mb M3 haplotype of an MHC-homozygous, Mauritian-origin cynomolgus macaque (Macaca fascicularis). The MHC homozygosity allowed us to assemble a single MHC haplotype unambiguously and avoid chimeric assemblies that hampered previous efforts to characterize this exceptionally complex genomic region in macaques. The high quality of this new assembly is exemplified by the identification of an extended cluster of six Mafa-AG genes that contains a recent duplication with a highly similar ∼48.5-kb block of sequence. The MHC class II region of this M3 haplotype is similar to the previously sequenced rhesus macaque haplotype and HLA class II haplotypes. The MHC class I region, in contrast, contains 13 MHC-B genes, four MHC-A genes, and three MHC-E genes (vs. 19 MHC-B, two MHC-A, and one MHC-E in the previously sequenced haplotype). These results provide an unambiguously assembled single contiguous cynomolgus macaque MHC haplotype with fully curated gene annotations that will inform infectious disease and transplantation research.

Compacted de Bruijn graphs are one of the most fundamental data structures in computational genomics. Colored compacted de Bruijn graphs are a variant built on a collection of sequences and associate to each k-mer the sequences in which it appears. We present GGCAT, a tool for constructing both types of graphs, based on a new approach merging the k-mer counting step with the unitig construction step, as well as on numerous practical optimizations. For compacted de Bruijn graph construction, GGCAT achieves speed-ups of 3× to 21× compared with the state-of-the-art tool Cuttlefish 2. When constructing the colored variant, GGCAT achieves speed-ups of 5× to 39× compared with the state-of-the-art tool BiFrost. Additionally, GGCAT is up to 480× faster than BiFrost for batch sequence queries on colored graphs.

The genomes of maize and other eukaryotes contain stable haplotypes in regions of low recombination. These regions, including centromeres, long heterochromatic blocks, and rDNA arrays, have been difficult to analyze with respect to their diversity and origin. Greatly improved genome assemblies are now available that enable comparative genomics over these and other nongenic spaces. Using 26 complete maize genomes, we developed methods to align intergenic sequences while excluding genes and regulatory regions. The centromere haplotypes (cenhaps) extend for megabases on either side of the functional centromere regions and appear as evolutionary strata, with haplotype divergence/coalescence times dating as far back as 450 thousand years ago (kya). Application of the same methods to other low recombination regions (heterochromatic knobs and rDNA) and all intergenic spaces revealed that deep coalescence times are ubiquitous across the maize pan-genome. Divergence estimates vary over a broad timescale with peaks at ∼16 and 300 kya, reflecting a complex history of gene flow among diverging populations and changes in population size associated with domestication. Cenhaps and other long haplotypes provide vivid displays of this ancient diversity.

Alternative polyadenylation (APA) enables a gene to generate multiple transcripts with different 3' ends, which is dynamic across different cell types or conditions. Many computational methods have been developed to characterize sample-specific APA using the corresponding RNA-seq data, but suffered from high error rate on both polyadenylation site (PAS) identification and quantification of PAS usage (PAU), and bias toward 3' untranslated regions. Here we developed a tool for APA identification and quantification (APAIQ) from RNA-seq data, which can accurately identify PAS and quantify PAU in a transcriptome-wide manner. Using 3' end-seq data as the benchmark, we showed that APAIQ outperforms current methods on PAS identification and PAU quantification, including DaPars2, Aptardi, mountainClimber, SANPolyA, and QAPA. Finally, applying APAIQ on 421 RNA-seq samples from liver cancer patients, we identified >540 tumor-associated APA events and experimentally validated two intronic polyadenylation candidates, demonstrating its capacity to unveil cancer-related APA with a large-scale RNA-seq data set.

Large-scale high-throughput sequencing data sets have been transformative for informing clinical variant interpretation and for use as reference panels for statistical and population genetic efforts. Although such resources are often treated as ground truth, we find that in widely used reference data sets such as the Genome Aggregation Database (gnomAD), some variants pass gold-standard filters, yet are systematically different in their genotype calls across genotype discovery approaches. The inclusion of such discordant sites in study designs involving multiple genotype discovery strategies could bias results and lead to false-positive hits in association studies owing to technological artifacts rather than a true relationship to the phenotype. Here, we describe this phenomenon of discordant genotype calls across genotype discovery approaches, characterize the error mode of wrong calls, provide a list of discordant sites identified in gnomAD that should be treated with caution in analyses, and present a metric and machine learning classifier trained on gnomAD data to identify likely discordant variants in other data sets. We find that different genotype discovery approaches have different sets of variants at which this problem occurs, but there are characteristic variant features that can be used to predict discordant behavior. Discordant sites are largely shared across ancestry groups, although different populations are powered for the discovery of different variants. We find that the most common error mode is that of a variant being heterozygous for one approach and homozygous for the other, with heterozygous in the genomes and homozygous reference in the exomes making up the majority of miscalls.

Density separation is a process routinely used to segregate minerals, organic matter, and even microplastics, from soils and sediments. Here we apply density separation to archaeological bone powders before DNA extraction to increase endogenous DNA recovery relative to a standard control extraction of the same powders. Using nontoxic heavy liquid solutions, we separated powders from the petrous bones of 10 individuals of similar archaeological preservation into eight density intervals (2.15 to 2.45 g/cm3, in 0.05 increments). We found that the 2.30 to 2.35 g/cm3 and 2.35 to 2.40 g/cm3 intervals yielded up to 5.28-fold more endogenous unique DNA than the corresponding standard extraction (and up to 8.53-fold before duplicate read removal), while maintaining signals of ancient DNA authenticity and not reducing library complexity. Although small 0.05 g/cm3 intervals may maximally optimize yields, a single separation to remove materials with a density above 2.40 g/cm3 yielded up to 2.57-fold more endogenous DNA on average, which enables the simultaneous separation of samples that vary in preservation or in the type of material analyzed. While requiring no new ancient DNA laboratory equipment and fewer than 30 min of extra laboratory work, the implementation of density separation before DNA extraction can substantially boost endogenous DNA yields without decreasing library complexity. Although subsequent studies are required, we present theoretical and practical foundations that may prove useful when applied to other ancient DNA substrates such as teeth, other bones, and sediments.

The information about when and where each gene is to be expressed is mainly encoded in the DNA sequence of enhancers, sequence elements that comprise binding sites (motifs) for different transcription factors (TFs). Most of the research on enhancer sequences has been focused on TF motif presence, whereas the enhancer syntax, that is, the flexibility of important motif positions and how the sequence context modulates the activity of TF motifs, remains poorly understood. Here, we explore the rules of enhancer syntax by a two-pronged approach in Drosophila melanogaster S2 cells: we (1) replace important TF motifs by all possible 65,536 eight-nucleotide-long sequences and (2) paste eight important TF motif types into 763 positions within 496 enhancers. These complementary strategies reveal that enhancers display constrained sequence flexibility and the context-specific modulation of motif function. Important motifs can be functionally replaced by hundreds of sequences constituting several distinct motif types, but these are only a fraction of all possible sequences and motif types. Moreover, TF motifs contribute with different intrinsic strengths that are strongly modulated by the enhancer sequence context (the flanking sequence, the presence and diversity of other motif types, and the distance between motifs), such that not all motif types can work in all positions. The context-specific modulation of motif function is also a hallmark of human enhancers, as we demonstrate experimentally. Overall, these two general principles of enhancer sequences are important to understand and predict enhancer function during development, evolution, and in disease.

Adipocytes contribute to metabolic disorders such as obesity, diabetes, and atherosclerosis. Prior characterizations of the transcriptional network driving adipogenesis have overlooked transiently acting transcription factors (TFs), genes, and regulatory elements that are essential for proper differentiation. Moreover, traditional gene regulatory networks provide neither mechanistic details about individual regulatory element-gene relationships nor temporal information needed to define a regulatory hierarchy that prioritizes key regulatory factors. To address these shortcomings, we integrate kinetic chromatin accessibility (ATAC-seq) and nascent transcription (PRO-seq) data to generate temporally resolved networks that describe TF binding events and resultant effects on target gene expression. Our data indicate which TF families cooperate with and antagonize each other to regulate adipogenesis. Compartment modeling of RNA polymerase density quantifies how individual TFs mechanistically contribute to distinct steps in transcription. The glucocorticoid receptor activates transcription by inducing RNA polymerase pause release, whereas SP and AP-1 factors affect RNA polymerase initiation. We identify Twist2 as a previously unappreciated effector of adipocyte differentiation. We find that TWIST2 acts as a negative regulator of 3T3-L1 and primary preadipocyte differentiation. We confirm that Twist2 knockout mice have compromised lipid storage within subcutaneous and brown adipose tissue. Previous phenotyping of Twist2 knockout mice and Setleis syndrome Twist2 -/- patients noted deficiencies in subcutaneous adipose tissue. This network inference framework is a powerful and general approach for interpreting complex biological phenomena and can be applied to a wide range of cellular processes.

Killer cell immunoglobulin like receptor (KIR) genes and human leukocyte antigen (HLA) genes play important roles in innate and adaptive immunity. They are highly polymorphic and cannot be genotyped with standard variant calling pipelines. Compared with HLA genes, many KIR genes are similar to each other in sequences and may be absent in the chromosomes. Therefore, although many tools have been developed to genotype HLA genes using common sequencing data, none of them work for KIR genes. Even specialized KIR genotypers could not resolve all the KIR genes. Here we describe T1K, a novel computational method for the efficient and accurate inference of KIR or HLA alleles from RNA-seq, whole-genome sequencing, or whole-exome sequencing data. T1K jointly considers alleles across all genotyped genes, so it can reliably identify present genes and distinguish homologous genes, including the challenging KIR2DL5A/KIR2DL5B genes. This model also benefits HLA genotyping, where T1K achieves high accuracy in benchmarks. Moreover, T1K can call novel single-nucleotide variants and process single-cell data. Applying T1K to tumor single-cell RNA-seq data, we found that KIR2DL4 expression was enriched in tumor-specific CD8+ T cells. T1K may open the opportunity for HLA and KIR genotyping across various sequencing applications.

DNA sequencing data continue to progress toward longer reads with increasingly lower sequencing error rates. We focus on the critical problem of mapping, or aligning, low-divergence sequences from long reads (e.g., Pacific Biosciences [PacBio] HiFi) to a reference genome, which poses challenges in terms of accuracy and computational resources when using cutting-edge read mapping approaches that are designed for all types of alignments. A natural idea would be to optimize efficiency with longer seeds to reduce the probability of extraneous matches; however, contiguous exact seeds quickly reach a sensitivity limit. We introduce mapquik, a novel strategy that creates accurate longer seeds by anchoring alignments through matches of k consecutively sampled minimizers (k-min-mers) and only indexing k-min-mers that occur once in the reference genome, thereby unlocking ultrafast mapping while retaining high sensitivity. We show that mapquik significantly accelerates the seeding and chaining steps-fundamental bottlenecks to read mapping-for both the human and maize genomes with [Formula: see text] sensitivity and near-perfect specificity. On the human genome, for both real and simulated reads, mapquik achieves a [Formula: see text] speedup over the state-of-the-art tool minimap2, and on the maize genome, mapquik achieves a [Formula: see text] speedup over minimap2, making mapquik the fastest mapper to date. These accelerations are enabled from not only minimizer-space seeding but also a novel heuristic [Formula: see text] pseudochaining algorithm, which improves upon the long-standing [Formula: see text] bound. Minimizer-space computation builds the foundation for achieving real-time analysis of long-read sequencing data.

Epigenetic modifications undergo drastic erasure and reestablishment after fertilization. This reprogramming is required for proper embryonic development and cell differentiation. In mammals, some histone modifications are not completely reprogrammed and play critical roles in later development. In contrast, in nonmammalian vertebrates, most histone modifications are thought to be more intensively erased and reestablished by the stage of zygotic genome activation (ZGA). However, histone modifications that escape reprogramming in nonmammalian vertebrates and their potential functional roles remain unknown. Here, we quantitatively and comprehensively analyzed histone modification dynamics during epigenetic reprogramming in Japanese killifish, medaka (Oryzias latipes) embryos. Our data revealed that H3K27ac, H3K27me3, and H3K9me3 escape complete reprogramming, whereas H3K4 methylation is completely erased during cleavage stage. Furthermore, we experimentally showed the functional roles of such retained modifications at early stages: (i) H3K27ac premarks promoters during the cleavage stage, and inhibition of histone acetyltransferases disrupts proper patterning of H3K4 and H3K27 methylation at CpG-dense promoters, but does not affect chromatin accessibility after ZGA; (ii) H3K9me3 is globally erased but specifically retained at telomeric regions, which is required for maintenance of genomic stability during the cleavage stage. These results expand the understanding of diversity and conservation of reprogramming in vertebrates, and unveil previously uncharacterized functions of histone modifications retained during epigenetic reprogramming.

Minimizers are ubiquitously used in data structures and algorithms for efficient searching, mapping, and indexing of high-throughput DNA sequencing data. Minimizer schemes select a minimum k-mer in every L-long subsequence of the target sequence, where minimality is with respect to a predefined k-mer order. Commonly used minimizer orders select more k-mers than necessary and therefore provide limited improvement in runtime and memory usage of downstream analysis tasks. The recently introduced universal k-mer hitting sets produce minimizer orders with fewer selected k-mers. Generating compact universal k-mer hitting sets is currently infeasible for k > 13, and thus, they cannot help in the many applications that require minimizer orders for larger k Here, we close the gap of efficient minimizer orders for large values of k by introducing decycling-set-based minimizer orders: new minimizer orders based on minimum decycling sets. We show that in practice these new minimizer orders select a number of k-mers comparable to that of minimizer orders based on universal k-mer hitting sets and can also scale to a larger k Furthermore, we developed a method that computes the minimizers in a sequence on the fly without keeping the k-mers of a decycling set in memory. This enables the use of these minimizer orders for any value of k We expect the new orders to improve the runtime and memory usage of algorithms and data structures in high-throughput DNA sequencing analysis.

For most biological and medical applications of single-cell transcriptomics, an integrative study of multiple heterogeneous single-cell RNA sequencing (scRNA-seq) data sets is crucial. However, present approaches are unable to integrate diverse data sets from various biological conditions effectively because of the confounding effects of biological and technical differences. We introduce single-cell integration (scInt), an integration method based on accurate, robust cell-cell similarity construction and unified contrastive biological variation learning from multiple scRNA-seq data sets. scInt provides a flexible and effective approach to transfer knowledge from the already integrated reference to the query. We show that scInt outperforms 10 other cutting-edge approaches using both simulated and real data sets, particularly in the case of complex experimental designs. Application of scInt to mouse developing tracheal epithelial data shows its ability to integrate development trajectories from different developmental stages. Furthermore, scInt successfully identifies functionally distinct condition-specific cell subpopulations in single-cell heterogeneous samples from a variety of biological conditions.

Topologically associating domains (TADs) have emerged as basic structural and functional units of genome organization and have been determined by many computational methods from Hi-C contact maps. However, the TADs obtained by different methods vary greatly, which makes the accurate determination of TADs a challenging issue and hinders subsequent biological analyses about their organization and functions. Obvious inconsistencies among the TADs identified by different methods indeed make the statistical and biological properties of TADs overly depend on the chosen method rather than on the data. To this end, we use the consensus structural information captured by these methods to define the TAD separation landscape for decoding the consensus domain organization of the 3D genome. We show that the TAD separation landscape could be used to compare domain boundaries across multiple cell types for discovering conserved and divergent topological structures, decipher three types of boundary regions with diverse biological features, and identify consensus TADs (ConsTADs). We illustrate that these analyses could deepen our understanding of the relationships between the topological domains and chromatin states, gene expression, and DNA replication timing.

The collection and sharing of genomic data are becoming increasingly commonplace in research, clinical, and direct-to-consumer settings. The computational protocols typically adopted to protect individual privacy include sharing summary statistics, such as allele frequencies, or limiting query responses to the presence/absence of alleles of interest using web services called Beacons. However, even such limited releases are susceptible to likelihood ratio-based membership-inference attacks. Several approaches have been proposed to preserve privacy, which either suppress a subset of genomic variants or modify query responses for specific variants (e.g., adding noise, as in differential privacy). However, many of these approaches result in a significant utility loss, either suppressing many variants or adding a substantial amount of noise. In this paper, we introduce optimization-based approaches to explicitly trade off the utility of summary data or Beacon responses and privacy with respect to membership-inference attacks based on likelihood ratios, combining variant suppression and modification. We consider two attack models. In the first, an attacker applies a likelihood ratio test to make membership-inference claims. In the second model, an attacker uses a threshold that accounts for the effect of the data release on the separation in scores between individuals in the data set and those who are not. We further introduce highly scalable approaches for approximately solving the privacy-utility tradeoff problem when information is in the form of either summary statistics or presence/absence queries. Finally, we show that the proposed approaches outperform the state of the art in both utility and privacy through an extensive evaluation with public data sets.

Insights into host-virus interactions during SARS-CoV-2 infection are needed to understand COVID-19 pathogenesis and may help to guide the design of novel antiviral therapeutics. N 6-Methyladenosine modification (m6A), one of the most abundant cellular RNA modifications, regulates key processes in RNA metabolism during stress response. Gene expression profiles observed postinfection with different SARS-CoV-2 variants show changes in the expression of genes related to RNA catabolism, including m6A readers and erasers. We found that infection with SARS-CoV-2 variants causes a loss of m6A in cellular RNAs, whereas m6A is detected abundantly in viral RNA. METTL3, the m6A methyltransferase, shows an unusual cytoplasmic localization postinfection. The B.1.351 variant has a less-pronounced effect on METTL3 localization and loss of m6A than did the B.1 and B.1.1.7 variants. We also observed a loss of m6A upon SARS-CoV-2 infection in air/liquid interface cultures of human airway epithelia, confirming that m6A loss is characteristic of SARS-CoV-2-infected cells. Further, transcripts with m6A modification are preferentially down-regulated postinfection. Inhibition of the export protein XPO1 results in the restoration of METTL3 localization, recovery of m6A on cellular RNA, and increased mRNA expression. Stress granule formation, which is compromised by SARS-CoV-2 infection, is restored by XPO1 inhibition and accompanied by a reduced viral infection in vitro. Together, our study elucidates how SARS-CoV-2 inhibits the stress response and perturbs cellular gene expression in an m6A-dependent manner.

The rates of mutation, recombination, and transposition are core parameters in models of evolution. They impact genetic diversity, responses to ongoing selection, and levels of genetic load. However, even for key evolutionary model species such as Drosophila melanogaster and Drosophila simulans, few estimates of these parameters are available, and we have little idea of how rates vary between individuals, sexes, or populations. Knowledge of this variation is fundamental for parameterizing models of genome evolution. Here, we provide direct estimates of mutation, recombination, and transposition rates and their variation in a West African and a European population of D. melanogaster and a European population of D. simulans Across 89 flies, we observe 58 single-nucleotide mutations, 286 crossovers, and 89 transposable element (TE) insertions. Compared to the European D. melanogaster, we find the West African population has a lower mutation rate (1.67 × 10-9 site-1 gen-1 vs. 4.86 × 10-9 site-1 gen-1) and a lower transposition rate (8.99 × 10-5 copy-1 gen-1 vs. 23.36 × 10-5 copy-1 gen-1), but a higher recombination rate (3.44 cM/Mb vs. 2.06 cM/Mb). The European D. simulans population has a similar mutation rate to European D. melanogaster, but a significantly higher recombination rate and a lower, but not significantly different, transposition rate. Overall, we find paternal-derived mutations are more frequent than maternal ones in both species. Our study quantifies the variation in rates of mutation, recombination, and transposition among different populations and sexes, and our direct estimates of these parameters in D. melanogaster and D. simulans will benefit future studies in population and evolutionary genetics.

Approximately 13% of the human genome at certain motifs have the potential to form noncanonical (non-B) DNA structures (e.g., G-quadruplexes, cruciforms, and Z-DNA), which regulate many cellular processes but also affect the activity of polymerases and helicases. Because sequencing technologies use these enzymes, they might possess increased errors at non-B structures. To evaluate this, we analyzed error rates, read depth, and base quality of Illumina, Pacific Biosciences (PacBio) HiFi, and Oxford Nanopore Technologies (ONT) sequencing at non-B motifs. All technologies showed altered sequencing success for most non-B motif types, although this could be owing to several factors, including structure formation, biased GC content, and the presence of homopolymers. Single-nucleotide mismatch errors had low biases in HiFi and ONT for all non-B motif types but were increased for G-quadruplexes and Z-DNA in all three technologies. Deletion errors were increased for all non-B types but Z-DNA in Illumina and HiFi, as well as only for G-quadruplexes in ONT. Insertion errors for non-B motifs were highly, moderately, and slightly elevated in Illumina, HiFi, and ONT, respectively. Additionally, we developed a probabilistic approach to determine the number of false positives at non-B motifs depending on sample size and variant frequency, and applied it to publicly available data sets (1000 Genomes, Simons Genome Diversity Project, and gnomAD). We conclude that elevated sequencing errors at non-B DNA motifs should be considered in low-read-depth studies (single-cell, ancient DNA, and pooled-sample population sequencing) and in scoring rare variants. Combining technologies should maximize sequencing accuracy in future studies of non-B DNA.

The Diversity Outbred (DO) mice and their inbred founders are widely used models of human disease. However, although the genetic diversity of these mice has been well documented, their epigenetic diversity has not. Epigenetic modifications, such as histone modifications and DNA methylation, are important regulators of gene expression and, as such, are a critical mechanistic link between genotype and phenotype. Therefore, creating a map of epigenetic modifications in the DO mice and their founders is an important step toward understanding mechanisms of gene regulation and the link to disease in this widely used resource. To this end, we performed a strain survey of epigenetic modifications in hepatocytes of the DO founders. We surveyed four histone modifications (H3K4me1, H3K4me3, H3K27me3, and H3K27ac), as well as DNA methylation. We used ChromHMM to identify 14 chromatin states, each of which represents a distinct combination of the four histone modifications. We found that the epigenetic landscape is highly variable across the DO founders and is associated with variation in gene expression across strains. We found that epigenetic state imputed into a population of DO mice recapitulated the association with gene expression seen in the founders, suggesting that both histone modifications and DNA methylation are highly heritable mechanisms of gene expression regulation. We illustrate how DO gene expression can be aligned with inbred epigenetic states to identify putative cis-regulatory regions. Finally, we provide a data resource that documents strain-specific variation in the chromatin state and DNA methylation in hepatocytes across nine widely used strains of laboratory mice.

Alternative splicing (AS) regulates gene expression and increases proteomic diversity for the fine tuning of stress responses in plants, but the exact mechanism through which AS functions in plant stress responses is not thoroughly understood. Here, we investigated how AS functions in poplar (Populus trichocarpa), a popular plant for bioremediation, in response to lead (Pb) stress. Using a proteogenomic analysis, we determine that Pb stress induced alterations in AS patterns that are characterized by an increased use of nonconventional splice sites and a higher abundance of Pb-responsive splicing factors (SFs) associated with Pb-responsive transcription factors. A strong Pb(II)-inducible chaperone protein, PtHSP70, that undergoes AS was further characterized. Overexpression of its two spliced isoforms, PtHSP70-AS1 and PtHSP70-AS2, in poplar and Arabidopsis significantly enhances the tolerance to Pb. Further characterization shows that both isoforms can directly bind to Pb(II), and PtHSP70-AS2 exhibits 10-fold higher binding capacities and a greater increase in expression under Pb stress, thereby reducing cellular toxicity through Pb(II) extrusion and conferring Pb tolerance. AS of PtHSP70 is found to be regulated by PtU1-70K, a Pb(II)-inducible core SF involved in 5'-splice site recognition. Because the same splicing pattern is also found in HSP70 orthologs in other plant species, AS of HSP70 may be a common regulatory mechanism to cope with Pb(II) toxicity. Overall, we have revealed a novel post-transcriptional machinery that mediates heavy metal tolerance in diverse plant species. Our findings offer new molecular targets and bioengineering strategies for phytoremediation and provide new insight for future directions in AS research.

In vitro studies indicate the neurodevelopmental disorder gene myelin transcription factor 1-like (MYT1L) suppresses non-neuronal lineage genes during fibroblast-to-neuron direct differentiation. However, MYT1L's molecular and cellular functions in the adult mammalian brain have not been fully characterized. Here, we found that MYT1L loss leads to up-regulated deep layer (DL) gene expression, corresponding to an increased ratio of DL/UL neurons in the adult mouse cortex. To define potential mechanisms, we conducted Cleavage Under Targets & Release Using Nuclease (CUT&RUN) to map MYT1L binding targets and epigenetic changes following MYT1L loss in mouse developing cortex and adult prefrontal cortex (PFC). We found MYT1L mainly binds to open chromatin, but with different transcription factor co-occupancies between promoters and enhancers. Likewise, multiomic data set integration revealed that, at promoters, MYT1L loss does not change chromatin accessibility but increases H3K4me3 and H3K27ac, activating both a subset of earlier neuronal development genes as well as Bcl11b, a key regulator for DL neuron development. Meanwhile, we discovered that MYT1L normally represses the activity of neurogenic enhancers associated with neuronal migration and neuronal projection development by closing chromatin structures and promoting removal of active histone marks. Further, we showed that MYT1L interacts with HDAC2 and transcriptional repressor SIN3B in vivo, providing potential mechanisms underlying repressive effects on histone acetylation and gene expression. Overall, our findings provide a comprehensive map of MYT1L binding in vivo and mechanistic insights into how MYT1L loss leads to aberrant activation of earlier neuronal development programs in the adult mouse brain.

methods are widely used to estimate species trees from genome-scale data. However, they can fail to produce accurate species trees when the input gene trees are highly discordant because of estimation error and biological processes, such as incomplete lineage sorting. Here, we introduce TREE-QMC, a new summary method that offers accuracy and scalability under these challenging scenarios. TREE-QMC builds upon weighted Quartet Max Cut, which takes weighted quartets as input and then constructs a species tree in a divide-and-conquer fashion, at each step forming a graph and seeking its max cut. The wQMC method has been successfully leveraged in the context of species tree estimation by weighting quartets by their frequencies in the gene trees; we improve upon this approach in two ways. First, we address accuracy by normalizing the quartet weights to account for "artificial taxa" introduced during the divide phase so subproblem solutions can be combined during the conquer phase. Second, we address scalability by introducing an algorithm to construct the graph directly from the gene trees; this gives TREE-QMC a time complexity of [Formula: see text], where n is the number of species and k is the number of gene trees, assuming the subproblem decomposition is perfectly balanced. These contributions enable TREE-QMC to be highly competitive in terms of species tree accuracy and empirical runtime with the leading quartet-based methods, even outperforming them on some model conditions explored in our simulation study. We also present the application of these methods to an avian phylogenomics data set.

Mapping the gene targets of chromatin-associated transcription regulators (TRs) is a major goal of genomics research. ChIP-seq of TRs and experiments that perturb a TR and measure the differential abundance of gene transcripts are a primary means by which direct relationships are tested on a genomic scale. It has been reported that there is a poor overlap in the evidence across gene regulation strategies, emphasizing the need for integrating results from multiple experiments. Although research consortia interested in gene regulation have produced a valuable trove of high-quality data, there is an even greater volume of TR-specific data throughout the literature. In this study, we show a workflow for the identification, uniform processing, and aggregation of ChIP-seq and TR perturbation experiments for the ultimate purpose of ranking human and mouse TR-target interactions. Focusing on an initial set of eight regulators (ASCL1, HES1, MECP2, MEF2C, NEUROD1, PAX6, RUNX1, and TCF4), we identified 497 experiments suitable for analysis. We used this corpus to examine data concordance, to identify systematic patterns of the two data types, and to identify putative orthologous interactions between human and mouse. We build upon commonly used strategies to forward a procedure for aggregating and combining these two genomic methodologies, assessing these rankings against independent literature-curated evidence. Beyond a framework extensible to other TRs, our work also provides empirically ranked TR-target listings, as well as transparent experiment-level gene summaries for community use.

Recent studies exploring the impact of methylation in tumor evolution suggest that although the methylation status of many of the CpG sites are preserved across distinct lineages, others are altered as the cancer progresses. Because changes in methylation status of a CpG site may be retained in mitosis, they could be used to infer the progression history of a tumor via single-cell lineage tree reconstruction. In this work, we introduce the first principled distance-based computational method, Sgootr, for inferring a tumor's single-cell methylation lineage tree and for jointly identifying lineage-informative CpG sites that harbor changes in methylation status that are retained along the lineage. We apply Sgootr on single-cell bisulfite-treated whole-genome sequencing data of multiregionally sampled tumor cells from nine metastatic colorectal cancer patients, as well as multiregionally sampled single-cell reduced-representation bisulfite sequencing data from a glioblastoma patient. We show that the tumor lineages constructed reveal a simple model underlying tumor progression and metastatic seeding. A comparison of Sgootr against alternative approaches shows that Sgootr can construct lineage trees with fewer migration events and with more in concordance with the sequential-progression model of tumor evolution, with a running time a fraction of that used in prior studies. Lineage-informative CpG sites identified by Sgootr are in inter-CpG island (CGI) regions, as opposed to intra-CGIs, which have been the main regions of interest in genomic methylation-related analyses.

The reconstruction of phylogenetic networks is an important but challenging problem in phylogenetics and genome evolution, as the space of phylogenetic networks is vast and cannot be sampled well. One approach to the problem is to solve the minimum phylogenetic network problem, in which phylogenetic trees are first inferred, and then the smallest phylogenetic network that displays all the trees is computed. The approach takes advantage of the fact that the theory of phylogenetic trees is mature, and there are excellent tools available for inferring phylogenetic trees from a large number of biomolecular sequences. A tree-child network is a phylogenetic network satisfying the condition that every nonleaf node has at least one child that is of indegree one. Here, we develop a new method that infers the minimum tree-child network by aligning lineage taxon strings in the phylogenetic trees. This algorithmic innovation enables us to get around the limitations of the existing programs for phylogenetic network inference. Our new program, named ALTS, is fast enough to infer a tree-child network with a large number of reticulations for a set of up to 50 phylogenetic trees with 50 taxa that have only trivial common clusters in about a quarter of an hour on average.

Recombination is a key molecular mechanism that has profound implications on both micro- and macroevolutionary processes. However, the determinants of recombination rate variation in holocentric organisms are poorly understood, in particular in Lepidoptera (moths and butterflies). The wood white butterfly (Leptidea sinapis) shows considerable intraspecific variation in chromosome numbers and is a suitable system for studying regional recombination rate variation and its potential molecular underpinnings. Here, we developed a large whole-genome resequencing data set from a population of wood whites to obtain high-resolution recombination maps using linkage disequilibrium information. The analyses revealed that larger chromosomes had a bimodal recombination landscape, potentially caused by interference between simultaneous chiasmata. The recombination rate was significantly lower in subtelomeric regions, with exceptions associated with segregating chromosome rearrangements, showing that fissions and fusions can have considerable effects on the recombination landscape. There was no association between the inferred recombination rate and base composition, supporting a limited influence of GC-biased gene conversion in butterflies. We found significant but variable associations between the recombination rate and the density of different classes of transposable elements, most notably a significant enrichment of short interspersed nucleotide elements in genomic regions with higher recombination rate. Finally, the analyses unveiled significant enrichment of genes involved in farnesyltranstransferase activity in recombination coldspots, potentially indicating that expression of transferases can inhibit formation of chiasmata during meiotic division. Our results provide novel information about recombination rate variation in holocentric organisms and have particular implications for forthcoming research in population genetics, molecular/genome evolution, and speciation.

Structural variations (SVs) are a major contributor to genetic diversity and phenotypic variations, but their prevalence and functions in domestic animals are largely unexplored. Here we generated high-quality genome assemblies for 15 individuals from genetically diverse sheep breeds using Pacific Biosciences (PacBio) high-fidelity sequencing, discovering 130.3 Mb nonreference sequences, from which 588 genes were annotated. A total of 149,158 biallelic insertions/deletions, 6531 divergent alleles, and 14,707 multiallelic variations with precise breakpoints were discovered. The SV spectrum is characterized by an excess of derived insertions compared to deletions (94,422 vs. 33,571), suggesting recent active LINE expansions in sheep. Nearly half of the SVs display low to moderate linkage disequilibrium with surrounding single-nucleotide polymorphisms (SNPs) and most SVs cannot be tagged by SNP probes from the widely used ovine 50K SNP chip. We identified 865 population-stratified SVs including 122 SVs possibly derived in the domestication process among 690 individuals from sheep breeds worldwide. A novel 168-bp insertion in the 5' untranslated region (5' UTR) of HOXB13 is found at high frequency in long-tailed sheep. Further genome-wide association study and gene expression analyses suggest that this mutation is causative for the long-tail trait. In summary, we have developed a panel of high-quality de novo assemblies and present a catalog of structural variations in sheep. Our data capture abundant candidate functional variations that were previously unexplored and provide a fundamental resource for understanding trait biology in sheep.

Because of diverged adaptative phenotypes, fish species of the genus Xiphophorus have contributed to a wide range of research for a century. Existing Xiphophorus genome assemblies are not at the chromosomal level and are prone to sequence gaps, thus hindering advancement of the intra- and inter-species differences for evolutionary, comparative, and translational biomedical studies. Herein, we assembled high-quality chromosome-level genome assemblies for three distantly related Xiphophorus species, namely, X. maculatus, X. couchianus, and X. hellerii Our overall goal is to precisely assess microevolutionary processes in the clade to ascertain molecular events that led to the divergence of the Xiphophorus species and to progress understanding of genetic incompatibility to disease. In particular, we measured intra- and inter-species divergence and assessed gene expression dysregulation in reciprocal interspecies hybrids among the three species. We found expanded gene families and positively selected genes associated with live bearing, a special mode of reproduction. We also found positively selected gene families are significantly enriched in nonpolymorphic transposable elements, suggesting the dispersal of these nonpolymorphic transposable elements has accompanied the evolution of the genes, possibly by incorporating new regulatory elements in support of the Britten-Davidson hypothesis. We characterized inter-specific polymorphisms, structural variants, and polymorphic transposable element insertions and assessed their association to interspecies hybridization-induced gene expression dysregulation related to specific disease states in humans.

SWI/SNF and NuRD are protein complexes that antagonistically regulate DNA accessibility. However, repression of their activities often leads to unanticipated changes in target gene expression (paradoxical), highlighting our incomplete understanding of their activities. Here we show that SWI/SNF and NuRD are in a tug-of-war to regulate PRC2 occupancy at lowly expressed and bivalent genes in mouse embryonic stem cells (mESCs). In contrast, at promoters of average or highly expressed genes, SWI/SNF and NuRD antagonistically modulate RNA polymerase II (Pol II) release kinetics, arguably owing to accompanying alterations in H3.3 and H2A.Z levels at promoter-flanking nucleosomes, leading to paradoxical changes in gene expression. Owing to this mechanism, the relative activities of the two remodelers potentiate gene promoters toward Pol II-dependent open or PRC2-dependent closed chromatin states. Our results highlight RNA Pol II occupancy as the key parameter in determining the direction of gene expression changes in response to SWI/SNF and NuRD inactivation at gene promoters in mESCs.

Across eukaryotes, gene regulation is manifested via chromatin states roughly distinguished as heterochromatin and euchromatin. The establishment, maintenance, and modulation of the chromatin states is mediated using several factors including chromatin modifiers. However, factors that avoid the intrusion of silencing signals into protein-coding genes are poorly understood. Here we show that a plant specific paralog of RNA polymerase (Pol) II, named Pol IV, is involved in avoidance of facultative heterochromatic marks in protein-coding genes, in addition to its well-established functions in silencing repeats and transposons. In its absence, H3K27 trimethylation (me3) mark intruded the protein-coding genes, more profoundly in genes embedded with repeats. In a subset of genes, spurious transcriptional activity resulted in small(s) RNA production, leading to post-transcriptional gene silencing. We show that such effects are significantly pronounced in rice, a plant with a larger genome with distributed heterochromatin compared with Arabidopsis Our results indicate the division of labor among plant-specific polymerases, not just in establishing effective silencing via sRNAs and DNA methylation but also in influencing chromatin boundaries.

Novel genes are essential for evolutionary innovations and differ substantially even between closely related species. Recently, multiple studies across many taxa showed that some novel genes arise de novo, that is, from previously noncoding DNA. To characterize the underlying mutations that allowed de novo gene emergence and their order of occurrence, homologous regions must be detected within noncoding sequences in closely related sister genomes. So far, most studies do not detect noncoding homologs of de novo genes because of incomplete assemblies and annotations, and long evolutionary distances separating genomes. Here, we overcome these issues by searching for de novo expressed open reading frames (neORFs), the not-yet fixed precursors of de novo genes that emerged within a single species. We sequenced and assembled genomes with long-read technology and the corresponding transcriptomes from inbred lines of Drosophila melanogaster, derived from seven geographically diverse populations. We found line-specific neORFs in abundance but few neORFs shared by lines, suggesting a rapid turnover. Gain and loss of transcription is more frequent than the creation of ORFs, for example, by forming new start and stop codons. Consequently, the gain of ORFs becomes rate limiting and is frequently the initial step in neORFs emergence. Furthermore, transposable elements (TEs) are major drivers for intragenomic duplications of neORFs, yet TE insertions are less important for the emergence of neORFs. However, highly mutable genomic regions around TEs provide new features that enable gene birth. In conclusion, neORFs have a high birth-death rate, are rapidly purged, but surviving neORFs spread neutrally through populations and within genomes.

Sketching methods offer computational biologists scalable techniques to analyze data sets that continue to grow in size. MinHash is one such technique to estimate set similarity that has enjoyed recent broad application. However, traditional MinHash has previously been shown to perform poorly when applied to sets of very dissimilar sizes. FracMinHash was recently introduced as a modification of MinHash to compensate for this lack of performance when set sizes differ. This approach has been successfully applied to metagenomic taxonomic profiling in the widely used tool sourmash gather. Although experimental evidence has been encouraging, FracMinHash has not yet been analyzed from a theoretical perspective. In this paper, we perform such an analysis to derive various statistics of FracMinHash, and prove that although FracMinHash is not unbiased (in the sense that its expected value is not equal to the quantity it attempts to estimate), this bias is easily corrected for both the containment and Jaccard index versions. Next, we show how FracMinHash can be used to compute point estimates as well as confidence intervals for evolutionary mutation distance between a pair of sequences by assuming a simple mutation model. We also investigate edge cases in which these analyses may fail to effectively warn the users of FracMinHash indicating the likelihood of such cases. Our analyses show that FracMinHash estimates the containment of a genome in a large metagenome more accurately and more precisely compared with traditional MinHash, and the point estimates and confidence intervals perform significantly better in estimating mutation distances.

Short tandem repeats (STRs) are a class of rapidly mutating genetic elements typically characterized by repeated units of 1-6 bp. We leveraged whole-genome sequencing data for 152 recombinant inbred (RI) strains from the BXD family of mice to map loci that modulate genome-wide patterns of new mutations arising during parent-to-offspring transmission at STRs. We defined quantitative phenotypes describing the numbers and types of germline STR mutations in each strain and performed quantitative trait locus (QTL) analyses for each of these phenotypes. We identified a locus on Chromosome 13 at which strains inheriting the C57BL/6J (B) haplotype have a higher rate of STR expansions than those inheriting the DBA/2J (D) haplotype. The strongest candidate gene in this locus is Msh3, a known modifier of STR stability in cancer and at pathogenic repeat expansions in mice and humans, as well as a current drug target against Huntington's disease. The D haplotype at this locus harbors a cluster of variants near the 5' end of Msh3, including multiple missense variants near the DNA mismatch recognition domain. In contrast, the B haplotype contains a unique retrotransposon insertion. The rate of expansion covaries positively with Msh3 expression-with higher expression from the B haplotype. Finally, detailed analysis of mutation patterns showed that strains carrying the B allele have higher expansion rates, but slightly lower overall total mutation rates, compared with those with the D allele, particularly at tetranucleotide repeats. Our results suggest an important role for inherited variants in Msh3 in modulating genome-wide patterns of germline mutations at STRs.

Small nucleolar RNAs (snoRNAs) are structured noncoding RNAs present in multiple copies within eukaryotic genomes. snoRNAs guide chemical modifications on their target RNA and regulate processes like ribosome assembly and splicing. Most human snoRNAs are embedded within host gene introns, the remainder being independently expressed from intergenic regions. We recently characterized the abundance of snoRNAs and their host gene across several healthy human tissues and found that the level of most snoRNAs does not correlate with that of their host gene, with the observation that snoRNAs embedded within the same host gene often differ drastically in abundance. To better understand the determinants of snoRNA expression, we trained machine learning models to predict whether snoRNAs are expressed or not in human tissues based on more than 30 collected features related to snoRNAs and their genomic context. By interpreting the models' predictions, we find that snoRNAs rely on conserved motifs, a stable global structure and terminal stem, and a transcribed locus to be expressed. We observe that these features explain well the varying abundance of snoRNAs embedded within the same host gene. By predicting the expression status of snoRNAs across several vertebrates, we notice that only one-third of all annotated snoRNAs are expressed per genome, as in humans. Our results suggest that ancestral snoRNAs disseminated within vertebrate genomes, sometimes leading to the development of new functions and a probable gain in fitness and thereby conserving features favorable to the expression of these few snoRNAs, the large remainder often degenerating into pseudogenes.

Seed-chain-extend with k-mer seeds is a powerful heuristic technique for sequence alignment used by modern sequence aligners. Although effective in practice for both runtime and accuracy, theoretical guarantees on the resulting alignment do not exist for seed-chain-extend. In this work, we give the first rigorous bounds for the efficacy of seed-chain-extend with k-mers in expectation Assume we are given a random nucleotide sequence of length ∼n that is indexed (or seeded) and a mutated substring of length ∼m ≤ n with mutation rate θ < 0.206. We prove that we can find a k = Θ(log n) for the k-mer size such that the expected runtime of seed-chain-extend under optimal linear-gap cost chaining and quadratic time gap extension is O(mn f (θ) log n), where f(θ) < 2.43 · θ holds as a loose bound. The alignment also turns out to be good; we prove that more than [Formula: see text] fraction of the homologous bases is recoverable under an optimal chain. We also show that our bounds work when k-mers are sketched, that is, only a subset of all k-mers is selected, and that sketching reduces chaining time without increasing alignment time or decreasing accuracy too much, justifying the effectiveness of sketching as a practical speedup in sequence alignment. We verify our results in simulation and on real noisy long-read data and show that our theoretical runtimes can predict real runtimes accurately. We conjecture that our bounds can be improved further, and in particular, f(θ) can be further reduced.

The mammalian suprachiasmatic nucleus (SCN), located in the ventral hypothalamus, synchronizes and maintains daily cellular and physiological rhythms across the body, in accordance with environmental and visceral cues. Consequently, the systematic regulation of spatiotemporal gene transcription in the SCN is vital for daily timekeeping. So far, the regulatory elements assisting circadian gene transcription have only been studied in peripheral tissues, lacking the critical neuronal dimension intrinsic to the role of the SCN as central brain pacemaker. By using histone-ChIP-seq, we identified SCN-enriched gene regulatory elements that associated with temporal gene expression. Based on tissue-specific H3K27ac and H3K4me3 marks, we successfully produced the first-ever SCN gene-regulatory map. We found that a large majority of SCN enhancers not only show robust 24-h rhythmic modulation in H3K27ac occupancy, peaking at distinct times of day, but also possess canonical E-box (CACGTG) motifs potentially influencing downstream cycling gene expression. To establish enhancer-gene relationships in the SCN, we conducted directional RNA-seq at six distinct times across the day and night, and studied the association between dynamically changing histone acetylation and gene transcript levels. About 35% of the cycling H3K27ac sites were found adjacent to rhythmic gene transcripts, often preceding the rise in mRNA levels. We also noted that enhancers encompass noncoding, actively transcribing enhancer RNAs (eRNAs) in the SCN, which in turn oscillate, along with cyclic histone acetylation, and correlate with rhythmic gene transcription. Taken together, these findings shed light on genome-wide pretranscriptional regulation operative in the central clock that confers its precise and robust oscillation necessary to orchestrate daily timekeeping in mammals.

In this paper, we developed a highly sensitive approach to detect interchromosomal rearrangements in cattle by searching for abnormal linkage disequilibrium patterns between markers located on different chromosomes in large paternal half-sib families genotyped as part of routine genomic evaluations. We screened 5571 families of artificial insemination sires from 15 breeds and revealed 13 putative interchromosomal rearrangements, 12 of which were validated by cytogenetic analysis and long-read sequencing. These consisted of one Robertsonian fusion, 10 reciprocal translocations, and the first case of insertional translocation reported in cattle. Taking advantage of the wealth of data available in cattle, we performed a series of complementary analyses to define the exact nature of these rearrangements, investigate their origins, and search for factors that may have favored their occurrence. We also evaluated the risks to the livestock industry and showed significant negative effects on several traits in the sires and in their balanced or aneuploid progeny compared with wild-type controls. Thus, we present the most comprehensive and thorough screen for interchromosomal rearrangements compatible with normal spermatogenesis in livestock species. This approach is readily applicable to any population that benefits from large genotype data sets, and will have direct applications in animal breeding. Finally, it also offers interesting prospects for basic research by allowing the detection of smaller and rarer types of chromosomal rearrangements than GTG banding, which are interesting models for studying gene regulation and the organization of genome structure.