Sustainable data analysis with Snakemake

Accelerating genomic workflows using NVIDIA Parabricks

BACKGROUND

As genome sequencing becomes better integrated into scientific research, government policy, and personalized medicine, the primary challenge for researchers is shifting from generating raw data to analyzing these vast datasets. Although much work has been done to reduce compute times using various configurations of traditional CPU computing infrastructures, Graphics Processing Units (GPUs) offer opportunities to accelerate genomic workflows by orders of magnitude. Here we benchmark one GPU-accelerated software suite called NVIDIA Parabricks on Amazon Web Services (AWS), Google Cloud Platform (GCP), and an NVIDIA DGX cluster. We benchmarked six variant calling pipelines, including two germline callers (HaplotypeCaller and DeepVariant) and four somatic callers (Mutect2, Muse, LoFreq, SomaticSniper).

RESULTS

We achieved up to 65 × acceleration with germline variant callers, bringing HaplotypeCaller runtimes down from 36 h to 33 min on AWS, 35 min on GCP, and 24 min on the NVIDIA DGX. Somatic callers exhibited more variation between the number of GPUs and computing platforms. On cloud platforms, GPU-accelerated germline callers resulted in cost savings compared with CPU runs, whereas some somatic callers were more expensive than CPU runs because their GPU acceleration was not sufficient to overcome the increased GPU cost.

CONCLUSIONS

Germline variant callers scaled well with the number of GPUs across platforms, whereas somatic variant callers exhibited more variation in the number of GPUs with the fastest runtimes, suggesting that, at least with the version of Parabricks used here, these workflows are less GPU optimized and require benchmarking on the platform of choice before being deployed at production scales. Our study demonstrates that GPUs can be used to greatly accelerate genomic workflows, thus bringing closer to grasp urgent societal advances in the areas of biosurveillance and personalized medicine.

Implementing landscape genetics in molecular epidemiology to determine drivers of vector-borne disease: A malaria case study

This study employs landscape genetics to investigate the environmental drivers of a deadly vector-borne disease, malaria caused by Plasmodium falciparum, in a more spatially comprehensive manner than any previous work. With 1804 samples from 44 sites collected in western Kenya in 2012 and 2013, we performed resistance surface analysis to show that Lake Victoria acts as a barrier to transmission between areas north and south of the Winam Gulf. In addition, Mantel correlograms clearly showed significant correlations between genetic and geographic distance over short distances (less than 70 km). In both cases, we used an identity-by-state measure of relatedness tailored to find highly related individual parasites in order to focus on recent gene flow that is more relevant to disease transmission. To supplement these results, we performed conventional population genetics analyses, including Bayesian clustering methods and spatial ordination techniques. These analyses revealed some differentiation on the basis of geography and elevation and a cluster of genetic similarity in the lowlands north of the Winam Gulf of Lake Victoria. Taken as a whole, these results indicate low overall genetic differentiation in the Lake Victoria region, but with some separation of parasite populations north and south of the Winam Gulf that is explained by the presence of the lake as a geographic barrier to gene flow. We recommend similar landscape genetics analyses in future molecular epidemiology studies of vector-borne diseases to extend and contextualize the results of traditional population genetics.

Lineage tracing reveals fate bias and transcriptional memory in human B cells

We combined single-cell transcriptomics and lineage tracing to understand fate choice in human B cells. Using the antibody sequences of B cells, we tracked clones during in vitro differentiation. Clonal analysis revealed a subset of IgM+ B cells which were more proliferative than other B-cell types. Whereas the population of B cells adopted diverse states during differentiation, clones had a restricted set of fates available to them; there were two times more single-fate clones than expected given population-level cell-type diversity. This implicated a molecular memory of initial cell states that was propagated through differentiation. We then identified the genes which had strongest coherence within clones. These genes significantly overlapped known B-cell fate determination programs, suggesting the genes which determine cell identity are most robustly controlled on a clonal level. Persistent clonal identities were also observed in human plasma cells from bone marrow, indicating that these transcriptional programs maintain long-term cell identities in vivo. Thus, we show how cell-intrinsic fate bias influences differentiation outcomes in vitro and in vivo.

Metagenomic Analysis of the Abundance and Composition of Antibiotic Resistance Genes in Hospital Wastewater in Benin, Burkina Faso, and Finland

Antibiotic resistance is a global threat to human health, with the most severe effect in low- and middle-income countries. We explored the presence of antibiotic resistance genes (ARGs) in the hospital wastewater (HWW) of nine hospitals in Benin and Burkina Faso, two low-income countries in West Africa, with shotgun metagenomic sequencing. For comparison, we also studied six hospitals in Finland. The highest sum of the relative abundance of ARGs in the 68 HWW samples was detected in Benin and the lowest in Finland. HWW resistomes and mobilomes in Benin and Burkina Faso resembled each other more than those in Finland. Many carbapenemase genes were detected at various abundances, especially in HWW from Burkina Faso and Finland. The bla_GES genes, the most widespread carbapenemase gene in the Beninese HWW, were also found in water intended for hand washing and in a puddle at a hospital yard in Benin. mcr genes were detected in the HWW of all three countries, with mcr-5 being the most common mcr gene. These and other mcr genes were observed in very high relative abundances, even in treated wastewater in Burkina Faso and a street gutter in Benin. The results highlight the importance of wastewater treatment, with particular attention to HWW. IMPORTANCE The global emergence and increased spread of antibiotic resistance threaten the effectiveness of antibiotics and, thus, the health of the entire population. Therefore, understanding the resistomes in different geographical locations is crucial in the global fight against the antibiotic resistance crisis. However, this information is scarce in many low- and middle-income countries (LMICs), such as those in West Africa. In this study, we describe the resistomes of hospital wastewater in Benin and Burkina Faso and, as a comparison, Finland. Our results help to understand the hitherto unrevealed resistance in Beninese and Burkinabe hospitals. Furthermore, the results emphasize the importance of wastewater management infrastructure design to minimize exposure events between humans, HWW, and the environment, preventing the circulation of resistant bacteria and ARGs between humans (hospitals and community) and the environment.

