AbstractThe annotation of genomic ranges of interest represents a recurring task for bioinformatics analyses. These ranges can originate from various sources, including peaks called for transcription factor binding sites (TFBS) or histone modification ChIP-seq experiments, chromatin structure and accessibility experiments (such as ATAC-seq), but also from other types of predictions that result in genomic ranges. While peak annotation primarily driven by ChiP-seq was extensively explored, many approaches remain simplistic (“most closely located TSS”), rely on fixed pre-built references, or require complex scripting tasks on behalf of the user. An adaptable, fast, and universal tool, capable to annotate genomic ranges in the respective biological context is critically missing. UROPA (Universal RObust Peak Annotator) is a command line based tool, intended for universal genomic range annotation. Based on a configuration file, different target features can be prioritized with multiple integrated queries. These can be sensitive for feature type, distance, strand specificity, feature attributes (e.g. protein_coding) or anchor position relative to the feature. UROPA can incorporate reference annotation files (GTF) from different sources (Gencode, Ensembl, RefSeq), as well as custom reference annotation files. Statistics and plots transparently summarize the annotation process. UROPA is implemented in Python and R.

Whole genome sequencing (WGS) is now playing a central role in the control and study of tuberculosis. Factors hindering genomic data interpretation can difficult our understanding of the pathogen biology and lead to incorrect clinical predictions. Analyzing an extensive Mycobacterium tuberculosis dataset comprising more than 1,500 sequencing samples from different published works, with additional samples sequenced in our laboratory, we find that contamination with non-target DNA is a common phenomenon among WGS studies. By using this data and in-silico simulations, we show that even subtle contaminations can produce dozens of false variants and large miscalculations of allele frequencies, often leading to errors that are very hard to detect and propagate through the analysis. In our dataset, 94% of the polymorphic positions were incorrectly identified due to contaminations. We exemplify the consequences of these errors in the context of clinical predictions for all the studies analyzed, and demonstrate that unexpected contaminations suppose a major pitfall in WGS studies. In addition, we present an approach based on the removal of contaminant reads that shows an outstanding performance analyzing both clean and contaminated data. Based in our findings, applicable to most of organisms, we urge for the implementation of contamination-aware analysis pipelines.

We report the development of the ReproGenomics Viewer (RGV), a multi- and cross-species working environment for the visualization, mining and comparison of published omics data sets for the reproductive science community. The system currently embeds 15 published data sets related to gametogenesis from nine model organisms. Data sets have been curated and conveniently organized into broad categories including biological topics, technologies, species and publications. RGV's modular design for both organisms and genomic tools enables users to upload and compare their data with that from the data sets embedded in the system in a cross-species manner. The RGV is freely available at http://rgv.genouest.org.

We report full-length draft de novo genome assemblies for 16 widely used inbred mouse strains and find extensive strain-specific haplotype variation. We identify and characterize 2,567 regions on the current mouse reference genome exhibiting the greatest sequence diversity. These regions are enriched for genes involved in pathogen defence and immunity and exhibit enrichment of transposable elements and signatures of recent retrotransposition events. Combinations of alleles and genes unique to an individual strain are commonly observed at these loci, reflecting distinct strain phenotypes. We used these genomes to improve the mouse reference genome, resulting in the completion of 10 new gene structures. Also, 62 new coding loci were added to the reference genome annotation. These genomes identified a large, previously unannotated, gene (Efcab3-like) encoding 5,874 amino acids. Mutant Efcab3-like mice display anomalies in multiple brain regions, suggesting a possible role for this gene in the regulation of brain development.Sequence assemblies for the genomes of 16 widely used inbred laboratory mouse strains highlight considerable strain-specific haplotype variation and allow for the identification of regions with the greatest sequence diversity between strains.

We introduce CUT&RUNTools as a flexible, general pipeline for facilitating the identification of chromatin-associated protein binding and genomic footprinting analysis from antibody-targeted CUT&RUN primary cleavage data. CUT&RUNTools extracts endonuclease cut site information from sequences of short-read fragments and produces single-locus binding estimates, aggregate motif footprints, and informative visualizations to support the high-resolution mapping capability of CUT&RUN. CUT&RUNTools is available at https://bitbucket.org/qzhudfci/cutruntools/.

