Microbial bioinformatics in one picture

Every analysis in this handbook turns DNA sequence into a biological statement. Knowing where each step sits in that chain makes it easier to choose methods, and to notice when a result is coming from the pipeline rather than from the biology.

Two questions, two kinds of data

Most microbial sequencing projects ask one of two questions. The first is who is there and what are they doing in a mixed community: a stool sample, a soil core, a fermenting food. This is the domain of microbiome analysis. It uses either marker-gene (amplicon) sequencing, usually of the 16S rRNA gene for bacteria and archaea or ITS and 18S for fungi and other eukaryotes, or shotgun metagenomics, where all DNA in the sample is sequenced. The second question is what is this organism: the complete genome of a single bacterial isolate, its taxonomy, its genes, its plasmids and how it differs from its relatives. This is isolate genomics.

The two fields share tools and ideas, such as quality control of reads, reference databases and the need for careful controls, but their statistics are very different. A microbiome study compares many samples and has to deal with sparse, compositional count tables. A genome project usually deals with a handful of genomes, and the hard part is getting each sequence right.

The chain of steps

For a microbiome study the chain runs: study design and sampling, DNA extraction, library preparation and sequencing, read quality control, feature construction (amplicon sequence variants or taxonomic profiles), taxonomy and function assignment, filtering and contaminant removal, normalisation and transformation, then the ecological and statistical questions: diversity, community differences, differentially abundant features, associations, networks and prediction. Chapters 2 to 13 follow this chain.

For an isolate genome it runs: culture and DNA extraction, sequencing with short reads, long reads or both, read quality control, assembly, polishing, quality assessment, taxonomic placement, annotation, then targeted questions such as typing, plasmids, resistance genes and comparison with reference strains. Chapters 14 to 21 follow that chain.

What "literature-based" means here

Each chapter explains the reasoning behind a step, the main methods and the trade-offs between them, with references to the papers that introduced or tested them. Where the field has not settled on one answer, for example on rarefaction or on which differential abundance test to trust, the chapter says so and describes what careful authors do instead of pretending there is a single correct choice. Boxes marked In practice give concrete recommendations; boxes marked Pitfall describe mistakes that are common in published work.

How the chapters map onto the tools

MicroRes follows the microbiome chain step by step. Its numbered steps (01 Import to 13 Biomarkers and prediction, then F1 to F6 for functional data) correspond to chapters 4 to 13. BactFlow follows the genome chain: its QC, Assembly and Post-assembly modules correspond to chapters 14 to 21. You can read the handbook on its own, but every chapter ends with a short note on where the step lives in the tools.

How to use these docs to learn

Each chapter starts with what you will learn, explains the idea before the tool, works through at least one small example with real numbers, lists the common mistakes and ends with questions to test yourself. Book references in the text and in the reference list name the chapter to read next: the docs are a map, the books are the territory. A sensible order for a newcomer is Part I (foundations), then the part for your data type, then Part IV on practice. The statistics are best learned from Xia, Sun & Chen, Statistical Analysis of Microbiome Data with R and Xia & Sun, Bioinformatic and Statistical Analysis of Microbiome Data for microbiomes, Efron & Hastie, Computer Age Statistical Inference for large-scale testing, and James et al., An Introduction to Statistical Learning for machine learning; the sequence side from Ismail, Bioinformatics: A Practical Guide to NGS Data Analysis, Low & Tammi, Practical Bioinformatics for Beginners and Christensen (ed.), Introduction to Bioinformatics in Microbiology.

Study design, controls and metadata

The most important decisions in a microbiome study are made before any DNA is extracted. Design determines which questions can be answered at all; controls determine whether the answer can be trusted.

Replication and the unit of analysis

The unit of replication is the thing that was independently assigned to a treatment or sampled from a population: a person, an animal, a pen, a field plot. Several samples from the same unit (repeated time points, technical replicates, several gut sections) are not independent and must be modelled as such, usually with a random effect or by averaging. Treating them as independent inflates the apparent sample size and produces false positives, a problem known as pseudoreplication.

Microbial communities vary a lot between individuals, so effects that are clear within an individual can be small compared with the differences between individuals. Longitudinal designs, where each subject is its own control, are therefore much more powerful than cross-sectional ones for many questions. Animal studies add a well-known complication: animals housed together share microbes through coprophagy and the environment, so the cage, not the animal, is often the true unit.

Confounders that matter in practice

Several factors are known to shift community composition strongly enough to create spurious associations: antibiotics and other medication (proton-pump inhibitors, metformin, laxatives), diet in the days before sampling, stool consistency and transit time, age, and host genotype in animal models. Technical factors matter just as much: the extraction kit, the sequencing run, the person who processed the samples, and the time a sample spent at room temperature. The safest approach is to randomise samples across extraction batches and sequencing runs, so that technical batches cannot line up with biological groups, and to record every one of these variables.

Negative and positive controls

Reagents and laboratory surfaces carry their own low-level microbial DNA. In samples with high biomass, such as stool, this "kitome" is swamped by real signal. In low-biomass samples, such as skin swabs, tissue biopsies, blood, milk or air, contaminants can make up a large share of the reads and have produced headline findings that later proved to be artefacts (Salter et al., 2014; Eisenhofer et al., 2019). Every extraction batch should include an extraction blank (no sample) and every PCR plate a no-template control. These controls are sequenced like real samples and used later to identify contaminants (chapter 6).

Positive controls are mock communities: defined mixtures of known strains or of their DNA, such as the ZymoBIOMICS standards. Sequencing a mock community shows how far your workflow distorts the true composition. Extraction bias against Gram-positive cells, primer mismatches and GC bias all show up as a gap between expected and observed abundances. A mock community is also the best way to compare two pipelines on your own data.

Relative and absolute abundance

Sequencing returns a sample of reads whose total depends on the library, not on how many microbes were in the sample. Counts are therefore relative: if one taxon doubles, every other taxon appears to shrink. Several methods recover absolute numbers. Flow cytometry measures total cell load per gram (Vandeputte et al., 2017), qPCR of the 16S gene measures total gene copies, and spike-ins add a known number of foreign cells or synthetic DNA before extraction. If the scientific question is about absolute change (did the antibiotic reduce the total bacterial load?), plan for one of these from the start.

Metadata

A sample without metadata is almost worthless. Record subject, time point, treatment, the confounders above, extraction batch, sequencing run and any deviation from protocol, in a single tidy table with one row per sample. Use controlled vocabularies, units in the column names and no merged cells. Standards such as MIxS from the Genomic Standards Consortium describe the minimum information expected when data are deposited in public archives, and following them from the start saves work at publication time (Wilkinson et al., 2016).

Sequencing technologies and read quality

Short reads are cheap and accurate; long reads span repeats and reveal structure. Most choices downstream follow from which kind of read you have and from how good it is.

Short reads

Illumina instruments produce reads of 100 to 300 bases, often as pairs from the two ends of a fragment. Per-base accuracy is high (typically better than 99.9 percent) and errors are mostly substitutions, which rise towards the end of each read. Short reads are the standard for amplicon sequencing, for shotgun metagenomic profiling and for polishing bacterial genomes. Their weakness is length: any repeat longer than a read or a fragment, such as rRNA operons, insertion sequences and duplicated genes, cannot be resolved, so short-read genome assemblies usually break into tens to hundreds of contigs.

Long reads

Oxford Nanopore (ONT) reads DNA as it passes through a protein pore and can produce reads from a few kilobases to more than a hundred kilobases. Older chemistries had error rates of several percent, dominated by insertions and deletions in homopolymers. The R10.4.1 pores with recent basecallers reach modal accuracies above 99 percent, enough to produce near-finished bacterial genomes without short-read polishing (Sereika et al., 2022). PacBio HiFi reads are 10 to 25 kb long with accuracy above 99.9 percent, from circular consensus sequencing of the same molecule. Older PacBio CLR reads are longer but much noisier. A long read that spans every repeat in a genome allows assembly into one circular chromosome plus plasmids.

Quality scores

Each base comes with a Phred quality score Q, where the probability that the base is wrong is

P(error) = 10^−Q/10

so Q20 means a 1 in 100 chance of error and Q30 a 1 in 1000 chance. Quality plots per position show where reads degrade; read-length histograms show whether long-read libraries contain the long molecules you paid for. Simple summary tools such as SeqKit report read counts, total bases, N50 and mean quality per file (Shen et al., 2016); NanoPack does the same with plots tailored to long reads (De Coster et al., 2018).

Ismail works through FASTQ encodings, quality plots and trimming in detail, and Low & Tammi give a hands-on version of the same checks (Ismail, Bioinformatics: A Practical Guide to NGS Data Analysis, ch. 1; Low & Tammi, Practical Bioinformatics for Beginners, ch. 3).

What read QC should do

• Remove adapters and primers. Adapter read-through happens whenever a fragment is shorter than the read. In amplicon data the primer sequences must also be removed because they do not represent the template.
• Trim or filter by quality. For short reads, trim low-quality tails (fastp, Trimmomatic). For long reads, filter whole reads below a mean quality or length threshold rather than trimming.
• Remove host reads in shotgun samples from host-associated sites by mapping to the host genome, both to save compute and for privacy in human studies.
• Check depth. Samples far below the run median are often failed libraries.

How much sequencing?

For 16S amplicons, 10,000 to 50,000 reads per sample usually captures the dominant community; rare taxa need more. For shotgun profiling of gut samples, a few million read pairs per sample are enough for species-level taxonomic profiles, while functional profiling and genome reconstruction from metagenomes benefit from 10 to 20 million or more. For isolate genomes, 50 to 100-fold short-read coverage or 30 to 50-fold long-read coverage is typical; very high long-read coverage can be subsampled to the longest, best reads.

From reads to features: amplicon sequencing

An amplicon study turns millions of reads into a table of counts: rows are sequence features, columns are samples. How those features are defined decides the resolution of everything that follows.

Choosing the marker and region

The 16S rRNA gene has nine hypervariable regions (V1 to V9) separated by conserved stretches where universal primers bind. Short-read studies sequence one or two regions, commonly V4 or V3–V4. No region separates all taxa equally well: V4 is robust and widely used but resolves few species, and primers differ in which groups they miss. Full-length 16S from long reads resolves many more taxa to species and sometimes strain level (Johnson et al., 2019). For fungi the ITS1 or ITS2 region is standard; its length varies a lot between taxa, which affects both PCR and read merging.

OTUs and ASVs

Early pipelines clustered reads into operational taxonomic units (OTUs) at 97 percent identity to absorb sequencing errors. Clustering merges real biological variants and gives results that depend on the dataset. Denoising methods instead learn the error profile of a run and infer the exact biological sequences present, called amplicon sequence variants (ASVs). DADA2 fits a parametric error model to each run and tests whether each unique sequence can be explained as an error from a more abundant one (Callahan et al., 2016). ASVs are reproducible across studies because they are sequences, not cluster labels, and are now the standard unit (Callahan, McMurdie & Holmes, 2017).

A typical DADA2 workflow

1. Remove primers; inspect quality profiles; trim reads where quality collapses, leaving enough overlap to merge pairs.
2. Filter reads with too many expected errors.
3. Learn error rates per run, then denoise forward and reverse reads separately.
4. Merge pairs, build the sequence table and remove chimeras, which are PCR artefacts joined from two parent sequences.
5. Track how many reads survive each step. Large losses at merging usually mean the reads were trimmed too short to overlap.

Xia & Sun devote three chapters to the shift from clustering to denoising: how OTU clustering grew out of numerical taxonomy, why a fixed 97 % threshold merges real variants and splits error clouds unevenly, and how DADA2 and Deblur instead model the error process so that every exact sequence variant is either real or an error (Xia & Sun, Bioinformatic and Statistical Analysis of Microbiome Data, ch. 6–8). Mitra's volume compares the main pipelines side by side (Mitra (ed.), Metagenomic Data Analysis, ch. 4).

Taxonomic assignment

Each ASV is assigned a lineage by comparison with a reference database, typically SILVA (Quast et al., 2013) or GTDB-derived 16S sets for bacteria and archaea, and UNITE for fungi. The naïve Bayesian classifier counts shared k-mers between the query and reference sequences and reports a bootstrap confidence at each rank (Wang et al., 2007). Classifiers trained on the exact amplified region perform better than those trained on full-length genes (Bokulich et al., 2018). Species-level calls from short amplicons should be treated with caution; genus is usually the finest reliable rank.