The omnitig framework can improve genome assembly contiguity in practice

Despite the long history of genome assembly research, there remains a large gap between the theoretical and practical work. There is practical software with little theoretical underpinning of accuracy on one hand and theoretical algorithms which have not been adopted in practice on the other. In this paper we attempt to bridge the gap between theory and practice by showing how the theoretical safe-and-complete framework can be integrated into existing assemblers in order to improve contiguity. The optimal algorithm in this framework, called the omnitig algorithm, has not been used in practice due to its complexity and its lack of robustness to real data. Instead, we pursue a simplified notion of omnitigs, giving an efficient algorithm to compute them and demonstrating their safety under certain conditions. We modify two assemblers (wtdbg2 and Flye) by replacing their unitig algorithm with the simple omnitig algorithm. We test our modifications using real HiFi data from the Drosophilia melanogaster and the Caenorhabditis elegans genome. Our modified algorithms lead to a substantial improvement in alignment-based contiguity, with negligible computational costs and either no or a small increase in the number of misassemblies.

PlantTribes2: Tools for comparative gene family analysis in plant genomics

Plant genome-scale resources are being generated at an increasing rate as sequencing technologies continue to improve and raw data costs continue to fall; however, the cost of downstream analyses remains large. This has resulted in a considerable range of genome assembly and annotation qualities across plant genomes due to their varying sizes, complexity, and the technology used for the assembly and annotation. To effectively work across genomes, researchers increasingly rely on comparative genomic approaches that integrate across plant community resources and data types. Such efforts have aided the genome annotation process and yielded novel insights into the evolutionary history of genomes and gene families, including complex non-model organisms. The essential tools to achieve these insights rely on gene family analysis at a genome-scale, but they are not well integrated for rapid analysis of new data, and the learning curve can be steep. Here we present PlantTribes2, a scalable, easily accessible, highly customizable, and broadly applicable gene family analysis framework with multiple entry points including user provided data. It uses objective classifications of annotated protein sequences from existing, high-quality plant genomes for comparative and evolutionary studies. PlantTribes2 can improve transcript models and then sort them, either genome-scale annotations or individual gene coding sequences, into pre-computed orthologous gene family clusters with rich functional annotation information. Then, for gene families of interest, PlantTribes2 performs downstream analyses and customizable visualizations including, (1) multiple sequence alignment, (2) gene family phylogeny, (3) estimation of synonymous and non-synonymous substitution rates among homologous sequences, and (4) inference of large-scale duplication events. We give examples of PlantTribes2 applications in functional genomic studies of economically important plant families, namely transcriptomics in the weedy Orobanchaceae and a core orthogroup analysis (CROG) in Rosaceae. PlantTribes2 is freely available for use within the main public Galaxy instance and can be downloaded from GitHub or Bioconda. Importantly, PlantTribes2 can be readily adapted for use with genomic and transcriptomic data from any kind of organism.

Mis-spliced transcripts generate de novo proteins in TDP-43-related ALS/FTD

Functional loss of TDP-43, an RNA-binding protein genetically and pathologically linked to ALS and FTD, leads to inclusion of cryptic exons in hundreds of transcripts during disease. Cryptic exons can promote degradation of affected transcripts, deleteriously altering cellular function through loss-of-function mechanisms. However, the possibility of de novo protein synthesis from cryptic exon transcripts has not been explored. Here, we show that mRNA transcripts harboring cryptic exons generate de novo proteins both in TDP-43 deficient cellular models and in disease. Using coordinated transcriptomic and proteomic studies of TDP-43 depleted iPSC-derived neurons, we identified numerous peptides that mapped to cryptic exons. Cryptic exons identified in iPSC models were highly predictive of cryptic exons expressed in brains of patients with TDP-43 proteinopathy, including cryptic transcripts that generated de novo proteins. We discovered that inclusion of cryptic peptide sequences in proteins altered their interactions with other proteins, thereby likely altering their function. Finally, we showed that these de novo peptides were present in CSF from patients with ALS. The demonstration of cryptic exon translation suggests new mechanisms for ALS pathophysiology downstream of TDP-43 dysfunction and may provide a strategy for novel biomarker development.

Integrative proteogenomics using ProteomeGenerator2

Recent advances in nucleic acid sequencing now permit rapid and genome-scale analysis of genetic variation and transcription, enabling population-scale studies of human biology, disease, and diverse organisms. Likewise, advances in mass spectrometry proteomics now permit highly sensitive and accurate studies of protein expression at the whole proteome-scale. However, most proteomic studies rely on consensus databases to match spectra to peptide and proteins sequences, and thus remain limited to the analysis of canonical protein sequences. Here, we develop ProteomeGenerator2 (PG2), based on the scalable and modular ProteomeGenerator framework. PG2 integrates genome and transcriptome sequencing to incorporate protein variants containing amino acid substitutions, insertions, and deletions, as well as non-canonical reading frames, exons, and other variants caused by genomic and transcriptomic variation. We benchmarked PG2 using synthetic data and genomic, transcriptomic, and proteomic analysis of human leukemia cells. PG2 can be integrated with current and emerging sequencing technologies, assemblers, variant callers, and mass spectral analysis algorithms, and is available open-source from https://github.com/kentsisresearchgroup/ProteomeGenerator2 .

Genomic selection for survival under naturally occurring Saprolegnia oomycete infection in farmed European whitefish Coregonus lavaretus