We introduce CUT&RUNTools as a flexible, general pipeline for facilitating the identification of chromatin-associated protein binding and genomic footprinting analysis from antibody-targeted CUT&RUN primary cleavage data. CUT&RUNTools extracts endonuclease cut site information from sequences of short read fragments and produces single-locus binding estimates, aggregate motif footprints, and informative visualizations to support the high-resolution mapping capability of CUT&RUN. CUT&RUNTools is available at https://bitbucket.org/qzhudfci/cutruntools/.

We here pioneer a low-cost assembly strategy for 20 Heliconiini genomes to characterize the evolutionary history of the rapidly radiating genus Heliconius. A bifurcating tree provides a poor fit to the data, and we therefore explore a reticulate phylogeny for Heliconius. We probe the genomic architecture of gene flow, and develop a new method to distinguish incomplete lineage sorting from introgression. We find that most loci with non-canonical histories arose through introgression, and are strongly underrepresented in regions of low recombination and low gene density. This is expected if introgressed alleles are more likely to be purged in such regions due to tighter linkage with incompatibility loci. Finally, we identify a hitherto unrecognized inversion, and show it is a convergent structural rearrangement that captures a known color pattern switch locus within the genus. Our multi-genome assembly approach enables an improved understanding of adaptive radiation.

We evaluated the performance of bcbio-nextgen pipeline for clinical whole exome sequencing (WES) using Genome in a Bottle (GIAB) NA12878 reference sample with Agilent Clinical Research Exome (CRE) v2 panel and Agilent SureSelect XT2 library preparation kit. The current literature on NA12878 reference materials predominantly describes the performance of WES-based pipelines for Illumina HiSeqX. Here, we evaluated the analytical performance of the bcbio-nextgen pipeline for WES for Illumina NextSeq500. We compared the performance of the bcbio-nextgen pipeline with the performance of other WES-based pipelines reported in literature. Our comparative analyses showed that the sensitivity and precision of our study compared favourably with the performance of other bioinformatics pipelines. Notably, our results demonstrate that analytical sensitivity of WES deteriorates outside of capture regions and suggests that any variants identified outside of capture regions are more likely to be lower confidence variants and therefore require confirmatory testing. Careful consideration of intersect between genomic regions of interest (ROI) and capture regions is required to achieve optimal performance of WES for clinical applications.

Typical analysis workflows map reads to a reference genome in order to detect genetic variants. Generating such alignments introduces references biases, in particular against insertion alleles absent in the reference and comes with substantial computational burden. In contrast, recent k-mer-based genotyping methods are fast, but struggle in repetitive or duplicated regions of the genome. We propose a novel algorithm, called PanGenie, that leverages a pangenome reference built from haplotype-resolved genome assemblies in conjunction with k-mer count information from raw, short-read sequencing data to genotype a wide spectrum of genetic variation. The given haplotypes enable our method to take advantage of linkage information to aid genotyping in regions poorly covered by unique k-mers and provides access to regions otherwise inaccessible by short reads. Compared to classic mapping-based approaches, our approach is more than 4× faster at 30× coverage and at the same time, reached significantly better genotype concordances for almost all variant types and coverages tested. Improvements are especially pronounced for large insertions (> 50bp), where we are able to genotype > 99.9% of all tested variants with over 90% accuracy at 30× short-read coverage, where the best competing tools either typed less than 60% of variants or reached accuracies below 70%. PanGenie now enables the inclusion of this commonly neglected variant type in downstream analyses.