Phylogeny

A phylogenetic tree of the ASVs is needed for phylogeny-aware measures such as Faith's phylogenetic diversity and UniFrac. Trees are usually built by aligning ASVs with MAFFT and inferring a tree with FastTree, or by inserting ASVs into a reference tree.

Taxonomic profiling of shotgun metagenomes

Shotgun sequencing reads every genome in the sample, not one gene. Profilers turn those reads into species abundances in two different ways, and the difference matters when you interpret the numbers.

Read classifiers

Kraken 2 assigns each read to the lowest common ancestor of all reference genomes that share its k-mers (Wood, Lu & Langmead, 2019). It is fast and classifies reads from any domain in the database, including host, fungi and viruses. Because many reads stop at genus or family when species share sequence, Bracken redistributes them to species using the expected read distribution of each genome (Lu et al., 2017). The output is a read share: the fraction of sequenced reads that came from each taxon. Larger genomes contribute more reads per cell, so read share is not the same as cell share.

Zhong Wang explains the algorithms behind both families, from k-mer matching and lowest common ancestor assignment to marker-gene and binning approaches (Zhong Wang, Introduction to Computational Metagenomics, ch. 6).

Marker-gene profilers

MetaPhlAn maps reads only to a curated set of clade-specific marker genes, about a million of them in MetaPhlAn 4, and estimates abundances from marker coverage (Blanco-Míguez et al., 2023). Because coverage is normalised by marker length, the result approximates the relative abundance of genomes, closer to a cell share. MetaPhlAn 4 defines species-level genome bins (SGBs) from large collections of metagenome-assembled genomes, so it can profile species that have never been cultured (Pasolli et al., 2019). Reads from organisms without markers, and from host or unknown sources, are simply not counted, so profiles sum to 100 percent of the known fraction.

Database choice is a modelling choice

Both approaches see only what is in their database. A taxon missing from the reference set is either invisible or attributed to its nearest relative. When comparing studies, use the same tool and database version; when reporting, give both. For environments far from the human gut, such as soils or fermented foods, database completeness is often the main limit, and assembly-based approaches (binning metagenome-assembled genomes) recover organisms that profilers cannot see (Quince et al., 2017).

Validation with a mock community

Running a mock community through the same pipeline shows whether expected species are recovered, whether unexpected ones appear (false positives) and how far abundances drift from the truth. Recovering all expected species with few false positives is a better test than matching the exact proportions, which extraction bias distorts anyway.

Filtering, quality control and contaminants

A raw feature table contains real taxa, off-target sequences, contaminants, artefacts and samples that failed. Filtering removes what cannot be biology, but every filter is a decision that must be reported.

Off-target and unassigned features

16S primers also amplify chloroplast and mitochondrial rRNA, which are bacterial in origin. In plant, food or host-tissue samples these can dominate. Features assigned to Chloroplast (order) or Mitochondria (family), and features unassigned at kingdom or phylum level, are normally removed before analysis.

Prevalence and abundance filters

Most features in a microbiome table are rare: seen in one or two samples at a few reads. They add noise, inflate the number of tests in differential abundance and destabilise network inference. Common filters keep features present in at least a minimum number or fraction of samples (for example 5 to 10 percent) and above a minimum mean relative abundance. These filters should be applied without looking at the group labels, otherwise they bias downstream tests. Alpha diversity is the exception: richness estimators need the rare features, so compute diversity on the unfiltered table.

Filtering is a modelling decision, not housekeeping. Xia, Sun & Chen show how prevalence filters reduce the number of tests and the share of zeros, which improves power, but also how a filter applied separately to each group creates differences that are not there (Xia, Sun & Chen, Statistical Analysis of Microbiome Data with R, ch. 7). Apply one filter to all samples, decide it before testing, and report it.

Sample-level checks

Samples with very low read counts after processing are usually failed libraries. Plot depth per sample and per group: if one group systematically has lower depth, depth becomes a confounder. Ordinations of all samples, including controls, also reveal outliers and mislabelled samples. Controls should cluster away from real samples.

Contaminant identification

Contaminant DNA comes from reagents, kits and the lab, and it behaves in two recognisable ways (Davis et al., 2018):

• Frequency. A contaminant is added in roughly constant amounts, so its relative frequency is inversely related to the total DNA concentration of the sample. Real taxa do not show this pattern.
• Prevalence. A contaminant is more prevalent in negative controls than in real samples.

The decontam method tests both patterns per feature and flags features with a score below a threshold. Flagged features should be inspected, not deleted blindly: common gut taxa can appear in blanks through well-to-well cross-contamination, and removing them from real samples would be wrong. Known reagent contaminant genera, such as Ralstonia, Bradyrhizobium and Sphingomonas, deserve suspicion in low-biomass samples (Salter et al., 2014).

Compositionality, normalisation and transformation

Sequencing counts carry information only about ratios between taxa. Normalisation corrects for unequal sequencing depth; transformation puts the data on a scale where standard statistics behave. They solve different problems and are best chosen separately.

Why the data are compositional

The total number of reads per sample is set by the sequencing run, not by the biology. Counts are therefore a random sample of a composition, and only ratios between parts are meaningful (Gloor et al., 2017). Two consequences follow. First, correlations between relative abundances are biased: when one taxon increases, the others must decrease, creating negative correlations that do not exist in the community. Second, a taxon can appear to increase only because others decreased. Aitchison's log-ratio framework (Aitchison, 1986) is the principled way around this.

The full theory, from Aitchison's simplex to CLR, ILR and ALR transformations and their use in tests, is in Xia, Sun & Chen, Statistical Analysis of Microbiome Data with R, ch. 10, and Xia & Sun, Bioinformatic and Statistical Analysis of Microbiome Data, ch. 14.

Normalisation: dealing with depth

• Total-sum scaling (TSS) divides by library size to give proportions. Simple, but sensitive to a few dominant taxa.
• Rarefying subsamples every sample to the same depth. It discards data and adds randomness, and is inadmissible for differential abundance testing (McMurdie & Holmes, 2014). It remains defensible for diversity estimates, where unequal effort biases richness directly (Weiss et al., 2017; Schloss, 2024).
• Size-factor methods estimate a scaling factor from features assumed not to change: the median-of-ratios (RLE) of DESeq2 (Love, Huber & Anders, 2014), TMM (Robinson & Oshlack, 2010), upper quartile, and GMPR, designed for zero-inflated microbiome data (Chen et al., 2018). CSS scales by a data-driven quantile (Paulson et al., 2013).
• Absolute scaling multiplies proportions by measured cell load or uses spike-ins (chapter 2). It is the only approach that recovers real changes in total abundance.

Transformation: getting a usable scale

Counts are skewed and their variance grows with the mean, so most methods based on Euclidean distance or linear models need a transformation first.

• Centred log-ratio (CLR): the log of each count divided by the geometric mean of the sample. Euclidean distances on CLR values give the Aitchison distance. Zeros must first be replaced, by a pseudocount or better by a Bayesian-multiplicative method.
• Robust CLR computes the geometric mean over non-zero values only and handles sparse data better (Martino et al., 2019).
• Hellinger (square root of proportions) makes Euclidean methods such as PCA and RDA suitable for abundance data (Legendre & Gallagher, 2001).
• Log counts per million and variance-stabilising transformation come from RNA-seq and work well when depth is moderately uniform.
• Presence/absence keeps only occurrence and is useful for rare biosphere questions.

Alpha diversity

Alpha diversity summarises a single community in one number: how many kinds of organisms it holds and how evenly they are spread. It is the most reported microbiome statistic and one of the easiest to get wrong.

Richness, evenness and Hill numbers

Richness counts the taxa present. Shannon entropy and the Simpson index also weight taxa by their abundance, so they reflect evenness as well. These indices sit on different scales, which makes them hard to compare. Hill numbers put them on one scale, the effective number of taxa (Hill, 1973):

^qD = ( Σ p_i^q )^1/(1−q)

With q = 0 this is richness, q = 1 (as a limit) is the exponential of Shannon entropy, and q = 2 is the inverse Simpson index. Plotting D across q gives a diversity profile: a community with many rare taxa has high D at q = 0 that falls steeply as q rises. Faith's phylogenetic diversity sums the branch lengths of the tree spanned by the taxa present and rewards communities that cover distant lineages (Faith, 1992).

Sampling effort

Observed richness grows with sequencing depth, so samples sequenced more deeply look more diverse. There are three common ways to compare at equal effort. Rarefying once to the smallest depth is simple but adds noise and discards data. Repeated rarefaction averages over many subsamples and reduces that noise (Schloss, 2024). Coverage-based standardisation compares samples at equal sample completeness rather than equal read count, which is fairer when communities differ in how many rare taxa they hold (Chao & Jost, 2012). Richness is the most sensitive to effort; Shannon and Simpson much less (Willis, 2019).

The missing-species problem has a famous history. Fisher was asked how many new butterfly species a collector would find in a second year of trapping in Malaya; his answer, refined by Good and Toulmin, is the alternating series above. Efron & Hastie use it to open their discussion of empirical Bayes (Efron & Hastie, Computer Age Statistical Inference, section 6.2). MicroRes reports these quantities for every sample in step 06 ("How many taxa are still unseen?"), together with the sample size a follow-up study needs (chapter 2).

Testing

Diversity values are one number per sample and can be compared with ordinary tests: Wilcoxon or Kruskal–Wallis for groups, linear models when there are covariates, and mixed models when subjects are measured repeatedly. Report effect sizes, not only p-values: a statistically significant difference of two taxa in a richness of 300 is rarely meaningful.

Beta diversity and ordination

Beta diversity asks how different two communities are. A distance between every pair of samples, an ordination to look at it, and a permutation test to judge it form the backbone of most microbiome papers.

Choosing a distance

The distance encodes what you consider "different":

• Bray–Curtis weights abundant taxa and ignores shared absences. It is the ecological default for abundance data but is affected by depth unless the data are normalised or rarefied.
• Jaccard uses presence/absence only and emphasises rare membership.
• Aitchison is the Euclidean distance between CLR-transformed samples. It respects compositionality and is scale-invariant. Robust Aitchison handles sparsity better (Martino et al., 2019).
• UniFrac measures the fraction of tree branch length unique to one sample. Unweighted UniFrac uses presence/absence; weighted UniFrac uses abundance (Lozupone & Knight, 2005).

Different distances can lead to different conclusions, and that is information. A treatment that shifts unweighted UniFrac but not Bray–Curtis changes which rare lineages are present without changing the dominant taxa.

Ordination is unsupervised learning: PCA, PCoA and NMDS all try to show high-dimensional distances in two dimensions, and the share of variance explained by each axis tells you how much they hide (James et al., An Introduction to Statistical Learning, ch. 12; Xia & Sun, Bioinformatic and Statistical Analysis of Microbiome Data, ch. 10).

Ordination

Principal coordinates analysis (PCoA) places samples in a few dimensions that preserve the chosen distance as well as possible. The share of variation on each axis tells you how much of the structure the plot shows. PCA is PCoA with Euclidean distance and works directly on CLR or Hellinger data. Non-metric multidimensional scaling (NMDS) preserves only the rank order of distances and reports a stress value; stress above about 0.2 means the picture is unreliable. Constrained ordinations (RDA, CCA, distance-based RDA) show only the variation explained by chosen metadata and are useful for asking how much a set of variables explains.

Testing: PERMANOVA and dispersion

PERMANOVA partitions the distance matrix among groups or covariates and assesses significance by permuting labels (Anderson, 2001). Its R² is the fraction of variation explained. In microbiome studies, single factors often explain 1 to 10 percent, which is still meaningful given how variable communities are. Two cautions apply. First, PERMANOVA is sensitive to differences in spread (dispersion) as well as location, so test dispersion with PERMDISP and report both. Second, repeated measures must be handled by restricting permutations within subject (strata); otherwise the test treats correlated samples as independent. Covariates are added sequentially, so their order matters unless marginal tests are requested.

Differential abundance

Which taxa differ between groups is the question most microbiome studies end on, and the one where methods disagree most. Using more than one method and reading the agreement is the honest answer.

Why it is hard