Saprolegnia oomycete infection causes serious economic losses and reduces fish health in aquaculture. Genomic selection based on thousands of DNA markers is a powerful tool to improve fish traits in selective breeding programs. Our goal was to develop a single nucleotide polymorphism (SNP) marker panel and to test its use in genomic selection for improved survival against Saprolegnia infection in European whitefish Coregonus lavaretus, the second most important farmed fish species in Finland. We used a double digest restriction site associated DNA (ddRAD) genotyping by sequencing method to produce a SNP panel, and we tested it analyzing data from a cohort of 1,335 fish, which were measured at different times for mortality to Saprolegnia oomycete infection and weight traits. We calculated the genetic relationship matrix (GRM) from the genome-wide genetic data, integrating it in multivariate mixed models used for the estimation of variance components and genomic breeding values (GEBVs), and to carry out Genome-Wide Association Studies for the presence of quantitative trait loci (QTL) affecting the phenotypes in analysis. We identified one major QTL on chromosome 6 affecting mortality to Saprolegnia infection, explaining 7.7% to 51.3% of genetic variance, and a QTL for weight on chromosome 4, explaining 1.8% to 5.4% of genetic variance. Heritability for mortality was 0.20 to 0.43 on the liability scale, and heritability for weight was 0.44 to 0.53. The QTL for mortality showed an additive allelic effect. We tested whether integrating the QTL for mortality as a fixed factor, together with a new GRM calculated excluding the QTL from the genetic data, would improve the accuracy estimation of GEBVs. This test was done through a cross-validation approach, which indicated that the inclusion of the QTL increased the mean accuracy of the GEBVs by 0.28 points, from 0.33 to 0.61, relative to the use of full GRM only. The area under the curve of the receiver-operator curve for mortality increased from 0.58 to 0.67 when the QTL was included in the model. The inclusion of the QTL as a fixed effect in the model increased the correlation between the GEBVs of early mortality with the late mortality, compared to a model that did not include the QTL. These results validate the usability of the produced SNP panel for genomic selection in European whitefish and highlight the opportunity for modeling QTLs in genomic evaluation of mortality due to Saprolegnia infection.

Cellsnake: a user-friendly tool for single-cell RNA sequencing analysis

BACKGROUND

Single-cell RNA sequencing (scRNA-seq) provides high-resolution transcriptome data to understand the heterogeneity of cell populations at the single-cell level. The analysis of scRNA-seq data requires the utilization of numerous computational tools. However, nonexpert users usually experience installation issues, a lack of critical functionality or batch analysis modes, and the steep learning curves of existing pipelines.

RESULTS

We have developed cellsnake, a comprehensive, reproducible, and accessible single-cell data analysis workflow, to overcome these problems. Cellsnake offers advanced features for standard users and facilitates downstream analyses in both R and Python environments. It is also designed for easy integration into existing workflows, allowing for rapid analyses of multiple samples.

CONCLUSION

As an open-source tool, cellsnake is accessible through Bioconda, PyPi, Docker, and GitHub, making it a cost-effective and user-friendly option for researchers. By using cellsnake, researchers can streamline the analysis of scRNA-seq data and gain insights into the complex biology of single cells.

Genetic Drift Shapes the Evolution of a Highly Dynamic Metapopulation

The dynamics of extinction and (re)colonization in habitat patches are characterizing features of dynamic metapopulations, causing them to evolve differently than large, stable populations. The propagule model, which assumes genetic bottlenecks during colonization, posits that newly founded subpopulations have low genetic diversity and are genetically highly differentiated from each other. Immigration may then increase diversity and decrease differentiation between subpopulations. Thus, older and/or less isolated subpopulations are expected to have higher genetic diversity and less genetic differentiation. We tested this theory using whole-genome pool-sequencing to characterize nucleotide diversity and differentiation in 60 subpopulations of a natural metapopulation of the cyclical parthenogen Daphnia magna. For comparison, we characterized diversity in a single, large, and stable D. magna population. We found reduced (synonymous) genomic diversity, a proxy for effective population size, weak purifying selection, and low rates of adaptive evolution in the metapopulation compared with the large, stable population. These differences suggest that genetic bottlenecks during colonization reduce effective population sizes, which leads to strong genetic drift and reduced selection efficacy in the metapopulation. Consistent with the propagule model, we found lower diversity and increased differentiation in younger and also in more isolated subpopulations. Our study sheds light on the genomic consequences of extinction-(re)colonization dynamics to an unprecedented degree, giving strong support for the propagule model. We demonstrate that the metapopulation evolves differently from a large, stable population and that evolution is largely driven by genetic drift.

A small bacteriophage protein determines the hierarchy over co-residential jumbo phage in Bacillus thuringiensis serovar israelensis

Bacillus thuringiensis serovar israelensis is the most widely used biopesticide against insects, including vectors of animal and human diseases. Among several extrachromosomal elements, this endospore-forming entomopathogen harbors two bacteriophages: a linear DNA replicon named GIL01 that does not integrate into the chromosome during lysogeny and a circular-jumbo prophage known as pBtic235. Here, we show that GIL01 hinders the induction of cohabiting prophage pBtic235. The GIL01-encoded small protein, gp7, which interacts with the host LexA repressor, is a global transcription regulator and represses the induction of pBtic235 after DNA damage to presumably allow GIL01 to multiply first. In a complex with host LexA in stressed cells, gp7 down-regulates the expression of more than 250 host and pBtic235 genes, many of which are involved in the cellular functions of genome maintenance, cell-wall transport, and membrane and protein stability. We show that gp7 homologs that are found exclusively in bacteriophages act in a similar fashion to enhance LexA's binding to DNA, while likely also affecting host gene expression. Our results provide evidence that GIL01 influences both its host and its co-resident bacteriophage.

Plastaumatic: Automating plastome assembly and annotation