Tools for mapping transcription factor binding sites are evaluated with a protocol to map gene regulatory from open chromatin data sets. Determining where transcription factors (TFs) bind in genomes provides insight into which transcriptional programs are active across organs, tissue types, and environmental conditions. Recent advances in high-throughput profiling of regulatory DNA have yielded large amounts of information about chromatin accessibility. Interpreting the functional significance of these data sets requires knowledge of which regulators are likely to bind these regions. This can be achieved by using information about TF-binding preferences, or motifs, to identify TF-binding events that are likely to be functional. Although different approaches exist to map motifs to DNA sequences, a systematic evaluation of these tools in plants is missing. Here, we compare four motif-mapping tools widely used in the Arabidopsis (Arabidopsis thaliana) research community and evaluate their performance using chromatin immunoprecipitation data sets for 40 TFs. Downstream gene regulatory network (GRN) reconstruction was found to be sensitive to the motif mapper used. We further show that the low recall of Find Individual Motif Occurrences, one of the most frequently used motif-mapping tools, can be overcome by using an Ensemble approach, which combines results from different mapping tools. Several examples are provided demonstrating how the Ensemble approach extends our view on transcriptional control for TFs active in different biological processes. Finally, a protocol is presented to effectively derive more complete cell type-specific GRNs through the integrative analysis of open chromatin regions, known binding site information, and expression data sets. This approach will pave the way to increase our understanding of GRNs in different cellular conditions.

Three-dimensional (3D) conformation of the chromatin is crucial to stringently regulate gene expression patterns and DNA replication in a cell-type specific manner. HiC is a key technique for measuring 3D chromatin interactions genome wide. Estimating and predicting the resolution of a library is an essential step in any HiC experimental design. Here, we present the mathematical concepts to estimate the resolution of a library and predict whether deeper sequencing would enhance the resolution. We have developed HiCRes, a docker pipeline, by applying these concepts to human and mouse HiC libraries.

Thousands of alternative exons are spliced out of messenger RNA to increase protein diversity. High‐throughput sequencing of short cDNA fragments (RNA‐seq) generates a genome‐wide snapshot of these post‐transcriptional processes. RNA‐seq reads yield insights into the regulation of alternative splicing by revealing the usage of known or unknown splice sites as well as the expression level of exons. Constitutive exons are never covered by split alignments, whereas alternative exonic parts are located within highly expressed splicing junctions. The ratio between reads including or excluding exons, also known as percent spliced in index (PSI), indicates how efficiently sequences of interest are spliced into transcripts. This protocol describes a method to calculate the PSI without prior knowledge of splicing patterns. It provides a quantitative, global assessment of exon usage that can be integrated with other tools that identify differential isoform processing. Novel, complex splicing events along a genetic locus can be visualized in an exon‐centric manner and compared across conditions. © 2015 by John Wiley & Sons, Inc.

This report describes the miRQuest - a novel middleware available in a Web server that allows the end user to do the miRNA research in a user-friendly way. It is known that there are many prediction tools for microRNA (miRNA) identification that use different programming languages and methods to realize this task. It is difficult to understand each tool and apply it to diverse datasets and organisms available for miRNA analysis. miRQuest can easily be used by biologists and researchers with limited experience with bioinformatics. We built it using the middleware architecture on a Web platform for miRNA research that performs two main functions: i) integration of different miRNA prediction tools for miRNA identification in a user-friendly environment; and ii) comparison of these prediction tools. In both cases, the user provides sequences (in FASTA format) as an input set for the analysis and comparisons. All the tools were selected on the basis of a survey of the literature on the available tools for miRNA prediction. As results, three different cases of use of the tools are also described, where one is the miRNA identification analysis in 30 different species. Finally, miRQuest seems to be a novel and useful tool; and it is freely available for both benchmarking and miRNA identification at http://mirquest.integrativebioinformatics.me/.