Microbiome counts are sparse (most entries are zero), over-dispersed, compositional and measured at unequal depth. A method must decide what "no change" means. Because only ratios are observed, a test on relative data can only say that a taxon changed relative to something: the rest of the community, a reference taxon or an assumed unchanged majority. When many taxa change in the same direction, methods that assume most taxa are unchanged break down.

Microbiome counts are also over-dispersed (the variance far exceeds the mean, unlike the Poisson model) and full of zeros. Faraway's chapter on count regression is the clearest introduction to Poisson, quasi-Poisson and negative binomial models (Faraway, Extending the Linear Model with R, ch. 5); Xia, Sun & Chen, Statistical Analysis of Microbiome Data with R, ch. 11–12, apply them to microbiome data, including zero-inflated and hurdle models. Benchmarks agree that no single method wins everywhere, which is why MicroRes runs several and reports the taxa they agree on.

The main families of methods

• Count models from RNA-seq. DESeq2 fits a negative binomial model with size factors and shrinks fold changes (Love, Huber & Anders, 2014). It is powerful but can be anti-conservative on very sparse data; size factors such as poscounts or GMPR help.
• Bias-corrected log-linear models. ANCOM-BC estimates and removes the sample-specific sampling fraction and reports fold changes with confidence intervals (Lin & Peddada, 2020). LinDA fits linear models to CLR data and corrects the compositional bias afterwards; it is fast and handles mixed models (Zhou et al., 2022).
• Log-ratio methods. ALDEx2 draws Monte Carlo instances from a Dirichlet distribution, applies CLR and tests each instance, so the uncertainty of low counts is propagated (Fernandes et al., 2014). It tends to be conservative.
• General linear models. MaAsLin2 fits per-feature models on transformed data with fixed and random effects (Mallick et al., 2021).

A comparison of 14 methods on 38 datasets found that they identified very different sets of significant taxa, and that ALDEx2 and ANCOM-II gave the most consistent results across datasets (Nearing et al., 2022). The practical lesson is to run several methods with different assumptions and give most weight to taxa they agree on.

Multiple testing and effect sizes

With hundreds of taxa tested, control the false discovery rate with the Benjamini–Hochberg procedure (Benjamini & Hochberg, 1995). Report the fold change and its uncertainty, the prevalence of the taxon in each group, and the method. A taxon present in 3 of 40 samples can be "significant" in some methods and is rarely worth a biological claim.

Associations and co-occurrence networks

Correlating taxa with each other or with host variables is easy to do and easy to over-interpret. Compositional methods remove the worst artefacts; careful reading removes most of the rest.

Taxa and host variables

Associations between taxa and continuous metadata, such as a metabolite, a clinical score or age, are usually tested with Spearman correlation on a transformed table, or with per-feature linear models that adjust for covariates. Partial correlations remove the shared influence of a third variable. Because every taxon is tested against every variable, multiple-testing correction is essential.

Taxon–taxon networks

Plain correlations between relative abundances are biased by compositionality (chapter 7). Two methods address this. SparCC estimates correlations between the underlying absolute abundances from log-ratio variances, assuming the network is sparse (Friedman & Alm, 2012). SPIEC-EASI goes further and estimates conditional dependence, through neighbourhood selection or the graphical lasso on CLR data, so that two taxa linked only through a third are not connected directly (Kurtz et al., 2015). SPIEC-EASI networks are sparser and more interpretable, at the cost of more computation and a need for many samples. Several hundred samples are a reasonable target, and networks built from 20 samples are mostly noise.

The problem is old: Karl Pearson warned in 1897 of "spurious correlation" between ratios that share a denominator. Proportions all share the sample total, so two unrelated taxa become negatively correlated when a third one blooms. SparCC and SPIEC-EASI estimate associations between log-ratios and add sparsity, assuming most taxa do not interact directly (Xia, Sun & Chen, Statistical Analysis of Microbiome Data with R, ch. 10). Graph concepts such as degree, hubs, modules and betweenness are explained for biological networks in Baxevanis, Bader & Wishart, Bioinformatics, ch. 13, and the ecology of microbial interactions in Zhong Wang, Introduction to Computational Metagenomics, ch. 10.

Reading a network

Common summaries are node degree, betweenness centrality, modules (clusters of densely connected taxa, for example from the Louvain algorithm (Blondel et al., 2008)) and the balance of positive and negative edges. Highly connected or central taxa are sometimes called keystone taxa, but a statistical hub is not the same as an ecological keystone (Berry & Widder, 2014). An edge means co-variation across samples: it can come from a shared niche, a shared response to diet or a direct interaction, and the network cannot tell these apart (Röttjers & Faust, 2018).

Comparing networks between groups and samples

Comparing summary numbers of two networks, such as density or modularity, says little about where they differ. Differential network analysis asks, for every pair of features, whether their association changes between the groups. With correlations this is a test of the difference between two Fisher-transformed coefficients (Fukushima, 2013; McKenzie et al., 2016). The pairs whose association changes form a differential network. Clustering that network gives modules of features that rewire together, for example a group of taxa that co-vary tightly in controls but fall apart after antibiotics. For each module, the change in connectivity can be tested by permuting group labels. That turns a picture into a statistic and highlights the modules that carry the difference. NetCoMi implements related comparisons for microbiome data (Peschel et al., 2021).

Networks can also be built between samples: each sample is linked to its most similar samples, and community detection splits the graph into clusters of samples with similar whole profiles. Cross-tabulating these clusters against a study variable shows whether groups form their own clusters or mix. This is a graph-based complement to ordination and PERMANOVA (chapter 9).

Differential associations need more samples than a single network: each group must support its own correlation estimates. With fewer than about 20 to 30 samples per group, few changes pass a false discovery threshold. Results are then best read as exploratory.

Machine learning and biomarkers

Predicting a phenotype from the microbiome tests whether the community carries information about it. Done carelessly, the same exercise produces impressive accuracies that vanish on new data.

Models that suit microbiome data

Microbiome datasets have many more features than samples. Regularised linear models handle this well: LASSO selects a small set of features, ridge shrinks all of them, and the elastic net mixes the two and copes better with correlated taxa (Zou & Hastie, 2005). Random forests capture non-linear effects and interactions and need little tuning (Breiman, 2001). In a systematic comparison, L2-regularised logistic regression performed about as well as more complex models on microbiome classification while being far easier to interpret (Topçuoğlu et al., 2020).

Validation without leakage

The model must be evaluated on samples it never saw during training, including during feature selection and tuning. Three rules prevent most leakage:

1. Do every data-dependent step inside the cross-validation loop. Filtering by prevalence is safe if it ignores labels; selecting features by their association with the outcome must happen within each training fold.
2. Keep related samples together. Repeated samples from one subject, or animals from one cage, must fall in the same fold (grouped cross-validation). Otherwise the model learns individuals, not the phenotype.
3. Tune in an inner loop. Hyperparameters chosen on the test folds give optimistic estimates; nested cross-validation avoids this.

Report the area under the ROC curve (or balanced accuracy) with its spread across folds, and compare with a model trained on shuffled labels. An external validation cohort is the strongest evidence.

James et al., An Introduction to Statistical Learning explains cross-validation and the bootstrap (ch. 5), regularised regression such as the lasso (ch. 6) and random forests (ch. 8) with the clearest diagrams in the literature; it is freely available from the authors.

Interpreting features

Model coefficients and importances show what the model used, not what causes the phenotype. Features that are selected in most folds are more trustworthy than those selected once. Linear discriminant analysis effect size (LEfSe) is a popular biomarker ranking (Segata et al., 2011), but it is a univariate screen, not a predictive model, and should not be reported as one.

Functional metagenomics

Taxonomy tells you who is there; function tells you what the community can do. Functional profiles add a layer of interpretation, and a layer of assumptions that should be stated.

Where functional profiles come from

• HUMAnN maps shotgun reads to the pangenomes of species detected by MetaPhlAn, then translates the remaining reads against UniRef proteins. It reports gene families and MetaCyc pathways, each split by the species that contribute them (stratified output) (Beghini et al., 2021).
• PICRUSt2 predicts gene content from 16S ASVs by placing them in a reference tree and inferring the genomes of their relatives (Douglas et al., 2020). It is a prediction, not a measurement, and is least reliable for environments poorly represented by sequenced genomes.
• eggNOG-mapper annotates predicted genes, usually from assemblies or metagenome-assembled genomes, with orthologous groups, KEGG orthologs, enzymes, CAZymes and GO terms (Cantalapiedra et al., 2021).

Cleaning and normalising

HUMAnN tables contain UNMAPPED, UNINTEGRATED and UNGROUPED rows that are bookkeeping, not functions; they are usually removed or reported separately. Gene families can be regrouped to KEGG orthologs or enzyme commission numbers for interpretation. Functional tables are compositional just like taxonomic ones, so the same normalisation logic applies (chapter 7). Copies per million (CPM) or relative abundance followed by a log transformation are common choices.

Testing functions

Per-function models with covariates and random effects, in the style of MaAsLin (Mallick et al., 2021), are standard. MaAsLin 3 separates abundance (how much, when present) from prevalence (whether present) and adds a correction for compositionality (Nickols et al., 2024). With thousands of functions, false discovery rate control is essential, and pathway-level results are easier to interpret than gene-family results.

Who carries a function

Stratified output links each function to its contributing taxa, so a change in a pathway can be traced to the species behind it. A function carried by many taxa is functionally redundant: losing one carrier changes little. Tian and colleagues define functional redundancy from the gap between taxonomic and functional diversity (Tian et al., 2020). Contributional diversity, the diversity of carriers for one function, distinguishes functions that depend on a single keystone carrier from those spread across the community. Bipartite species–function networks make these patterns visible; their nestedness (NODF) shows whether specialist carriers hold subsets of the functions of generalists (Almeida-Neto et al., 2008).

A functional profile from DNA describes potential; the community may not express it. The step from genes to activity is covered by metatranscriptomics and by function-based screens, in which DNA from the community is cloned and tested for an activity directly, as described by Charles, Liles & Sessitsch in Functional Metagenomics: Tools and Applications. For the computational side, Zhong Wang, Introduction to Computational Metagenomics, ch. 7, covers gene prediction and annotation in metagenomes and ch. 8 the recovery of genomes (binning), which is what makes "who does what" possible beyond reference databases.

Read quality control for isolate genomes

An isolate genome is only as good as the reads behind it. A few minutes spent on read statistics saves hours of puzzling over a broken assembly.

What to look at first

For every sample, check the number of reads, total bases, mean and N50 read length, and mean quality. Divide total bases by the expected genome size to get coverage. Short-read assemblers need roughly 50 to 100-fold coverage; long-read assemblers do well at 30 to 50-fold of good reads, and more is not always better because very deep data can be subsampled to the longest, most accurate reads.

• Illumina: quality should stay high across most of the read. Adapter content shows short inserts. Duplication can indicate low input DNA.
• Nanopore: the read-length distribution tells you whether the extraction preserved long DNA. A peak at a few hundred bases usually means sheared DNA or a failed barcode. Mean quality depends on chemistry and basecaller model.
• PacBio: HiFi reads cluster tightly in length and are very accurate; CLR reads are longer and noisier and need different assembler settings.

The same checks as in chapter 3 apply, with one addition: an isolate should contain one genome. A k-mer profile with one coverage peak and a read classifier that assigns nearly all reads to one species confirm it. Ismail, Bioinformatics: A Practical Guide to NGS Data Analysis, ch. 1, walks through FastQC and trimming, and Buffalo, Bioinformatics Data Skills, ch. 10, explains the FASTQ format and quality encodings at the level needed to write your own checks.

Pooling, filtering and trimming

Nanopore runs often write many small FASTQ files per barcode, which are concatenated into one file per sample. Unclassified reads and barcodes with very few reads should be inspected: a near-empty barcode is a failed sample, not a small genome. Long reads are filtered by minimum length (for example 1 kb) and mean quality, keeping the longest reads up to a target coverage. Short reads are adapter- and quality-trimmed with fastp or Trimmomatic (Chen et al., 2018; Bolger et al., 2014).

Contamination at the read level

A quick taxonomic screen of a subset of reads (for example with Kraken 2) reveals mixed cultures or contamination before assembly. Two species in one culture produce an assembly roughly twice the expected size and high contamination scores later (chapter 17).

Genome assembly

Assembly reconstructs a genome from overlapping reads. Repeats are what make it difficult, and read length is what solves them.