Plastome sequence data is most often extracted from plant whole genome sequencing data and need to be assembled and annotated separately from the nuclear genome sequence. In projects comprising multiple genomes, it is labour intense to individually process the plastomes as it requires many steps and software. This study developed Plastaumatic - an automated pipeline for both assembly and annotation of plastomes, with the scope of the researcher being able to load whole genome sequence data with minimal manual input, and therefore a faster runtime. The main structure of the current automated pipeline includes trimming of adaptor and low-quality sequences using fastp, de novo plastome assembly using NOVOPlasty, standardization and quality checking of the assembled genomes through a custom script utilizing BLAST+ and SAMtools, annotation of the assembled genomes using AnnoPlast, and finally generating the required files for NCBI GenBank submissions. The pipeline is demonstrated with 12 potato accessions and three soybean accessions.

A pseudovirus system enables deep mutational scanning of the full SARS-CoV-2 spike

A major challenge in understanding SARS-CoV-2 evolution is interpreting the antigenic and functional effects of emerging mutations in the viral spike protein. Here we describe a new deep mutational scanning platform based on non-replicative pseudotyped lentiviruses that directly quantifies how large numbers of spike mutations impact antibody neutralization and pseudovirus infection. We demonstrate this new platform by making libraries of the Omicron BA.1 and Delta spikes. These libraries each contain ~7000 distinct amino-acid mutations in the context of up to ~135,000 unique mutation combinations. We use these libraries to map escape mutations from neutralizing antibodies targeting the receptor binding domain, N-terminal domain, and S2 subunit of spike. Overall, this work establishes a high-throughput and safe approach to measure how ~10 5 combinations of mutations affect antibody neutralization and spike-mediated infection. Notably, the platform described here can be extended to the entry proteins of many other viruses.

Emergence of Tn1999.7, a New Transposon in blaOXA-48-Harboring Plasmids Associated with Increased Plasmid Stability

OXA-48 is the most common carbapenemase in Enterobacterales in Germany and many other European countries. Depending on the genomic location of bla_OXA-48, OXA-48-producing isolates vary in phenotype and intra- and interspecies transferability of bla_OXA-48. In most bacterial isolates, bla_OXA-48 is located on one of seven variants of Tn1999 (Tn1999.1 to Tn1999.6 and invTn1999.2). Here, a novel Tn1999 variant, Tn1999.7, is described, which was identified in 11 clinical isolates from 2016 to 2020. Tn1999.7 differs from Tn1999.1 by the insertion of the 8,349-bp Tn3 family transposon Tn7442 between the lysR gene and bla_OXA-48 open reading frame. Tn7442 carries genes coding for a restriction endonuclease and a DNA methyltransferase as cargo, forming a type III restriction modification system. Tn1999.7 was carried on an ~71-kb IncL plasmid in 9/11 isolates. In one isolate, Tn1999.7 was situated on an ~76-kb plasmid, harboring an additional insertion sequence in the plasmid backbone. In one isolate, the plasmid size is only ~63 kb due to a deletion adjacent to Tn7442 that extends into the plasmid backbone. Mean conjugation rates of the Tn1999.7-harboring plasmids in J53 ranged from 4.47 × 10^-5 to 2.03 × 10^-2, similar to conjugation rates of other pOXA-48-type IncL plasmids. The stability of plasmids with Tn1999.7 was significantly higher than that of a Tn1999.2-harboring plasmid in vitro. This increase in stability could be related to the insertion of a restriction-modification system, which can promote postsegregational killing. The increased plasmid stability associated with Tn1999.7 could contribute to the further spread of OXA-48.

Global and precise identification of functional miRNA targets in mESCs by integrative analysis

MicroRNA (miRNA) loaded Argonaute (AGO) complexes regulate gene expression via direct base pairing with their mRNA targets. Previous works suggest that up to 60% of mammalian transcripts might be subject to miRNA‐mediated regulation, but it remains largely unknown which fraction of these interactions are functional in a specific cellular context. Here, we integrate transcriptome data from a set of miRNA‐depleted mouse embryonic stem cell (mESC) lines with published miRNA interaction predictions and AGO‐binding profiles. Using this integrative approach, combined with molecular validation data, we present evidence that < 10% of expressed genes are functionally and directly regulated by miRNAs in mESCs. In addition, analyses of the stem cell‐specific miR‐290‐295 cluster target genes identify TFAP4 as an important transcription factor for early development. The extensive datasets developed in this study will support the development of improved predictive models for miRNA‐mRNA functional interactions.

Detection and Validation of Circular DNA Fragments Using Nanopore Sequencing

Occurrence of extra-chromosomal circular DNA is a phenomenon frequently observed in tumor cells, and the presence of such DNA has been recognized as a marker of adverse outcome across cancer types. We here describe a computational workflow for identification of DNA circles from long-read sequencing data. The workflow is implemented based on the Snakemake workflow management system. Its key step uses a graph-theoretic approach to identify putative circular fragments validated on simulated reads. We then demonstrate robustness of our approach using nanopore sequencing of selectively enriched circular DNA by highly sensitive and specific recovery of plasmids and the mitochondrial genome, which is the only circular DNA in normal human cells. Finally, we show that the workflow facilitates detection of larger circular DNA fragments containing extrachromosomal copies of the MYCN oncogene and the respective breakpoints, which is a potentially useful application in disease monitoring of several cancer types.

A computational workflow for predicting cancer neo-antigens

Neo-antigens presented on cell surface play a pivotal role in the success of immunotherapies. Peptides derived from mutant proteins are thought to be the primary source of neo-antigens presented on the surface of cancer cells. Mutation data from cancer genome sequencing is often used to predict cancer neo-antigens. However, this strategy is associated with significant false positives as many coding mutations may not be expressed at the protein level. Hence, we describe a computational workflow to integrate genomic and proteomic data to predictpotential neo-antigens.

Adding hydrogen atoms to molecular models via fragment superimposition

BACKGROUND