These findings demonstrate that manipulation of T-cell EZH2 in cellular therapies may yield cellular products able to withstand solid tumor metabolic–deficient environments. T-cell exhaustion in cancer is linked to poor clinical outcomes, where evidence suggests T-cell metabolic changes precede functional exhaustion. Direct competition between tumor-infiltrating lymphocytes (TIL) and cancer cells for metabolic resources often renders T cells dysfunctional. Environmental stress produces epigenome remodeling events within TIL resulting from loss of the histone methyltransferase EZH2. Here, we report an epigenetic mechanism contributing to the development of metabolic exhaustion in TIL. A multiomics approach revealed a Cdkn2a.Arf-mediated, p53-independent mechanism by which EZH2 inhibition leads to mitochondrial dysfunction and the resultant exhaustion. Reprogramming T cells to express a gain-of-function EZH2 mutant resulted in an enhanced ability of T cells to inhibit tumor growth in vitro and in vivo. Our data suggest that manipulation of T-cell EZH2 within the context of cellular therapies may yield lymphocytes that are able to withstand harsh tumor metabolic environments and collateral pharmacologic insults. Significance: These findings demonstrate that manipulation of T-cell EZH2 in cellular therapies may yield cellular products able to withstand solid tumor metabolic–deficient environments. Graphical Abstract: http://cancerres.aacrjournals.org/content/canres/80/21/4707/F1.large.jpg. Graphical Abstract

The telomere to telomere (T2T) genome project discovered and mapped ∼240 million additional base pairs of primarily telomeric and centromeric reads. Much of this sequence was comprised of satellite sequences and large segmental duplications. We evaluated the extent to which human bona fide microRNAs (miRNAs) may be found in additional paralogous genomic loci or if previously undescribed microRNAs are present in these newly sequenced regions of the human genome. New genomic regions of the T2T project spanning ∼240 million bp of sequence were obtained and evaluated by blastn for the human miRNAs contained in MirGeneDB2.0 (N=556) and miRBase (N = 1917) along with all species of MirGeneDB2.0 miRNAs (N=10,899). Additionally, bowtie was used to compare unmapped reads from >4,000 primary cell samples to the new T2T sequence. Based on sequence and structure, no bona fide miRNAs were identified. Ninety-seven miRNAs of questionable authenticity (frequently known repeat elements) were identified from the miRBase dataset across the newly described regions of the human genome. These 97 represent only 51 miRNA families due to paralogy of highly similar miRNAs such as 24 members of the hsa-mir-548 family. Altogether, this data strongly supports our having identified widely expressed bona fide miRNAs in the human genome and move us further toward the completion of human miRNA discovery.

The success of clinical genomics using next generation sequencing (NGS) requires the accurate and consistent identification of personal genome variants. Assorted variant calling methods have been developed, which show low concordance between their calls. Hence, a systematic comparison of the variant callers could give important guidance to NGS-based clinical genomics. Recently, a set of high-confident variant calls for one individual (NA12878) has been published by the Genome in a Bottle (GIAB) consortium, enabling performance benchmarking of different variant calling pipelines. Based on the gold standard reference variant calls from GIAB, we compared the performance of thirteen variant calling pipelines, testing combinations of three read aligners—BWA-MEM, Bowtie2, and Novoalign—and four variant callers—Genome Analysis Tool Kit HaplotypeCaller (GATK-HC), Samtools mpileup, Freebayes and Ion Proton Variant Caller (TVC), for twelve data sets for the NA12878 genome sequenced by different platforms including Illumina2000, Illumina2500, and Ion Proton, with various exome capture systems and exome coverage. We observed different biases toward specific types of SNP genotyping errors by the different variant callers. The results of our study provide useful guidelines for reliable variant identification from deep sequencing of personal genomes.

The rapidly declining costs of sequencing human genomes and exomes are providing deeper insights into genomic variation than previously possible. Growing sequence datasets are uncovering large numbers of rare single-nucleotide variants (SNVs) in coding regions, many of which may even be unique to single individuals. The rarity of such variants makes it difficult to use conventional variant-phenotype associations as a means of predicting their potential impacts. As such, protein structures may help to provide the needed means for inferring otherwise difficult-to-discern rare SNV-phenotype associations. Previous efforts have sought to quantify the effects of SNVs on structures by evaluating their impacts on global stability. However, local perturbations can severely impact functionality (such as catalysis, allosteric regulation, interactions and specificity) without strongly disrupting global stability. Here, we describe a workflow in which localized frustration (which quantifies unfavorable residue-residue interactions) is employed as a metric to investigate such effects. We apply frustration to study the impacts of a large number of SNVs available throughout a number of next-generation sequencing datasets. Most of our observations are intuitively consistent: we observe that disease-associated SNVs have a strong proclivity to induce strong changes in localized frustration, and rare variants tend to disrupt local interactions to a larger extent than do common variants. Furthermore, we observe that somatic SNVs associated with oncogenes induce stronger perturbations at the surface, whereas those associated with tumor suppressor genes (TSGs) induce stronger perturbations in the interior. These findings are consistent with the notion that gain-of-function (for oncogenes) and loss-of-function events (for TSGs) may act through changes in regulatory interactions and basic functionality, respectively