Graphs

Short-read assemblers build a de Bruijn graph: reads are broken into k-mers, and k-mers that overlap by k−1 bases are linked. A unique region becomes a simple path; a repeat becomes a node where paths cross. SPAdes uses several k-mer sizes and paired-end information to resolve as much as possible (Bankevich et al., 2012). Long-read assemblers mostly use overlap or repeat graphs. Flye builds a repeat graph from approximate overlaps and resolves repeats using reads that span them (Kolmogorov et al., 2019). When a read is longer than every repeat, the graph collapses to one circle per replicon.

A worked de Bruijn graph, with the repeat that breaks it, is in chapter 31 (fungal assembly); the logic is the same for bacteria, whose main repeats are rRNA operons (about 5 kb, often 5–10 copies) and insertion sequences (1–2.5 kb). Reads longer than the longest repeat resolve the whole chromosome, which is why a single Nanopore or HiFi run now routinely gives complete bacterial genomes. Ismail, Bioinformatics: A Practical Guide to NGS Data Analysis, ch. 3, compares greedy, overlap–layout–consensus and de Bruijn assemblers; Low & Tammi, Practical Bioinformatics for Beginners, ch. 6, compares short-read, long-read and hybrid assemblies of one genome.

Choosing an assembler

Unicycler builds a short-read assembly graph and uses long reads to bridge its repeats, which gives complete genomes with high base accuracy when both data types are good (Wick et al., 2017). Benchmarks of long-read assemblers show that each has characteristic failure modes, such as missing small plasmids or duplicating the start and end of a circular contig, and that combining several assemblies into a consensus (Trycycler) gives the most reliable finished genomes (Wick & Holt, 2019; Wick et al., 2021).

Small plasmids

Small plasmids (a few kilobases) are often lost by long-read assemblers because size selection removes their short molecules, or because they are filtered out as low-length reads. If plasmids matter, check the short-read assembly or rerun with relaxed read-length filters.

Polishing and circularisation

A long-read assembly gets the structure right. Polishing fixes the remaining base-level errors, and circularisation gives the genome a clean, conventional start.

Why polish

Residual errors in long-read assemblies are mostly small insertions and deletions in homopolymers. Each indel in a coding region shifts the reading frame and splits a gene into fragments, which inflates gene counts and creates false pseudogenes. Polishing aligns reads back to the draft and corrects positions where the consensus disagrees with the draft.

• Long-read polishing with Medaka uses a neural network trained for a specific pore and basecaller model, so the model setting must match how the reads were basecalled. With R10.4.1 data and current basecallers, long-read-only genomes can reach very high accuracy (Sereika et al., 2022).
• Short-read polishing aligns Illumina reads (for example with BWA (Li & Durbin, 2009)) and fixes remaining errors. It must be careful in repeats, where short reads map ambiguously and can introduce new errors.

Circularisation and rotation

A circular replicon assembled as a linear contig may have overlapping ends. Circularisation trims the overlap, and rotation sets the start position, by convention at the dnaA gene for chromosomes and at repA for plasmids where possible. A consistent start makes genomes easier to compare and annotate.

Assembly quality: QUAST and CheckM

Two different questions decide whether an assembly is fit for use: is it contiguous, and is it complete and clean?

Contiguity with QUAST

QUAST reports the number of contigs, total length, largest contig, N50, L50 and GC content (Gurevich et al., 2013). N50 is the length such that contigs of that size or longer contain half the assembly. A finished bacterial genome has one contig per replicon and an N50 close to the chromosome size; short-read assemblies typically have N50 values of tens to hundreds of kilobases. Total length and GC content should match what is known for the species. With a reference genome, QUAST also counts misassemblies and mismatches.

Completeness and contamination with CheckM

CheckM uses sets of marker genes that are present in a single copy in nearly all genomes of a lineage (Parks et al., 2015). Missing markers lower completeness; duplicated markers raise contamination. Lineage-specific marker sets, chosen by placing the genome in a reference tree, are more precise than universal sets. CheckM2 replaces fixed marker sets with a machine-learning model and works better for novel or reduced genomes (Chklovski et al., 2023).

These thresholds were defined for genomes from metagenomes (Bowers et al., 2017); an isolate genome should comfortably exceed the high-quality bar. High contamination in an isolate usually means a mixed culture.

Genome taxonomy: ANI and GTDB

With a whole genome in hand, species identity is decided by genome-wide similarity, not by a single marker gene.

Average nucleotide identity

Average nucleotide identity (ANI) compares two genomes over their shared regions. Across many thousands of genome pairs, ANI values show a clear gap around 95 percent, which matches traditional species boundaries well (Jain et al., 2018). Two genomes above about 95 percent ANI are generally the same species. Fast estimators such as Mash use k-mer sketches to screen large collections before an exact calculation (Ondov et al., 2016).

Christensen (ed.), Introduction to Bioinformatics in Microbiology, ch. 7, explains the move from 16S identity and DNA–DNA hybridisation to genome-based species definitions; ch. 6 covers the phylogenetic methods behind GTDB.

The Genome Taxonomy Database

GTDB is a genome-based taxonomy built from a phylogeny of conserved single-copy proteins and normalised so that ranks correspond to comparable evolutionary depth (Parks et al., 2022). It reorganised many groups whose traditional names did not reflect their phylogeny, so GTDB names sometimes differ from those in older literature and from NCBI. GTDB-Tk places a genome in the GTDB reference tree and assigns species by ANI to the closest reference (Chaumeil et al., 2022). A genome with no species-level match is reported to genus, which may indicate a new species.

Genome annotation

Annotation finds the genes and other features in a genome and gives them names. The quality of every downstream biological claim depends on it.

Structural and functional annotation

Structural annotation locates features: protein-coding sequences (CDS), predicted by Prodigal in most bacterial pipelines, plus rRNAs, tRNAs, tmRNA, non-coding RNAs, CRISPR arrays and origins of replication. Functional annotation assigns a product name and identifiers to each CDS by comparison with protein databases. Prokka made this fast and standard for a decade (Seemann, 2014). Bakta uses a large precomputed database of UniProt protein clusters, identifies known proteins by exact hashing, and returns database cross-references for a high share of genes (Schwengers et al., 2021).

Bacterial gene finders such as Prodigal score open reading frames with the codon and start-site statistics of the genome itself, a simple form of the hidden Markov models that eukaryotic gene finders need (chapter 33). Baxevanis, Bader & Wishart, Bioinformatics, ch. 5, and Zvelebil & Baum, Understanding Bioinformatics, ch. 9–10, explain both, including how accuracy is measured. Christensen (ed.), Introduction to Bioinformatics in Microbiology, ch. 2–3, describe the databases used for functional annotation and how far a function can be transferred from a distant homolog.

Reading the result

A typical bacterial genome has about one gene per kilobase. Counts far above that often mean frameshift errors fragmenting genes (chapter 16); counts far below suggest a fragmented assembly. Missing rRNA operons in a short-read assembly are usually an assembly artefact, because identical operon copies collapse. A large fraction of "hypothetical protein" annotations is normal, especially for less studied genera.

Typing, plasmids and resistance

Once a genome is annotated, the practical questions are usually about lineage, mobile elements and genes of concern.

Multilocus sequence typing

MLST assigns an allele number to each of about seven housekeeping genes; the combination defines a sequence type (ST). STs are stable, comparable across laboratories and widely used in epidemiology. Schemes and allele definitions are curated in PubMLST (Jolley, Bray & Maiden, 2018). An "untyped" result usually means no scheme exists for the species or a new allele is present.

MLST assigns an allele number to each of seven housekeeping gene sequences and a sequence type (ST) to each combination of alleles: two isolates share an ST only if all seven alleles are identical. Core-genome MLST does the same with hundreds to thousands of genes and resolves outbreaks. Christensen (ed.), Introduction to Bioinformatics in Microbiology, ch. 11, traces typing from MLST to whole-genome methods.

Plasmids and viruses

Plasmids and prophages carry many of the genes that make strains different, including resistance and virulence genes. In a complete assembly, small circular contigs are plasmid candidates. geNomad classifies sequences as chromosome, plasmid or virus using marker genes and a neural network, and annotates them (Camargo et al., 2024). It works on whole genomes and on metagenome assemblies.

Antimicrobial resistance genes

Resistance genes are found by comparison with curated databases such as NCBI's AMRFinderPlus database or CARD (Feldgarden et al., 2021). A gene's presence predicts potential resistance but not the phenotype: expression, point mutations and the genetic context matter. Report the database version and identity thresholds, and confirm important findings with phenotypic testing.

Variant calling and comparative genomics

Comparing a strain with its parent or with close relatives reveals the mutations behind an evolved trait, an outbreak link or a lab adaptation.

Small variants

Reads from the strain of interest are aligned to a reference genome, ideally a complete genome of the parent strain, with BWA for short reads or minimap2 for long reads (Li & Durbin, 2009; Li, 2018). A variant caller then reports single-nucleotide polymorphisms (SNPs) and short indels with their read support. FreeBayes is a widely used haplotype-based caller (Garrison & Marth, 2012). For haploid bacteria, true variants should be supported by nearly all reads at a position; variants supported by a fraction of reads point to mixed populations, mapping errors in repeats or contamination.

Variant calling is two problems: mapping each read to the right place, then deciding whether a difference is real. Ismail, Bioinformatics: A Practical Guide to NGS Data Analysis, ch. 2 and 4, cover both, including base and mapping qualities, duplicate marking and the statistics of genotype calling; Buffalo, Bioinformatics Data Skills, ch. 11, explains the SAM/BAM format you will meet at every step.

Structural variants

Insertions, deletions, inversions and duplications larger than about 50 bases are detected from split and discordant alignments. Long reads make this much more reliable because a single read can span the whole event (Sedlazeck et al., 2018). Comparing two complete assemblies directly is an alternative that avoids reference bias.

Gene content

Across several genomes, gene presence and absence (the pangenome) shows which genes are shared and which are strain-specific. A heatmap of gene families across genomes, clustered by similarity, quickly shows groups of strains with similar accessory content.

Reproducibility and workflows

An analysis is reproducible when someone else, or you in a year, can rerun it and get the same result. That depends less on good intentions than on how the work is set up.

Pin software and databases

Results change with tool versions and, even more, with reference database versions. Record both for every step. Package managers such as conda with Bioconda channels make exact versions installable (Grüning et al., 2018); containers (Docker, Apptainer) freeze the whole environment. A one-line note such as "SILVA 138.1, DADA2 1.30, MetaPhlAn 4.1 with database vJun23" is worth more than a paragraph of prose.

Use a workflow manager

Workflow managers describe each step, its inputs and outputs, and its software environment, then run the steps in the right order, in parallel where possible, and resume after a failure without redoing finished work. Nextflow is widely used in genomics, including by the nf-core community, and runs the same workflow on a laptop, a cluster or the cloud (Di Tommaso et al., 2017). Snakemake follows a similar idea with Python-like rules.

Buffalo, Bioinformatics Data Skills is the best single book on this side of the work: project layout and documentation (ch. 2), version control with Git (ch. 5), Unix data tools (ch. 7) and writing robust pipelines (ch. 12). Its rule of thumb applies everywhere: anything you did by hand you will need to redo, so script it.

Keep data, code and results apart

Raw data should be read-only and stored once. Code lives in version control. Results are generated and can always be deleted and remade. Random seeds should be fixed for anything stochastic: rarefaction, permutation tests, cross-validation splits and random forests.

Share data properly

Deposit raw reads in an archive of the International Nucleotide Sequence Database Collaboration (ENA, SRA or DDBJ), with sample metadata that follows the MIxS checklists. Genomes go to the same archives with their annotation. The FAIR principles (findable, accessible, interoperable and reusable) summarise what reviewers and funders now expect (Wilkinson et al., 2016).

Reporting a microbiome or genome study

A good methods section lets a reader judge every choice that could have changed the result. A good results section separates what the data show from what they suggest.

Methods checklist for microbiome studies

• Sample collection, storage and time to freezing; DNA extraction kit and any bead-beating.
• Marker region and primers, or shotgun library type; sequencing platform and read length.
• Negative and positive controls and how they were used.
• Read processing and feature construction (tool and version), reference database and version.
• Filtering thresholds and how many features and reads each removed.
• Normalisation and transformation for each analysis; how depth was handled for alpha diversity.
• Distances, ordination and tests (with permutations and strata); differential abundance methods, covariates, random effects and FDR control.
• For prediction: model, cross-validation scheme (grouped, nested), metrics and validation cohort.