Most experimentally determined structures of biomolecules lack annotated hydrogen positions due to their low electron density. However, thorough structure analysis and simulations require knowledge about the positions of hydrogen atoms. Existing methods for their prediction are either limited to a certain range of molecules or only work effectively on small compounds.

RESULTS

We present a novel algorithm that compiles fragments of molecules with known hydrogen atom positions into a library. Using this library the method is able to predict hydrogen positions for molecules with similar moieties. We show that the method is able to accurately assign hydrogen atoms to most organic compounds including biomacromolecules, if a sufficiently large library is used.

CONCLUSIONS

We bundled the algorithm into the open-source Python package and command line program Hydride. Since usually no additional parametrization is necessary for the problem at hand, the software works out-of-box for a wide range of molecular systems usually within a few seconds of computation time. Hence, we believe that Hydride could be a valuable tool for structural biologists and biophysicists alike.

Perplexity: evaluating transcript abundance estimation in the absence of ground truth

Background

There has been rapid development of probabilistic models and inference methods for transcript abundance estimation from RNA-seq data. These models aim to accurately estimate transcript-level abundances, to account for different biases in the measurement process, and even to assess uncertainty in resulting estimates that can be propagated to subsequent analyses. The assumed accuracy of the estimates inferred by such methods underpin gene expression based analysis routinely carried out in the lab. Although hyperparameter selection is known to affect the distributions of inferred abundances (e.g. producing smooth versus sparse estimates), strategies for performing model selection in experimental data have been addressed informally at best.

Results

We derive perplexity for evaluating abundance estimates on fragment sets directly. We adapt perplexity from the analogous metric used to evaluate language and topic models and extend the metric to carefully account for corner cases unique to RNA-seq. In experimental data, estimates with the best perplexity also best correlate with qPCR measurements. In simulated data, perplexity is well behaved and concordant with genome-wide measurements against ground truth and differential expression analysis. Furthermore, we demonstrate theoretically and experimentally that perplexity can be computed for arbitrary transcript abundance estimation models.

Conclusions

Alongside the derivation and implementation of perplexity for transcript abundance estimation, our study is the first to make possible model selection for transcript abundance estimation on experimental data in the absence of ground truth.

miND (miRNA NGS Discovery pipeline): a small RNA-seq analysis pipeline and report generator for microRNA biomarker discovery studies

In contrast to traditional methods like real-time polymerase chain reaction, next-generation sequencing (NGS), and especially small RNA-seq, enables the untargeted investigation of the whole small RNAome, including microRNAs (miRNAs) but also a multitude of other RNA species. With the promising application of small RNAs as biofluid-based biomarkers, small RNA-seq is the method of choice for an initial discovery study. However, the presentation of specific quality aspects of small RNA-seq data varies significantly between laboratories and is lacking a common (minimal) standard. The miRNA NGS Discovery pipeline (miND) aims to bridge the gap between wet lab scientist and bioinformatics with an easy to setup configuration sheet and an automatically generated comprehensive report that contains all essential qualitative and quantitative results that should be reported. Besides the standard steps like preprocessing, mapping, visualization, and quantification of reads, the pipeline also incorporates differential expression analysis when given the appropriate information regarding sample groups. Although miND has a focus on miRNAs, other RNA species like tRNAs, piRNA, snRNA, or snoRNA are included and mapping statistics are available for further analysis. miND has been developed and tested on a multitude of data sets with various RNA sources (tissue, plasma, extracellular vesicles, urine, etc.) and different species. miND is a Snakemake based pipeline and thus incorporates all advantages using a flexible workflow management system. Reference databases are downloaded, prepared and built with an included (but separate) workflow and thus can easily be updated to the most recent version but also stored for reproducibility. In conclusion, the miND pipeline aims to streamline the bioinformatics processing of small RNA-seq data by standardizing the processing from raw data to a final, comprehensive and reproducible report.

Wei2GO: weighted sequence similarity-based protein function prediction

BACKGROUND

Protein function prediction is an important part of bioinformatics and genomics studies. There are many different predictors available, however most of these are in the form of web-servers instead of open-source locally installable versions. Such local versions are necessary to perform large scale genomics studies due to the presence of limitations imposed by web servers such as queues, prediction speed, and updatability of databases.

METHODS

This paper describes Wei2GO: a weighted sequence similarity and python-based open-source protein function prediction software. It uses DIAMOND and HMMScan sequence alignment searches against the UniProtKB and Pfam databases respectively, transfers Gene Ontology terms from the reference protein to the query protein, and uses a weighing algorithm to calculate a score for the Gene Ontology annotations.

RESULTS

Wei2GO is compared against the Argot2 and Argot2.5 web servers, which use a similar concept, and DeepGOPlus which acts as a reference. Wei2GO shows an increase in performance according to precision and recall curves, Fmax scores, and Smin scores for biological process and molecular function ontologies. Computational time compared to Argot2 and Argot2.5 is decreased from several hours to several minutes.

AVAILABILITY

Wei2GO is written in Python 3, and can be found at https://gitlab.com/mreijnders/Wei2GO.

BUSCO: Assessing Genomic Data Quality and Beyond