The quality of RNA sequencing data relies on specific priming by the primer used for reverse transcription (RT-primer). Non-specific annealing of the RT-primer to the RNA template can generate reads with incorrect cDNA ends and can cause misinterpretation of data (RT mispriming). This kind of artifact in RNA-seq based technologies is underappreciated and currently no adequate tools exist to computationally remove them from published datasets. We show that mispriming can occur with as little as 2 bases of complementarity at the 3’ end of the primer followed by intermittent regions of complementarity. We also provide a computational pipeline that identifies cDNA reads produced from RT mispriming, allowing users to filter them out from any aligned dataset. Using this analysis pipeline, we identify thousands of mispriming events in a dozen published datasets from diverse technologies including short RNA-seq, total/mRNA-seq, HITS-CLIP and GRO-seq. We further show how RT-mispriming can lead to misinterpretation of data. In addition to providing a solution to computationally remove RT-misprimed reads, we also propose an experimental solution to avoid RT-mispriming by performing RNA-seq using thermostable group II intron derived reverse transcriptase (TGIRT-seq).

The origin of Tibetans remains one of the most contentious puzzles in history, anthropology, and genetics. Analyses of deeply sequenced (30×-60×) genomes of 38 Tibetan highlanders and 39 Han Chinese lowlanders, together with available data on archaic and modern humans, allow us to comprehensively characterize the ancestral makeup of Tibetans and uncover their origins. Non-modern human sequences compose ∼6% of the Tibetan gene pool and form unique haplotypes in some genomic regions, where Denisovan-like, Neanderthal-like, ancient-Siberian-like, and unknown ancestries are entangled and elevated. The shared ancestry of Tibetan-enriched sequences dates back to ∼62,000-38,000 years ago, predating the Last Glacial Maximum (LGM) and representing early colonization of the plateau. Nonetheless, most of the Tibetan gene pool is of modern human origin and diverged from that of Han Chinese ∼15,000 to ∼9,000 years ago, which can be largely attributed to post-LGM arrivals. Analysis of ∼200 contemporary populations showed that Tibetans share ancestry with populations from East Asia (∼82%), Central Asia and Siberia (∼11%), South Asia (∼6%), and western Eurasia and Oceania (∼1%). Our results support that Tibetans arose from a mixture of multiple ancestral gene pools but that their origins are much more complicated and ancient than previously suspected. We provide compelling evidence of the co-existence of Paleolithic and Neolithic ancestries in the Tibetan gene pool, indicating a genetic continuity between pre-historical highland-foragers and present-day Tibetans. In particular, highly differentiated sequences harbored in highlanders' genomes were most likely inherited from pre-LGM settlers of multiple ancestral origins (SUNDer) and maintained in high frequency by natural selection.

The major DNA sequencing technologies in use today produce either highly-accurate short reads or noisy long reads. We developed a protocol based on single-molecule, circular consensus sequencing (CCS) to generate highly-accurate (99.8%) long reads averaging 13.5 kb and applied it to sequence the well-characterized human HG002/NA24385. We optimized existing tools to comprehensively detect variants, achieving precision and recall above 99.91% for SNVs, 95.98% for indels, and 95.99% for structural variants. We estimate that 2,434 discordances are correctable mistakes in the high-quality Genome in a Bottle benchmark. Nearly all (99.64%) variants are phased into haplotypes, which further improves variant detection. De novo assembly produces a highly contiguous and accurate genome with contig N50 above 15 Mb and concordance of 99.998%. CCS reads match short reads for small variant detection, while enabling structural variant detection and de novo assembly at similar contiguity and markedly higher concordance than noisy long reads.