Methods checklist for genome studies

• Culture and DNA extraction; sequencing platform, chemistry and basecaller model.
• Read filtering; assembler, version and key settings; polishing tools and rounds.
• Assembly statistics (contigs, N50, size, GC) and completeness and contamination.
• Taxonomic method and database release; annotation tool and database version.
• Accession numbers for reads and assemblies.

Writing results

Give effect sizes with uncertainty: PERMANOVA R², fold changes with confidence intervals, AUC with its spread. State how many methods agreed on a differentially abundant taxon. Distinguish association from causation, and genetic potential from activity. Negative results reported with their power are useful to the field.

Report what was estimated, not only what was significant: effect sizes with intervals, the number of tests and the false discovery rate. When thousands of features are tested, a q-value of 0.05 means that about 5 % of the reported list is expected to be false; Efron & Hastie explain this and the related local false discovery rate in plain terms (Efron & Hastie, Computer Age Statistical Inference, ch. 15).

What RNA-seq measures, and how to design the study

RNA-seq turns the messenger RNA in a sample into millions of short sequences and counts them. Every later step, from normalisation to pathways, rests on two things decided before the first read exists: what RNA you captured, and how the samples were arranged.

From genes to reads

A gene is transcribed into RNA; only a small part of a cell's RNA is protein-coding messenger RNA (mRNA, about 2–4 %), while most is ribosomal RNA. In eukaryotes the primary transcript contains exons and introns; introns are spliced out (most begin with GT and end with AG), and one gene can produce several isoforms by alternative splicing. Bacterial genes have no introns and are often organised in operons, several genes transcribed together from one promoter. These differences shape everything that follows: eukaryotic reads must be aligned across splice junctions, and bacterial libraries cannot be built by catching poly(A) tails (Ismail, Bioinformatics: A Practical Guide to Next Generation Sequencing Data Analysis, ch. 5.1).

The library is made by enriching mRNA, reverse-transcribing it into complementary DNA (cDNA), fragmenting it, adding adapters and sequencing. Each read is a random fragment of a random transcript. A gene's count is therefore proportional to how many copies of its transcript were present multiplied by its length: a long transcript breaks into more fragments than a short one at the same expression.

Capturing the right RNA

• Poly(A) selection uses oligo-dT beads to keep mature eukaryotic mRNA. It is efficient and clean, but it misses non-polyadenylated RNA and performs poorly with partly degraded RNA.
• rRNA depletion removes ribosomal RNA instead. It is required for bacteria, whose mRNA has no poly(A) tail, and useful for degraded samples and for non-coding RNA (Christensen, Introduction to Bioinformatics in Microbiology, ch. 10).
• Stranded libraries (for example the dUTP method) record which DNA strand a transcript came from. This matters wherever genes overlap on opposite strands, which is common in compact bacterial and fungal genomes.
• RNA quality is summarised by the RNA integrity number (RIN). In good total RNA the large ribosomal band is about twice as intense as the small one; a smear below the bands signals degradation.

Replicates, depth and power

A biological replicate is an independent organism, culture or patient; a technical replicate is the same RNA prepared or sequenced twice. Differential expression asks whether a difference between groups is larger than the difference between independent individuals, so only biological replicates answer it. Three per group is a bare minimum; six or more is safer for modest effects, and beyond about 20 million reads per sample, more replicates add far more power than more depth (Christensen, ch. 10; Schurch et al. 2016).

Baxevanis and colleagues recommend a small pilot experiment: it gives a realistic estimate of variability, which makes a proper power calculation possible before the full study is funded (Baxevanis, Bader & Wishart, Bioinformatics, ch. 10).

Batches and confounding

A batch effect is any systematic difference introduced by handling: extraction day, kit lot, technician, sequencing run or lane. Batch effects are harmless when every batch contains every group, because the model can estimate and remove them. They are fatal when a batch coincides with a group (all controls on Monday, all treated on Tuesday): then no method can separate treatment from Monday (Baxevanis et al., Bioinformatics, ch. 10).

Transcriptome

All RNA molecules present in a sample at the time of collection.

Isoform

One of several transcripts produced from a gene by alternative splicing or start sites.

Library size

The number of reads assigned to genes in a sample; also called sequencing depth.

Biological replicate

An independent unit (organism, culture, patient) measured once.

Batch effect

A systematic technical difference between groups of samples processed together.

Stranded library

A library that keeps the information about which DNA strand was transcribed.

From reads to counts

The goal of the first half of an RNA-seq analysis is one table: genes in rows, samples in columns, a number in each cell. There are three ways to get there, and choosing between them is a question of what you need to know.

Quality control and trimming

Raw reads are first inspected with FastQC (one report per file) and summarised across samples with MultiQC. The usual problems are adapter sequence at read ends (the fragment was shorter than the read), falling base quality towards the 3′ end, over-represented sequences (often ribosomal RNA) and very high duplication. Tools such as fastp, Cutadapt or Trimmomatic remove adapters and low-quality tails; fastp is several times faster on large datasets (Low & Tammi, Practical Bioinformatics for Beginners, ch. 8). Duplication is normal in RNA-seq: a highly expressed gene produces many identical fragments, so duplicates are not removed as they are in genome sequencing.

Three routes to a count table

1. Splice-aware alignment to a genome

STAR and HISAT2 align reads to the genome and can split a read across an intron. STAR first finds the longest exact matches of each read ("seeds"), then stitches seeds that lie near each other into a spliced alignment (Ismail, ch. 5.3.2). Alignment gives positions, so it supports splicing analysis, variant calling and the discovery of unannotated genes; it needs a genome and an annotation (GTF) and, for STAR, a lot of memory.

2. Pseudo-alignment or selective alignment to a transcriptome

salmon and kallisto skip base-by-base alignment. They determine which transcripts are compatible with each read and then use an expectation-maximisation algorithm to share reads that fit several isoforms. They are fast, model multi-mapping reads and correct for fragment-level biases (GC content, sequence context). Transcript estimates are summed to genes for gene-level testing (tximport). For gene-level differential expression they are as accurate as alignment-based counting.

3. De novo transcriptome assembly

Without a reference, Trinity or rnaSPAdes assemble reads into transcripts. This is the route for non-model organisms; it is slower and produces redundant transcripts that need clustering (Low & Tammi, ch. 8).

Counting reads per gene