Evaluation of the quality of genomic "data products" such as genome assemblies or gene sets is of critical importance in order to recognize possible issues and correct them during the generation of new data. It is equally essential to guide subsequent or comparative analyses with existing data, as the correct interpretation of the results necessarily requires knowledge about the quality level and reliability of the inputs. Using datasets of near universal single-copy orthologs derived from OrthoDB, BUSCO can estimate the completeness and redundancy of genomic data by providing biologically meaningful metrics based on expected gene content. These can complement technical metrics such as contiguity measures (e.g., number of contigs/scaffolds, and N50 values). Here, we describe the use of the BUSCO tool suite to assess different data types that can range from genome assemblies of single isolates and assembled transcriptomes and annotated gene sets to metagenome-assembled genomes where the taxonomic origin of the species is unknown. BUSCO is the only tool capable of assessing all these types of sequences from both eukaryotic and prokaryotic species. The protocols detail the various BUSCO running modes and the novel workflows introduced in versions 4 and 5, including the batch analysis on multiple inputs, the auto-lineage workflow to run assessments without specifying a dataset, and a workflow for the evaluation of (large) eukaryotic genomes. The protocols further cover the BUSCO setup, guidelines to interpret the results, and BUSCO "plugin" workflows for performing common operations in genomics using BUSCO results, such as building phylogenomic trees and visualizing syntenies. © 2021 The Authors. Current Protocols published by Wiley Periodicals LLC. Basic Protocol 1: Assessing an input sequence with a BUSCO dataset specified manually Basic Protocol 2: Assessing an input sequence with a dataset automatically selected by BUSCO Basic Protocol 3: Assessing multiple inputs Alternate Protocol: Decreasing analysis runtime when assessing a large number of small genomes with BUSCO auto-lineage workflow and Snakemake Support Protocol 1: BUSCO setup Support Protocol 2: Visualizing BUSCO results Support Protocol 3: Building phylogenomic trees.

Pseudonatural Products Occur Frequently in Biologically Relevant Compounds

A new methodology for classifying fragment combinations and characterizing pseudonatural products (PNPs) is described. The source code is based on open-source tools and is organized as a Python package. Tasks can be executed individually or within the context of scalable, robust workflows. First, structures are standardized and duplicate entries are filtered out. Then, molecules are probed for the presence of predefined fragments. For molecules with more than one match, fragment combinations are classified. The algorithm considers the pairwise relative position of fragments within the molecule (fused atoms, linkers, intermediary rings), resulting in 18 different possible fragment combination categories. Finally, all combinations for a given molecule are assembled into a fragment combination graph, with fragments as nodes and combination types as edges. This workflow was applied to characterize PNPs in the ChEMBL database via comparison of fragment combination graphs with natural product (NP) references, represented by the Dictionary of Natural Products. The Murcko fragments extracted from 2000 structures previously described were used to define NP fragments. The results indicate that ca. 23% of the biologically relevant compounds listed in ChEMBL comply to the PNP definition and that, therefore, PNPs occur frequently among known biologically relevant small molecules. The majority (>95%) of PNPs contain two to four fragments, mainly (>95%) distributed in five different combination types. These findings may provide guidance for the design of new PNPs.

uPIC-M: Efficient and Scalable Preparation of Clonal Single Mutant Libraries for High-Throughput Protein Biochemistry

New high-throughput biochemistry techniques complement selection-based approaches and provide quantitative kinetic and thermodynamic data for thousands of protein variants in parallel. With these advances, library generation rather than data collection has become rate-limiting. Unlike pooled selection approaches, high-throughput biochemistry requires mutant libraries in which individual sequences are rationally designed, efficiently recovered, sequence-validated, and separated from one another, but current strategies are unable to produce these libraries at the needed scale and specificity at reasonable cost. Here, we present a scalable, rapid, and inexpensive approach for creating User-designed Physically Isolated Clonal-Mutant (uPIC-M) libraries that utilizes recent advances in oligo synthesis, high-throughput sample preparation, and next-generation sequencing. To demonstrate uPIC-M, we created a scanning mutant library of SpAP, a 541 amino acid alkaline phosphatase, and recovered 94% of desired mutants in a single iteration. uPIC-M uses commonly available equipment and freely downloadable custom software and can produce a 5000 mutant library at 1/3 the cost and 1/5 the time of traditional techniques.

Ribo-DT: An automated pipeline for inferring codon dwell times from ribosome profiling data

Protein synthesis is an energy consuming process characterised as a pivotal and highly regulated step in gene expression. The net protein output is dictated by a combination of translation initiation, elongation and termination rates that have remained difficult to measure. Recently, the development of ribosome profiling has enabled the inference of translation parameters through modelling, as this method informs on the ribosome position along the mRNA. Here, we present an automated, reproducible and portable computational pipeline to infer relative single-codon and codon-pair dwell times as well as gene flux from raw ribosome profiling sequencing data. As a case study, we applied our workflow to a publicly available yeast ribosome profiling dataset consisting of 57 independent gene knockouts related to RNA and tRNA modifications. We uncovered the effects of those modifications on translation elongation and codon selection during decoding. In particular, knocking out mod5 and trm7 increases codon-specific dwell times which indicates their potential tRNA targets, and highlights effects of nucleotide modifications on ribosome decoding rate.

Strategy and Performance Evaluation of Low-Frequency Variant Calling for SARS-CoV-2 Using Targeted Deep Illumina Sequencing

The ongoing COVID-19 pandemic, caused by SARS-CoV-2, constitutes a tremendous global health issue. Continuous monitoring of the virus has become a cornerstone to make rational decisions on implementing societal and sanitary measures to curtail the virus spread. Additionally, emerging SARS-CoV-2 variants have increased the need for genomic surveillance to detect particular strains because of their potentially increased transmissibility, pathogenicity and immune escape. Targeted SARS-CoV-2 sequencing of diagnostic and wastewater samples has been explored as an epidemiological surveillance method for the competent authorities. Currently, only the consensus genome sequence of the most abundant strain is taken into consideration for analysis, but multiple variant strains are now circulating in the population. Consequently, in diagnostic samples, potential co-infection(s) by several different variants can occur or quasispecies can develop during an infection in an individual. In wastewater samples, multiple variant strains will often be simultaneously present. Currently, quality criteria are mainly available for constructing the consensus genome sequence, and some guidelines exist for the detection of co-infections and quasispecies in diagnostic samples. The performance of detection and quantification of low-frequency variants using whole genome sequencing (WGS) of SARS-CoV-2 remains largely unknown. Here, we evaluated the detection and quantification of mutations present at low abundances using the mutations defining the SARS-CoV-2 lineage B.1.1.7 (alpha variant) as a case study. Real sequencing data were in silico modified by introducing mutations of interest into raw wild-type sequencing data, or by mixing wild-type and mutant raw sequencing data, to construct mixed samples subjected to WGS using a tiling amplicon-based targeted metagenomics approach and Illumina sequencing. As anticipated, higher variation and lower sensitivity were observed at lower coverages and allelic frequencies. We found that detection of all low-frequency variants at an abundance of 10, 5, 3, and 1%, requires at least a sequencing coverage of 250, 500, 1500, and 10,000×, respectively. Although increasing variability of estimated allelic frequencies at decreasing coverages and lower allelic frequencies was observed, its impact on reliable quantification was limited. This study provides a highly sensitive low-frequency variant detection approach, which is publicly available at https://galaxy.sciensano.be, and specific recommendations for minimum sequencing coverages to detect clade-defining mutations at certain allelic frequencies. This approach will be useful to detect and quantify low-frequency variants in both diagnostic (e.g., co-infections and quasispecies) and wastewater [e.g., multiple variants of concern (VOCs)] samples.

George K, Lasek-Nesselquist E, Sikkema RS, Oude Munnink B, Koopmans M, Brytting M, Sudha rani V, Pavani S, Smura T, Heim A, Kurkela S, Umair M, Salman M, Bartolini B, Rueca M, Drosten C, Wolff T, Silander O, Eggink D, Reusken C, Vennema H, Park A, Carrington C, Sahadeo N, Carr M, Gonzalez G, de Oliveira T, Faria N, Rambaut A, Kraemer MUG. Tracking the international spread of SARS-CoV-2 lineages B.1.1.7 and B.1.351/501Y-V2 with grinch

Late in 2020, two genetically-distinct clusters of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) with mutations of biological concern were reported, one in the United Kingdom and one in South Africa. Using a combination of data from routine surveillance, genomic sequencing and international travel we track the international dispersal of lineages B.1.1.7 and B.1.351 (variant 501Y-V2). We account for potential biases in genomic surveillance efforts by including passenger volumes from location of where the lineage was first reported, London and South Africa respectively. Using the software tool grinch (global report investigating novel coronavirus haplotypes), we track the international spread of lineages of concern with automated daily reports, Further, we have built a custom tracking website (cov-lineages.org/global_report.html) which hosts this daily report and will continue to include novel SARS-CoV-2 lineages of concern as they are detected.

Cloud Computing Enabled Big Multi-Omics Data Analytics

High-throughput experiments enable researchers to explore complex multifactorial diseases through large-scale analysis of omics data. Challenges for such high-dimensional data sets include storage, analyses, and sharing. Recent innovations in computational technologies and approaches, especially in cloud computing, offer a promising, low-cost, and highly flexible solution in the bioinformatics domain. Cloud computing is rapidly proving increasingly useful in molecular modeling, omics data analytics (eg, RNA sequencing, metabolomics, or proteomics data sets), and for the integration, analysis, and interpretation of phenotypic data. We review the adoption of advanced cloud-based and big data technologies for processing and analyzing omics data and provide insights into state-of-the-art cloud bioinformatics applications.

LGAAP: Leishmaniinae Genome Assembly and Annotation Pipeline

We present the LGAAP computational pipeline, which was successfully used to assemble six genomes of the parasite subfamily Leishmaniinae to chromosome-scale completeness from a combination of long- and short-read sequencing data. LGAAP is open source, and we suggest that it may easily be ported for assembly of any genome of comparable size (∼35 Mb).

Chromosome-Scale Assembly of the Complete Genome Sequence of Leishmania (Mundinia) martiniquensis, Isolate LSCM1, Strain LV760

Leishmania (Mundinia) martiniquensis is a kinetoplastid parasite that was first isolated in 1995 on Martinique. We report the first complete genome for Leishmania martiniquensis from Asia, isolate LSCM1, strain LV760, which was sequenced using combined short-read and long-read technologies. This will facilitate greater understanding of the evolution of the geographically dispersed subgenus Mundinia.

Construction of Condition-Specific Gene Regulatory Network Using Kernel Canonical Correlation Analysis

Gene expression profile or transcriptome can represent cellular states, thus understanding gene regulation mechanisms can help understand how cells respond to external stress. Interaction between transcription factor (TF) and target gene (TG) is one of the representative regulatory mechanisms in cells. In this paper, we present a novel computational method to construct condition-specific transcriptional networks from transcriptome data. Regulatory interaction between TFs and TGs is very complex, specifically multiple-to-multiple relations. Experimental data from TF Chromatin Immunoprecipitation sequencing is useful but produces one-to-multiple relations between TF and TGs. On the other hand, co-expression networks of genes can be useful for constructing condition transcriptional networks, but there are many false positive relations in co-expression networks. In this paper, we propose a novel method to construct a condition-specific and combinatorial transcriptional network, applying kernel canonical correlation analysis (kernel CCA) to identify multiple-to-multiple TF-TG relations in certain biological condition. Kernel CCA is a well-established statistical method for computing the correlation of a group of features vs. another group of features. We, therefore, employed kernel CCA to embed TFs and TGs into a new space where the correlation of TFs and TGs are reflected. To demonstrate the usefulness of our network construction method, we used the blood transcriptome data for the investigation on the response to high fat diet in a human and an arabidopsis data set for the investigation on the response to cold/heat stress. Our method detected not only important regulatory interactions reported in previous studies but also novel TF-TG relations where a module of TF is regulating a module of TGs upon specific stress.