After alignment, featureCounts or HTSeq-count assign each read to a gene using the annotation. A read overlapping two genes is ambiguous and, in the default "union" mode, not counted; a read aligning to several places is multi-mapping and usually skipped (Ismail, ch. 5.3.4). With a stranded library, telling the counter the right strand setting matters: the wrong setting can silently lose most reads. Pseudo-aligners make this choice for you (salmon's -l A detects the library type).

Is the library usable?

• Mapping rate. Typically above 70–80 % to a good reference; far lower suggests contamination, rRNA, a wrong reference or a wrong strand setting.
• Uniquely versus multi-mapped reads and the number of splice junctions (STAR's Log.final.out).
• Gene-body coverage. Reads should cover transcripts evenly; a pile-up at the 3′ end means degraded RNA (Qualimap, RSeQC).
• Composition. Share of rRNA and mitochondrial reads: high values point to failed depletion or stressed cells.

Adapter

Synthetic sequence added to fragments; appears in a read when the fragment is shorter than the read length.

Splice-aware aligner

An aligner that can split a read across an intron (STAR, HISAT2).

Pseudo-alignment

Assigning a read to compatible transcripts without base-level alignment (kallisto, salmon).

Multi-mapping read

A read that fits several genomic locations or transcripts equally well.

GTF/GFF

Annotation files that list gene, transcript and exon coordinates.

Normalisation and exploring samples

Raw counts mix three things: how much a gene is expressed, how deeply the sample was sequenced, and what else was in the sample. Normalisation removes the last two so that the first can be compared. Then, before any test, you look at the samples.

Units for comparing within a sample

With k reads on a gene of length l in a library of N reads (Ismail, ch. 5.3.5):

FPKM is RPKM counted in fragments (read pairs). TPM divides by length first and then rescales, so TPM always sums to one million in every sample; RPKM does not, which makes RPKM values from different samples subtly incomparable.

Comparing between samples: composition bias

Scaling by library size assumes the total amount of RNA is shared out the same way in every sample. When a few genes become very highly expressed in one condition, they take reads from every other gene, and simple scaling makes unchanged genes look down-regulated. Two robust methods assume instead that most genes do not change:

• TMM (edgeR) compares each sample with a reference sample, computes log-ratios M for all genes, trims the most extreme 30 % of M and 5 % of the average abundances, and takes a precision-weighted mean of what remains (Ismail, ch. 5.3.5.4; Robinson & Oshlack 2010).
• Median of ratios (DESeq2) divides each count by the gene's geometric mean across samples and takes the median of these ratios per sample (Ismail, ch. 5.3.5.5).

Filtering genes before testing

Genes with almost no reads carry no information and add to the multiple-testing burden. edgeR's filterByExpr keeps genes with a minimum count-per-million in at least as many samples as the smallest group (Ismail, ch. 5.3.7.4).

A scale for plots

Count variance grows with the mean, so on the raw scale a few highly expressed genes dominate every distance, and on a naive log scale the noisy low counts dominate. The variance-stabilising transformation (VST) or regularised log of DESeq2, and log-CPM with a prior count in edgeR, give roughly constant variance across the range. Use them for PCA, clustering and heatmaps, never for the test itself.

Looking at the samples

PCA on the 500 most variable genes is the first plot to make. Colour it by condition, then by every technical variable: if samples separate by extraction date or sequencing run, a batch effect is present (Baxevanis et al., Bioinformatics, ch. 10). Hierarchical clustering joins the two closest samples, recomputes distances, and repeats until one tree remains; single, complete and average linkage differ only in how the distance between groups is defined (Baxevanis et al., ch. 10; Zvelebil & Baum, Understanding Bioinformatics, ch. 16). A heatmap of row-scaled expression with its dendrograms shows which genes drive the groups.

Sanity checks catch errors no statistic can: expression of Y-chromosome genes (or XIST) must match the recorded sex. Baxevanis and colleagues describe a dataset in which samples clustered by sex rather than by tissue, revealing mislabelled samples.

Checking sample identity: the sex check

Sample swaps happen in every large study. A cheap test uses genes whose expression depends on sex: XIST is expressed from the inactive X in females, and genes on the Y chromosome (RPS4Y1, DDX3Y, KDM5D, UTY, EIF1AY) only in males. Plotting one against the other separates the samples into two clouds; a sample in the wrong cloud for its recorded sex was probably swapped or mislabelled (Baxevanis, Bader & Wishart, Bioinformatics, ch. 10). RNA-realm draws this plot in step 02 for human and mouse data and flags disagreements with a sex column in the metadata. The demo study contains one deliberately swapped label.

Differential expression, step by step

Differential expression asks, for every gene at once, whether the difference between groups is larger than the difference between replicates. The statistics look complicated, but they rest on four ideas: counts are over-dispersed, genes can borrow strength from each other, the design must be written down, and thousands of tests need a different kind of error control.

From binomial to negative binomial

Each read either comes from gene g or not, so the count of a gene in one library is binomial; because any single gene takes a tiny share of millions of reads, this becomes a Poisson distribution, whose variance equals its mean (Ismail, ch. 5.3.6). The Poisson describes technical sampling only. Independent biological replicates vary more, so the variance is modelled as

where φ is the dispersion. Its square root is the biological coefficient of variation (BCV): φ = 0.1 means replicates differ by about 32 % (Ismail, ch. 5.3.7; Xia, Sun & Chen, Statistical Analysis of Microbiome Data with R, ch. 11).

Borrowing strength between genes

With three replicates per group, a gene's own variance estimate is very noisy. The solution is empirical Bayes shrinkage: fit a trend of dispersion against mean expression across all genes, then pull each gene's estimate part-way towards the trend; genes with less information are pulled further. edgeR computes a common, a trended and a gene-wise ("tagwise") dispersion; DESeq2 and limma use the same principle (Ismail, ch. 5.3.7.6).

Efron and Hastie explain why this works with a classic example: predicting baseball players' season batting averages from their first 90 at-bats. Shrinking every player's early average towards the overall mean (the James–Stein estimator) reduces the total error substantially, because each player borrows indirect evidence from the others. They also show the cost: a few genuinely extreme players are pulled too far (Efron & Hastie, Computer Age Statistical Inference, ch. 7). The same is true for genes, which is why shrunken fold changes (apeGLM in DESeq2) are best for ranking, while a gene with a large, well-supported effect keeps most of it.

Writing down the design

Each gene is fitted with a generalised linear model on the log scale, with the library size as an offset (Ismail, ch. 5.3.6; Faraway, Extending the Linear Model with R, ch. 5):

Three families of methods

• DESeq2: negative binomial GLM, empirical Bayes dispersions, Wald test (or likelihood-ratio test for several levels), automatic handling of outlier counts (Cook's distance) and independent filtering.
• edgeR: negative binomial GLM; the quasi-likelihood F-test (glmQLFit) adds a second, gene-specific variability and is recommended for most designs (Ismail, ch. 5.3.7.8).
• limma-voom: transforms counts to log-CPM, estimates the mean–variance trend and gives each observation a precision weight, then uses limma's moderated t-statistics, originally developed for microarrays (Gentleman et al., Bioinformatics and Computational Biology Solutions Using R and Bioconductor, ch. 23).

All three control false discoveries well in designed experiments with few replicates. In large human cohorts (dozens of samples per group) their distributional assumptions can fail; the Wilcoxon rank-sum test then controls the false discovery rate better (Li et al. 2022).

Twenty thousand tests at once

At p < 0.05, testing 20,000 unchanged genes yields about 1,000 false positives. Controlling the family-wise error rate (Bonferroni) prevents even one, but is too strict to find anything. The false discovery rate (FDR) controls instead the expected proportion of false discoveries among the genes you report (Efron & Hastie, ch. 15).

The Benjamini–Hochberg recipe: order the N p-values, find the largest i with p_(i) ≤ i·q/N, and reject hypotheses 1 to i.

Reading the results

• Volcano plot: log fold change against −log₁₀ p. Top corners are large, well-supported changes.
• MA plot: log fold change against mean expression. Before shrinkage the low-count genes fan out; after shrinkage they collapse towards zero (Ismail, ch. 5.3.8).
• p-value histogram: a flat floor with a peak near zero is healthy. A hill in the middle or a peak near one points to a missing covariate or a wrong model.
• Look at single genes. Plot normalised counts per sample for the top genes; a "hit" driven by one sample is a warning.

Beyond the adjusted p-value: the local false discovery rate

The Benjamini–Hochberg q-value controls the false discovery rate of the whole list above a threshold. Efron's local false discovery rate answers a sharper question for each gene: what is the probability that this particular gene is null? It needs two densities: f(z), the distribution of all genes' z-values, and f₀(z), the distribution of null genes. Then lfdr(z) = π₀ f₀(z) / f(z), where π₀ is the share of null genes (Efron & Hastie, Computer Age Statistical Inference, ch. 15).

Efron's key observation is that the null need not be the theoretical N(0, 1). Correlation between genes, hidden batches or unmodelled covariates widen or shift the bulk of the z-values; estimating the null from that bulk ("empirical null", by fitting a parabola to the log density in the middle, called central matching) gives more honest calls.

RNA-realm estimates the empirical null in step 05, draws the histogram with f(z) and π₀f₀(z), and adds an lfdr column to the results table.

Dispersion (φ)

The extra-Poisson variance of a gene between biological replicates; √φ is the biological coefficient of variation.

Empirical Bayes

Estimating a prior from all genes and using it to stabilise each gene's estimate.

Design matrix

The table of 0/1 and numeric columns that encodes groups and covariates in the model.

Contrast

The comparison you test, for example treated versus control.

FDR

The expected proportion of false discoveries among reported genes.

LFC shrinkage

Pulling noisy log fold changes of low-information genes towards zero.

From genes to biology: pathways and co-expression

A list of 400 differentially expressed genes is not yet a result. Gene-set methods ask which processes those genes belong to; co-expression methods ask which genes rise and fall together. Both turn individual genes into biology.

Over-representation analysis (ORA)

ORA asks whether a gene set (for example "interferon response") contains more of your significant genes than chance would give. The probability comes from the hypergeometric distribution, as when drawing coloured balls from an urn without replacement.

The background must be the genes that could have been called significant (those tested), not the whole genome: lowly expressed tissue-specific genes would otherwise make every expressed pathway look enriched (Timmons et al. 2015).

Gene set enrichment analysis (GSEA)

ORA throws away everything below the significance cut-off. GSEA uses the whole list ranked by the test statistic. Walking down the list, a running score rises at each gene in the set and falls otherwise; the maximum deviation is the enrichment score. A set whose genes cluster at the top (or bottom) of the ranking is enriched even if no single gene is significant, which is how coordinated small shifts in, say, mitochondrial genes are found. Significance comes from permuting gene sets of the same size.

Gene-set libraries

Start with the 50 MSigDB Hallmark sets: curated, largely non-redundant summaries of well-defined processes. Then use GO Biological Process or Reactome for detail. Gene Ontology terms overlap heavily, so dozens of related "significant" terms often describe one signal. KEGG pathways are free for academic use only (Baxevanis et al., Bioinformatics, ch. 13 Biological Networks and Pathways).

Clustering genes by their profiles

Genes that respond together are often regulated together. Compeau and Pevzner tell this story through yeast: when Saccharomyces cerevisiae runs out of glucose it inverts its metabolism and starts consuming the ethanol it produced, the diauxic shift, and hundreds of genes change expression in coordinated waves. Clustering their expression profiles over time recovers groups of genes with shared functions (Compeau & Pevzner, Bioinformatics Algorithms: An Active Learning Approach, ch. 8 How Did Yeast Become a Wine Maker?).

• k-means chooses k centres and alternates between assigning each gene to its nearest centre and moving each centre to the mean of its genes (Lloyd's algorithm). It needs k in advance and finds roughly spherical clusters.
• Hierarchical clustering repeatedly merges the closest genes or clusters into a tree; cutting the tree gives clusters at any resolution.
• Distance matters more than the algorithm: Pearson correlation (1 − r) groups genes with the same shape of response regardless of level, which is usually what you want.

Co-expression networks (WGCNA)

WGCNA builds a network where every pair of genes is connected with a weight from their correlation, raised to a soft-threshold power β that keeps strong correlations and fades weak ones. Topological overlap then measures how many neighbours two genes share, and hierarchical clustering of 1 − overlap gives modules. Each module is summarised by its eigengene (its first principal component), which is correlated with sample traits; genes most correlated with the eigengene are the hubs.

Modules need many samples: with fewer than 15, correlations between genes are too noisy, and modules are hypotheses at best.

Gene set

A named list of genes sharing a function, location or pathway.

ORA

Over-representation analysis: a hypergeometric test of significant genes in a set.

GSEA

Gene set enrichment analysis on the full ranked gene list.

Leading edge

The genes of a set that drive its GSEA enrichment score.

Module eigengene

The first principal component of a co-expression module's expression.

Guide: an RNA-seq study in RNA-realm

A walk through a complete bulk RNA-seq analysis in RNA-realm, from FASTQ files or a count table to a methods section, with the decision that matters at each step and the chapter that explains it.

Step 0 · Reads → counts (optional)

Upload paired (_R1/_R2) or single-end FASTQ files; they go up in pieces, so tens of gigabytes are fine. Choose the GENCODE human or mouse transcriptome, or upload your organism's transcript FASTA. RNA-realm trims with fastp, quantifies with salmon (selective alignment, GC and sequence bias correction) and sums transcripts to genes. Open the MultiQC report: mapping rate, adapter content and duplication tell you whether each library is usable (chapter 25). Then load the results with a sample sheet.

Step 1 · Import counts

Or start from counts: one table with all samples, or one file per sample from featureCounts, STAR, HTSeq, salmon, kallisto or RSEM. The first column of the sample sheet must contain the sample names exactly as in the count file; every other column is a variable. Human and mouse genes get symbols, gene types and chromosomes automatically.

Step 2 · Quality control

Check library sizes, the gene-type composition, the mitochondrial share and the sample correlation heatmap. Flagged outliers are not removed: look at them in the PCA first (chapter 26).

Step 3 · Filter and normalise

Choose the grouping variable for filterByExpr; keep TMM as a robust default, or median of ratios for DESeq2 size factors. Keep VST for plots.

Step 4 · Explore

Colour the PCA by condition and then by every technical variable. The "drivers" heatmap gives the R² of each variable on each component; a technical variable with high R² on PC1 or PC2 belongs in the model. "Remove a batch for display" changes only plots.

Step 5 · Differential expression

Choose the variable, the reference level and the test level; tick covariates. Run DESeq2 and limma-voom together; genes found by both form the consensus. With eight or more independent samples per group, add the Wilcoxon test. Read the p-value histogram before the gene list, and plot single genes (chapter 27).

Step 6 · Pathways

Start with the Hallmark collection, then GO or Reactome. Click a term for its running enrichment score and leading-edge genes (chapter 28).

Step 7 · Co-expression

With fifteen or more samples, find modules, look at which module follows your condition and which follows batch, and read each module's hub genes and best Hallmark term.

Step 8 · Report

The report writes a methods paragraph for exactly the steps you ran, with linked references. The Download menu exports every table.

Fungal biology for genome analysts

Fungal genomes make sense only in the light of how fungi live: they digest food outside their cells, grow as networks of hyphae that may hold many nuclei, switch between haploid and diploid stages, and make chemical weapons. Each of these habits leaves a mark that a genome analyst has to recognise.

What fungi are

Fungi include moulds, yeasts and mushrooms. About 100,000 species have been described, and estimates of the real number reach 1.5 million or more; environmental sequencing suggests that over 90 % of fungal species remain undiscovered. Fungi are the closest relatives of animals among eukaryotes (Madigan et al., Brock Biology of Microorganisms, ch. 18).

Fungi are chemoorganotrophs that feed by secreting enzymes: they digest polymers such as cellulose, chitin and proteins outside the cell and absorb the monomers. Decomposition of wood is a hallmark: white-rot basidiomycetes degrade both cellulose and lignin, while brown-rot fungi take the cellulose and leave the lignin (Brock, ch. 18.8). In a genome this lifestyle appears as large families of carbohydrate-active enzymes (CAZymes) and secreted proteins.

The major groups

• Microsporidia: tiny obligate parasites with drastically reduced genomes; Encephalitozoon has only 2.9 Mb and about 2,000 genes, fewer than E. coli, and lacks genes for whole metabolic pathways.
• Chytridiomycota: the earliest-diverging fungi, with motile, flagellated zoospores; Batrachochytrium dendrobatidis is driving amphibian declines.
• Zygomycota / Mucoromycota and Glomeromycota: coenocytic hyphae (bread mould Rhizopus); the arbuscular mycorrhizal fungi that colonise most plant roots.
• Ascomycota: the largest group, including yeasts (Saccharomyces, Candida) and moulds (Aspergillus, Penicillium, Fusarium); sexual spores form in a sac (ascus).
• Basidiomycota: mushrooms, rusts, smuts and the pathogen Cryptococcus; sexual spores form on a club-shaped basidium.

(Brock, ch. 18.9–18.13)

Hyphae, nuclei and ploidy: why they matter for sequencing

Most fungi grow as hyphae, tubular cells that extend at the tip and form a mycelium. Hyphae can be divided by cross-walls (septate) or be one continuous cell with hundreds of nuclei (coenocytic). In many basidiomycetes, two compatible haploid mycelia fuse and the resulting dikaryon keeps two genetically different nuclei in every cell for most of its life (Brock, ch. 18.13). For a genome project this means:

• DNA from a dikaryon or a diploid yeast contains two haplotypes, which appear as heterozygosity and can make an assembler produce two copies of the same region;
• DNA from a coenocytic or multinucleate isolate may mix slightly different nuclei;
• single-spore or single-cell isolation is the best way to obtain a clean haploid genome.

Mating types: a built-in repeat

In baker's yeast, the active MAT locus determines mating type a or α. Two silent copies of the mating-type genes sit elsewhere in the genome, and a haploid cell can switch type by copying a silent cassette into MAT (Brock, ch. 18.12). Three near-identical copies of the same sequence are exactly the kind of repeat that breaks short-read assemblies. Mating-type loci are also of interest in their own right: their structure (one locus or two, two alleles or many) defines whether a species is heterothallic or homothallic.

What fungal genomes look like

Genome history matters too. Compeau and Pevzner recount how the ancestor of Saccharomyces underwent a whole-genome duplication; most duplicates were later lost, leaving about 13 % of genes in pairs, such as the alcohol dehydrogenases ADH1 and ADH2 that let yeast make ethanol and then consume it (Compeau & Pevzner, Bioinformatics Algorithms, ch. 8). Such ancient duplicates explain why even a perfect yeast assembly shows a few "duplicated" genes in BUSCO (chapter 32).

Chemistry: secondary metabolites

Fungi produce antibiotics, toxins and pigments: penicillin from Penicillium, statins, aflatoxins from Aspergillus. The genes for one compound are usually clustered in the genome as a biosynthetic gene cluster, which is why tools such as antiSMASH search for clusters rather than single genes (chapter 34).

Hypha / mycelium

A tubular fungal filament; the network of hyphae is the mycelium.

Coenocytic

Hyphae without cross-walls, with many nuclei in a shared cytoplasm.

Dikaryon

A cell or mycelium with two genetically different haploid nuclei.

MAT locus

The genomic locus that determines mating type.

ITS

Internal transcribed spacers between ribosomal RNA genes; the official fungal barcode.

CAZyme

Carbohydrate-active enzyme; families of enzymes that build or break down sugars.

Profiling and assembling a fungal genome

Before assembling, count k-mers: they tell you the genome size, how heterozygous it is and how much of it repeats. Then choose an assembler that matches your reads, and judge the result with numbers you understand.

Reading a k-mer spectrum

Cut every read into overlapping words of length k (k = 21 is typical) and count how often each distinct k-mer occurs. Plot the number of distinct k-mers against their count, and the genome speaks:

• a spike at count 1–3 is sequencing errors, which create k-mers that occur once;
• the main peak sits at the k-mer coverage of the genome;
• a second peak at half that coverage means a heterozygous diploid (or a dikaryon): k-mers that differ between the two haplotypes are covered half as often;
• a long tail at high counts is repeats.

Three ways to assemble

Ismail groups assembly algorithms into three families (Ismail, Bioinformatics: A Practical Guide to Next Generation Sequencing Data Analysis, ch. 3):

• Greedy: repeatedly merge the two reads with the best overlap. Simple, but a repeat can lead it down the wrong path for good.
• Overlap–layout–consensus (OLC): reads are nodes, overlaps are edges, and the assembly is a path visiting every read once (a Hamiltonian path). Used for long reads (Canu; Flye uses a related repeat graph).
• de Bruijn graph: reads are cut into k-mers; each (k−1)-mer is a node and each k-mer an edge, and the assembly is a path using every edge once (an Eulerian path), which is computationally easy. Used for short reads (SPAdes, ABySS).

Long reads, short reads and polishing

Long reads (Oxford Nanopore, PacBio HiFi) span most fungal repeats, including transposons and rDNA units, and routinely give chromosome-scale contigs; Flye and hifiasm are standard choices. Short reads alone (SPAdes) give hundreds of contigs but accurate bases. A hybrid strategy assembles long reads and then polishes the consensus with short reads (Pilon, Polypolish), correcting the small insertions and deletions that remain in long-read assemblies, especially in homopolymers (Low & Tammi, Practical Bioinformatics for Beginners, ch. 6). Polishing should converge: the number of changes should fall sharply from one round to the next.

Contiguity: N50 and L50

Sort contigs from longest to shortest and add their lengths until half of the total is reached. The length of the contig that crosses the half-way mark is the N50; the number of contigs needed is the L50 (Ismail, ch. 3).

k-mer

A substring of length k from a read.

Coverage (depth)

The average number of reads covering each base of the genome.

Contig / scaffold

A contiguous assembled sequence; scaffolds join contigs with gaps of estimated size.

de Bruijn graph

A graph of overlapping k-mers whose Eulerian path spells the genome.

Polishing

Correcting residual base errors in an assembly using aligned reads.

N50 / L50

The contig length, and count, at which half the assembly is reached.

Is the assembly complete, clean and correctly named?

Three questions decide whether an assembly can be trusted: does it contain the whole genome, does it contain only this genome, and is the organism what you think it is?

Completeness with BUSCO

BUSCO searches the assembly for genes that are present in a single copy in almost every member of a lineage (for example 758 genes for all fungi, 2,137 for Saccharomycetes). The narrower the lineage, the more genes and the more sensitive the test (Low & Tammi, Practical Bioinformatics for Beginners, ch. 6). Each gene is reported as:

• Complete, single copy (S): found once, full length. This should be the large majority.
• Complete, duplicated (D): found more than once. A few percent can be real (yeast's ancient whole-genome duplication, chapter 30); many points to an unpurged diploid assembly or contamination with a relative.
• Fragmented (F): only part found, often because of assembly errors or genes split across contigs.
• Missing (M): not found at all.

Run BUSCO twice: on the genome (in genome mode, using the miniprot aligner) and on the predicted proteins. A gap between the two (complete in the genome, missing among proteins) points to gene-prediction problems, not a missing genome.

Telomeres: counting complete chromosomes

Most fungi end their chromosomes with the telomeric repeat TTAGGG (budding yeasts use irregular TG₁₋₃ repeats). A contig with telomeric repeats at both ends is a complete chromosome. Counting them is the simplest proof of a telomere-to-telomere assembly; a contig with a telomere at one end is a chromosome arm.

Contamination: the GC/depth (blob) plot

Plot every contig by its GC content against its read depth, with dot size proportional to length:

• the nuclear genome forms one cloud at the expected depth and GC;
• the mitochondrial genome sits apart, usually at much higher depth (many copies per cell) and lower GC;
• bacterial contaminants appear as a separate cloud at their own GC and depth;
• small contigs at half depth in a diploid can be haplotype duplicates.

Tiara classifies each contig as nuclear, mitochondrial, plastid or prokaryotic with a neural network trained on k-mer composition, which makes these clouds easy to name.

Identity: the ITS barcode

The internal transcribed spacers (ITS1 and ITS2) between the ribosomal RNA genes are the official DNA barcode of fungi: variable enough to separate species, flanked by conserved genes that are easy to find. ITSx extracts the region; a BLAST search against UNITE or NCBI's curated fungal ITS reference set confirms the species. Identification by ITS is as good as the reference database, and some genera (for example Fusarium, Aspergillus section-level) need secondary markers such as TEF1 or β-tubulin.

BUSCO

Benchmarking Universal Single-Copy Orthologs: a completeness test based on conserved genes.

Telomere-to-telomere

A chromosome assembled from one end to the other without gaps.

Blob plot

A GC-versus-coverage plot of contigs used to spot contamination.

Mitogenome

The mitochondrial genome; circular, high-copy, usually AT-rich in fungi.

Repeats and gene prediction

Finding genes in a eukaryotic genome means recognising a grammar: a start codon, exons, splice sites, introns, a stop codon. Repeats must be hidden first, or the predictor will happily call every transposon a gene.

Repeats and masking

Fungal genomes carry transposable elements (retrotransposons with long terminal repeats, DNA transposons), simple sequence repeats, rDNA arrays and segmental duplications. Repeat content ranges from a few percent in yeasts to more than half the genome in some plant pathogens, where effector genes often sit in repeat-rich, fast-evolving compartments. RepeatModeler2 builds a species-specific library of repeat families de novo; RepeatMasker then marks every copy. Soft-masking writes repeats in lower case (ACGTacgtacgtACGT) instead of replacing them with N, so gene predictors can avoid starting genes in repeats yet still extend a real gene across a short repeat.

The grammar of a gene

A eukaryotic protein-coding gene is a sequence of states: intergenic region, start codon (usually in a Kozak-like context), initial exon, donor splice site (GT), intron, acceptor splice site (AG), internal exons, terminal exon, stop codon (Baxevanis, Bader & Wishart, Bioinformatics, ch. 5). A hidden Markov model (HMM) represents each element as a state with its own sequence statistics (codon usage in exons, base composition in introns, position-specific signals at splice sites) and transition probabilities between states. Given a genome sequence, the Viterbi algorithm finds the most probable path through the states, which is the predicted gene structure (Zvelebil & Baum, Understanding Bioinformatics, ch. 10).

The parameters come from a training set of genes whose structure is known, which is why the species model matters: AUGUSTUS ships models for Saccharomyces, Aspergillus, Fusarium, Neurospora and others, and a model trained on a distant species misplaces introns. GENSCAN introduced this design for human genes; AUGUSTUS, GeneMark-ES (self-training) and SNAP follow it (Baxevanis et al., ch. 5).

Ab initio versus evidence-based

• Ab initio (intrinsic) predictors use only the genome sequence and the model. They are fast and need no extra data, but miss untypical genes and misjudge untranslated regions.
• Evidence-based (extrinsic) predictors add RNA-seq alignments (which show exactly where introns are), protein alignments from related species, or both. BRAKER3 and funannotate combine evidence with AUGUSTUS and GeneMark to produce markedly better gene models.

RNA-seq of the same strain is the single most valuable addition to a fungal annotation project; the RNA-realm chapters (24–29) describe how to produce it.

How good are the predictions?

Accuracy is measured against a curated reference at the nucleotide, exon or gene level, counting true positives (TP), false positives (FP) and false negatives (FN) (Baxevanis et al., ch. 5).

In practice, without a curated reference, check BUSCO in protein mode, the distribution of protein lengths (a spike of very short proteins suggests fragmented models) and the share of genes with InterPro domains.

Non-coding genes

tRNA genes are found with tRNAscan-SE, which combines a fast first-pass search with covariance models of the tRNA secondary structure; ribosomal RNA genes with barrnap or RNAmmer.

Codon usage: the statistics behind gene finding

The 61 sense codons encode 20 amino acids, so most amino acids have synonymous codons. Genomes do not use them equally: base composition (mutation bias) pushes towards GC- or AT-ending codons, and in highly expressed genes selection favours codons matching abundant tRNAs. Gene finders exploit this bias to tell coding from non-coding sequence (Zvelebil & Baum, Understanding Bioinformatics, ch. 10).

Wright's effective number of codons (ENC) summarises the bias of a gene in one number, from 20 (one codon per amino acid) to 61 (all synonymous codons equally used). Plotted against GC at synonymous third positions (GC3s), genes whose codon choice is explained by composition alone fall on the curve ENC = 2 + s + 29 / (s² + (1 − s)²); genes well below it are shaped by selection. Relative synonymous codon usage (RSCU) shows which codon is preferred for each amino acid.

FungFlow draws the ENC–GC3s plot and the RSCU table in the Genes results.

Soft-masking

Writing repeat sequence in lower case so that tools can treat it differently without losing it.

HMM

A probabilistic model of hidden states (exon, intron…) that emit observed bases.

Splice donor / acceptor

The 5′ (GT) and 3′ (AG) ends of an intron.

Training set

Genes of known structure used to estimate a predictor's parameters.

Sensitivity / precision

The share of true features found / the share of predictions that are true.

What the genes do: functional annotation of fungi

A list of 10,000 predicted proteins becomes biology when you know which ones digest plant walls, which make toxins, which are secreted to manipulate a host, and how this compares with relatives.

Transferring function by orthology

Most proteins are annotated by similarity to proteins whose function is known. Orthologs (genes separated by speciation) usually keep their function; paralogs (separated by duplication) often diverge. eggNOG-mapper assigns each protein to a precomputed orthologous group and transfers its annotations: gene name, GO terms, EC numbers, KEGG pathways and COG functional category. Protein domains (Pfam, InterPro) add evidence at the level of building blocks (Baxevanis et al., Bioinformatics, ch. 7 Predictive Methods Using Protein Sequences).

CAZymes: reading a lifestyle from enzymes

Because fungi digest food outside their cells (chapter 30), the enzymes they secrete are a direct readout of what they eat. The CAZy database groups them into glycoside hydrolases (GH), glycosyltransferases (GT), polysaccharide lyases (PL), carbohydrate esterases (CE), auxiliary activities (AA, mostly oxidative enzymes, including those that attack lignin) and carbohydrate-binding modules (CBM). dbCAN annotates them with HMMs and sequence similarity.

Secondary metabolism: biosynthetic gene clusters

Genes for a fungal natural product usually lie together: a core enzyme (polyketide synthase, non-ribosomal peptide synthetase, terpene synthase), tailoring enzymes, a transporter and often a regulator. antiSMASH finds these clusters, predicts their class and compares them with clusters of known products. Penicillin, made by Penicillium, is the classic example (Brock, ch. 18.9). A typical filamentous ascomycete has 30–70 clusters, most of whose products are unknown; this "cryptic" chemistry is why fungal genomes are mined for new drugs.

Secretome and effectors

Secreted proteins carry an N-terminal signal peptide (predicted by SignalP or similar tools) and lack transmembrane domains. Plant- and animal-pathogenic fungi secrete effectors, often small, cysteine-rich proteins that manipulate the host; EffectorP predicts them from sequence features. Predictions are a shortlist for experiments, not proof of function.

Putting it in context

Annotation becomes meaningful by comparison: gene families expanded or lost relative to relatives (orthogroups from OrthoFinder), clusters unique to one strain, CAZyme profiles that set a pathogen apart from its saprotrophic sister species (Bergman, Comparative Genomics).

Ortholog / paralog

Homologous genes separated by speciation / by duplication.

COG category

A one-letter functional class (for example G, carbohydrate metabolism) from the COG/eggNOG system.

BGC

Biosynthetic gene cluster: co-located genes that make one natural product.

Secretome

The set of proteins an organism secretes.

Effector

A secreted protein that alters host cells to benefit the pathogen.

Guide: a fungal genome in FungFlow

How to go from sequencing reads of one fungal isolate to an assembled, checked and annotated genome in FungFlow, and how to read each part of the report, with the chapter that explains it.

1. Reads, or an assembly

In Reads, upload Illumina pairs, Nanopore or PacBio HiFi FASTQ files. If you already have an assembly, upload the FASTA in Assembly instead: FungFlow then starts at quality control. For a new isolate the best value today is Nanopore R10 (30–50×) plus Illumina (40×) for polishing, or HiFi alone (chapter 31).

2. Set up the run

In New run, mark each file (Illumina R1, R2, long reads, assembly). Give the expected genome size if you know it, the ploidy (haploid for most filamentous ascomycetes), the narrowest BUSCO lineage that contains your fungus, the AUGUSTUS gene model of the closest species, and the telomere repeat (TTAGGG for most fungi). De novo repeat modelling is slower but gives much better masking. Runs queue on the server and continue when you close the page; follow them in Runs.

3. Read the report

• Read quality: yield, long-read N50, Q30 of short reads (chapter 31).
• Genome profile: the k-mer estimate of size and heterozygosity, independent of assembly; compare it with the assembly size.
• Assembly: size, contigs, N50 and the cumulative length curve; Pilon rounds should drop to few changes.
• Completeness: BUSCO of the genome and of the proteins (chapter 32).
• Contigs & identity: the GC/depth plot and Tiara classes, contigs with telomeres at both ends, and the ITS barcode with a link to NCBI BLAST.
• Repeats, genes, function: masked fraction, gene and protein statistics, tRNAs, COG categories, CAZyme classes and biosynthetic gene clusters (chapters 33–34).

4. Before you publish

Remove contaminant contigs, check that BUSCO duplication is low, and consider RNA-seq evidence for gene models (BRAKER3 or funannotate) if genes are your focus. Submit the genome to NCBI or ENA with the run's methods text and the tool versions from the log.

The library: books and the state of the art

These docs are a map. This chapter points to the territory: the textbooks that teach each module properly, the order in which to read them, and the papers that define the current state of the art. Every entry is in the reference list at the end, with a link.

A learning path

Foundations shared by all modules

Christensen (ed.), Introduction to Bioinformatics in Microbiology (2018)

A course-sized textbook written for microbiologists: databases, alignment, phylogeny, identification, amplicon and shotgun metagenomics, transcriptomics and typing, each in one chapter with exercises. The best first book for MicroRes and BactFlow users.

Low & Tammi (eds.), Practical Bioinformatics for Beginners (2023)

Hands-on chapters from raw-read quality control to assembly, transcriptomics, metagenomics and machine learning, with commands to run.

Ismail, Bioinformatics: A Practical Guide to NGS Data Analysis (2023)

Each chapter explains the algorithm and then runs the tools: quality control, mapping, de novo assembly, variant discovery, RNA-seq and metagenomics.

Baxevanis, Bader & Wishart (eds.), Bioinformatics, 4th ed. (2020)

The standard reference for concepts: sequence similarity, genome annotation, expression analysis, networks and pathways.

Zvelebil & Baum, Understanding Bioinformatics (2008)

The theory behind the tools, with derivations: hidden Markov models, gene detection, clustering and statistics.

Compeau & Pevzner, Bioinformatics Algorithms (2015)

Algorithms taught through biological questions; volume II's chapter 8 uses yeast's switch to wine-making to teach clustering of expression data.

Buffalo, Bioinformatics Data Skills (2015)

The craft: Unix, file formats, Git, R and robust pipelines. Read it before your first large project.

Statistics

James, Witten, Hastie & Tibshirani, An Introduction to Statistical Learning (2nd ed., 2021)

Regression, classification, resampling, regularisation, trees, unsupervised learning and multiple testing, with labs; free from the authors' website.

Efron & Hastie, Computer Age Statistical Inference (2016)

How statistics changed with computing: the bootstrap, empirical Bayes, the missing-species problem, false discovery rates and the empirical null. Used in MicroRes (unseen taxa) and RNA-realm (local fdr).

Faraway, Extending the Linear Model with R (2nd ed., 2016)

Generalised linear and mixed models; chapter 5 on count regression explains over-dispersion and the negative binomial.

Gentleman et al. (eds.), Bioinformatics and Computational Biology Solutions Using R and Bioconductor (2005)

The origins of limma and of multiple-testing practice in genomics, from the people who wrote the software.

Lee (ed.), Statistical Bioinformatics (2010)

A collection of statistical methods for high-throughput biological data.

MicroRes: microbiome analysis

Books. Xia, Sun & Chen, Statistical Analysis of Microbiome Data with R (2018), is the statistician's text: power and sample size (ch. 5), diversity (ch. 6–9), compositional analysis (ch. 10) and models for over-dispersed and zero-inflated counts (ch. 11–12). Its successor, Xia & Sun, Bioinformatic and Statistical Analysis of Microbiome Data (2023), adds the path from raw reads (QIIME 2, DADA2, ASVs) to longitudinal mixed models. Zhong Wang's Introduction to Computational Metagenomics (2022) explains the algorithms of shotgun metagenomics; Mitra's Metagenomic Data Analysis (2023) gives protocols; Streit & Daniel's Metagenomics: Methods and Protocols and Izard & Rivera's Metagenomics for Microbiology cover the laboratory side; Charles, Liles & Sessitsch's Functional Metagenomics covers function-based screening; Paradis's Analysis of Phylogenetics and Evolution with R covers trees.

The state of the art. Microbiome data are compositional and must be analysed as such (Gloor et al. 2017). Exact sequence variants replaced OTUs (Callahan et al. 2016, 2017). For shotgun data, MetaPhlAn 4 extended marker-based profiling to unnamed species (Blanco-Míguez et al. 2023), and HUMAnN 3 profiles function (Beghini et al. 2021). Diversity should be compared at equal sample coverage (Chao & Jost 2012); whether to rarefy is still debated, with recent evidence in favour of repeated rarefaction (Schloss 2024). For differential abundance, a benchmark of 14 methods on 38 datasets found that they disagree widely, and recommended consensus across methods (Nearing et al. 2022); bias-corrected methods such as ANCOM-BC2 (Lin & Peddada 2024) and linear models such as MaAsLin 2 (Mallick et al. 2021) are current standards. Networks use compositionally aware estimators (Kurtz et al. 2015). The STORMS checklist sets the reporting standard (Mirzayi et al. 2021).

BactFlow: bacterial genomes

Books. Christensen's chapters on assembly and annotation (ch. 2), identification (ch. 7) and typing (ch. 11); Ismail on quality control, mapping, assembly and variants (ch. 1–4); Low & Tammi on comparing assemblies (ch. 6); Bergman's Comparative Genomics for whole-genome comparison.

The state of the art. With Nanopore R10.4 chemistry, long reads alone now give near-finished bacterial genomes (Sereika et al. 2022); consensus assembly from several assemblers improves reliability further (Wick et al. 2021), and short-read polishing that respects repeats remains useful (Wick & Holt 2022, Polypolish). Deep-learning variant callers outperform classical ones on bacterial Nanopore data (Hall et al. 2024). Annotation moved to fast, database-rich tools such as Bakta (Schwengers et al. 2021); completeness and contamination to machine learning (CheckM2; Chklovski et al. 2023); taxonomy to the Genome Taxonomy Database and GTDB-Tk (Parks et al. 2022; Chaumeil et al. 2022) with fast ANI (Jain et al. 2018; Shaw & Yu 2023). Pangenomes are built with error-aware tools such as Panaroo (Tonkin-Hill et al. 2020), and resistance genes found with curated databases such as AMRFinderPlus (Feldgarden et al. 2021).

RNA-realm: transcriptomics

Books. Ismail's RNA-seq chapter (ch. 5) covers quantification, normalisation and differential expression with worked commands; Christensen (ch. 10) and Low & Tammi (ch. 8) give shorter introductions; Baxevanis et al. (ch. 10) cover design, batch effects and exploration; Gentleman et al. explain limma and multiple testing; Efron & Hastie (ch. 15) the false discovery rate and empirical null; Compeau & Pevzner (vol. II, ch. 8) and James et al. (ch. 12) clustering.

The state of the art. A community survey still frames good practice (Conesa et al. 2016). Lightweight, bias-aware quantification (Salmon; Patro et al. 2017) imported at gene level (Soneson et al. 2015) is standard. DESeq2 (Love et al. 2014) with heavy-tailed fold-change shrinkage (Zhu et al. 2019) and limma-voom (Law et al. 2014) remain the reference methods; with large human cohorts, however, they can produce excess false positives and rank tests are safer (Li et al. 2022). At least six biological replicates per group are recommended (Schurch et al. 2016). Batch effects in counts can be adjusted with ComBat-seq (Zhang et al. 2020). Gene-set analysis uses fast pre-ranked GSEA (Korotkevich et al. 2021) on curated sets such as the Hallmarks (Liberzon et al. 2015), and co-expression networks follow WGCNA (Langfelder & Horvath 2008).

FungFlow: fungal genomes

Books. Brock Biology of Microorganisms (ch. 18) for fungal biology; Compeau & Pevzner for yeast genome duplication; Ismail (ch. 3) and Low & Tammi (ch. 6) for assembly; Baxevanis et al. (ch. 5) and Zvelebil & Baum (ch. 9–10) for gene prediction.

The state of the art. Fungal genomics is organised around large sequencing programmes such as the 1000 Fungal Genomes project and the MycoCosm portal (Grigoriev et al. 2014); Stajich (2017) reviews what genomes revealed about the kingdom. Genome profiling uses k-mers (GenomeScope 2; Ranallo-Benavidez et al. 2020); assembly uses long reads with Flye (Kolmogorov et al. 2019) or, for HiFi data, hifiasm (Cheng et al. 2021). Completeness is measured with BUSCO (Manni et al. 2021). Repeats are modelled with RepeatModeler2 (Flynn et al. 2020); gene prediction now combines RNA-seq and protein evidence (BRAKER3; Gabriel et al. 2024). Function is annotated with dbCAN3 for CAZymes (Zheng et al. 2023), antiSMASH 7 for secondary metabolism (Blin et al. 2023), SignalP 6 for secretion (Teufel et al. 2022) and EffectorP 3 for effectors (Sperschneider & Dodds 2022). The evolution of wood decay was reconstructed from CAZyme repertoires across 31 genomes (Floudas et al. 2012).

Keeping up

• Benchmarks beat claims: community challenges such as CAMI compare metagenomic tools on the same data (Meyer et al. 2022), and independent benchmarks such as Nearing et al. (2022) are worth more than a method's own paper.
• Read method papers' supplementary benchmarks and check the data they were tested on resemble yours.
• Pin versions: a method that is state of the art today will be updated next year; the methods text each Bioinformaniac tool writes records the versions.