George K, Lasek-Nesselquist E, Sikkema RS, Oude Munnink B, Koopmans M, Brytting M, Sudha rani V, Pavani S, Smura T, Heim A, Kurkela S, Umair M, Salman M, Bartolini B, Rueca M, Drosten C, Wolff T, Silander O, Eggink D, Reusken C, Vennema H, Park A, Carrington C, Sahadeo N, Carr M, Gonzalez G, de Oliveira T, Faria N, Rambaut A, Kraemer MUG. Tracking the international spread of SARS-CoV-2 lineages B.1.1.7 and B.1.351/501Y-V2

Late in 2020, two genetically-distinct clusters of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) with mutations of biological concern were reported, one in the United Kingdom and one in South Africa. Using a combination of data from routine surveillance, genomic sequencing and international travel we track the international dispersal of lineages B.1.1.7 and B.1.351 (variant 501Y-V2). We account for potential biases in genomic surveillance efforts by including passenger volumes from location of where the lineage was first reported, London and South Africa respectively. Using the software tool grinch (global report investigating novel coronavirus haplotypes), we track the international spread of lineages of concern with automated daily reports, Further, we have built a custom tracking website (cov-lineages.org/global_report.html) which hosts this daily report and will continue to include novel SARS-CoV-2 lineages of concern as they are detected.

Scalable Declarative HEP Analysis Workflows for Containerised Compute Clouds

We describe a novel approach for experimental High-Energy Physics (HEP) data analyses that is centred around the declarative rather than imperative paradigm when describing analysis computational tasks. The analysis process can be structured in the form of a Directed Acyclic Graph (DAG), where each graph vertex represents a unit of computation with its inputs and outputs, and the graph edges describe the interconnection of various computational steps. We have developed REANA, a platform for reproducible data analyses, that supports several such DAG workflow specifications. The REANA platform parses the analysis workflow and dispatches its computational steps to various supported computing backends (Kubernetes, HTCondor, Slurm). The focus on declarative rather than imperative programming enables researchers to concentrate on the problem domain at hand without having to think about implementation details such as scalable job orchestration. The declarative programming approach is further exemplified by a multi-level job cascading paradigm that was implemented in the Yadage workflow specification language. We present two recent LHC particle physics analyses, ATLAS searches for dark matter and CMS jet energy correction pipelines, where the declarative approach was successfully applied. We argue that the declarative approach to data analyses, combined with recent advancements in container technology, facilitates the portability of computational data analyses to various compute backends, enhancing the reproducibility and the knowledge preservation behind particle physics data analyses.

Bioimage analysis workflows: community resources to navigate through a complex ecosystem

Workflows are the keystone of bioimage analysis, and the NEUBIAS (Network of European BioImage AnalystS) community is trying to gather the actors of this field and organize the information around them.  One of its most recent outputs is the opening of the F1000Research NEUBIAS gateway, whose main objective is to offer a channel of publication for bioimage analysis workflows and associated resources. In this paper we want to express some personal opinions and recommendations related to finding, handling and developing bioimage analysis workflows.  The emergence of "big data" in bioimaging and resource-intensive analysis algorithms make local data storage and computing solutions a limiting factor. At the same time, the need for data sharing with collaborators and a general shift towards remote work, have created new challenges and avenues for the execution and sharing of bioimage analysis workflows. These challenges are to reproducibly run workflows in remote environments, in particular when their components come from different software packages, but also to document them and link their parameters and results by following the FAIR principles (Findable, Accessible, Interoperable, Reusable) to foster open and reproducible science. In this opinion paper, we focus on giving some directions to the reader to tackle these challenges and navigate through this complex ecosystem, in order to find and use workflows, and to compare workflows addressing the same problem. We also discuss tools to run workflows in the cloud and on High Performance Computing resources, and suggest ways to make these workflows FAIR.

The National Ecological Observatory Network's soil metagenomes: assembly and basic analysis

The National Ecological Observatory Network (NEON) annually performs shotgun metagenomic sequencing to sample genes within soils at 47 sites across the United States. NEON serves as a valuable educational resource, thanks to its open data policies and programming tutorials, but there is currently no introductory tutorial for performing analyses with the soil shotgun metagenomic dataset. Here, we describe a workflow for processing raw soil metagenome sequencing reads using the Sunbeam bioinformatics pipeline. The workflow includes cleaning and processing raw reads, taxonomic classification, assembly into contigs, annotation of predicted genes using custom protein databases, and exporting assemblies to the KBase platform for downstream analysis. This workflow is designed to be robust to annual data releases from NEON, and the underlying Snakemake framework can manage complex software dependencies. The workflow presented here aims to increase the accessibility of NEON's shotgun metagenome data, which can provide important clues about soil microbial communities and their ecological roles.

RADSex: A computational workflow to study sex determination using restriction site-associated DNA sequencing data

The study of sex determination and sex chromosome organization in nonmodel species has long been technically challenging, but new sequencing methodologies now enable precise and high-throughput identification of sex-specific genomic sequences. In particular, restriction site-associated DNA sequencing (RAD-Seq) is being extensively applied to explore sex determination systems in many plant and animal species. However, software specifically designed to search for and visualize sex-biased markers using RAD-Seq data is lacking. Here, we present RADSex, a computational analysis workflow designed to study the genetic basis of sex determination using RAD-Seq data. RADSex is simple to use, requires few computational resources, makes no prior assumptions about the type of sex-determination system or structure of the sex locus, and offers convenient visualization through a dedicated R package. To demonstrate the functionality of RADSex, we re-analysed a published data set of Japanese medaka, Oryzias latipes, where we uncovered a previously unknown Y chromosome polymorphism. We then used RADSex to analyse new RAD-Seq data sets from 15 fish species spanning multiple taxonomic orders. We identified the sex determination system and sex-specific markers in six of these species, five of which had no known sex-markers prior to this study. We show that RADSex greatly facilitates the study of sex determination systems in nonmodel species thanks to its speed of analyses, low resource usage, ease of application and visualization options. Furthermore, our analysis of new data sets from 15 species provides new insights on sex determination in fish.