nf-core/eager

Edit

Introduction

This document describes the output produced by the pipeline. Most of the plots are taken from the MultiQC report, which summarises results at the end of the pipeline.

The directories listed below will be created in the results directory after the pipeline has finished. All paths are relative to the top-level results directory.

Pipeline overview

The pipeline is built using Nextflow and processes data using the following steps:

• FastQC - Raw read QC
• MultiQC - Aggregate report describing results and QC from the whole pipeline
• Pipeline information - Report metrics generated during the workflow execution

FastQC

FastQC gives general quality metrics about your sequenced reads. It provides information about the quality score distribution across your reads, per base sequence content (%A/T/G/C), adapter contamination and overrepresented sequences. For further reading and documentation see the FastQC help pages.

FastQC results in the raw/ directory contain the results of running FastQC on the ‘raw’ input files (with no processing by the pipeline). The processed/ directory is for reads run through either fastp or AdapterRemoval preprocessing, and therefore will only be present if preprocessing is executed. These FastQC results occur in both html (the report) and .zip format (raw data).

MultiQC

MultiQC is a visualization tool that generates a single HTML report summarising all samples in your project. Most of the pipeline QC results are visualised in the report and further statistics are available in the report data directory.

Results generated by MultiQC collate pipeline QC from supported tools e.g. FastQC. The pipeline has special steps which also allow the software versions to be reported in the MultiQC output for future traceability. For more information about how to use MultiQC reports, see http://multiqc.info.

Pipeline information

Nextflow provides excellent functionality for generating various reports relevant to the running and execution of the pipeline. This will allow you to troubleshoot errors with the running of the pipeline, and also provide you with other information such as launch commands, run times and resource usage.

Reference Indexing

Depending on what is supplied by the user, and if --save_reference is supplied, this directory will contain various reference index files based on the reference FASTA file given to --fasta.

It is highly recommend to move these files to a central location or cache directory on your machine to facilitate resume of the indices across different pipeline runs. In many cases indexing the reference genome for alignment can be the longest step of a pipeline run, therefore re-using indices in future runs (supplied to the pipeline with flags such as --fasta_fai, --fasta_dict, etc. or added to the reference sheet provided to --fasta) can greatly speed up analyses on other samples.

Preprocessing

Falco

Falco is an emulation and drop-in replacement of the popular FastQC software to check large sequencing reads for common problems. It is technically more performant, and can be useful in cases where you have large data or issues with FastQC’s Java memory allocation. For further reading and documentation see the Falco GitHub repository.

Falco results in the raw/ directory contain the results of running Falco on the ‘raw’ input files (with no processing by the pipeline). The processed/ directory is for reads run through either fastp or AdapterRemoval preprocessing, and therefore will only be present if preprocessing is executed. These Falco results occur in both .html (the report) and .txt format (raw data), and should provide the same information as with FastQC.

AdapterRemoval

AdapterRemoval removes residual adapter sequences from single-end (SE) or paired-end (PE) FASTQ reads, optionally trimming Ns and low qualities bases and/or collapsing overlapping paired-end mates into one read.

The resulting FASTQ files will only be present in your results directory if you run with --preprocessing_savepreprocessedreads. The .settings log files should always be present if AdapterRemoval has been run.

fastp

fastp is an all-in-one FASTQ preprocessor (QC/adapters/trimming/filtering/splitting/merging etc.)

The resulting FASTQ files will only be present in your results directory if you run with --preprocessing_savepreprocessedreads. The various log files should always be present if fastp has been run.

Mapping

BWA

BWA is a software package for mapping low-divergent sequences against a large reference genome, such as the human genome. BWA-backtrack (a.k.a bwa aln) is designed for Illumina sequence reads up to 100bp, while BWA-MEM (bwa mem) is optimised for longer sequences ranged from 70bp to 1Mbp and split alignment.

Bowtie 2

Bowtie 2 is an ultrafast and memory-efficient tool for aligning sequencing reads to long reference sequences. It is particularly good at aligning reads of about 50 up to 100s of characters to relatively long (e.g. mammalian) genomes. Bowtie 2 indexes the genome with an FM Index (based on the Burrows-Wheeler Transform or BWT) to keep its memory footprint small and supports gapped, local, and paired-end alignment modes.

Host Removal

The nf-core/eager host removal step removes any reads from the original input FASTQ files that mapped to the provided reference. This step is performed by a custom python script extract_map_reads.py.

The resulting FASTQ files can be used for data sharing purposes where the host DNA should not be made publically available for ethical reasons (for example when you do not have permission to analyse the host DNA but only analysis on the metagenomic content were agreed) or for due to identifiability reasons (for example, for individuals that died in recent times and could be identified based on their DNA alone).

Alternatively you could use the resulting files for manual metagenomic screen outside of nf-core/eager (e.g. when your preferred taxonomic classifier isn’t supported by nf-core/eager).

Mapping statistics

EndorSpy

endorS.py calculates percent on target (aka Endogenous DNA) from samtools flagstat files and prints to screen. The ‘Percent on target (%)’ reported will be different depending on the combination of samtools flagstat inputs provided. This program also calculates clonality (aka ‘cluster factor’) and percent duplicates when the flagstat file after duplicate removal is provided.

Definitions and calculations:

We provide up to 3 different estimates of percent on target:

Percent on target (%): percent of mapping reads mapping to a specific reference. This is calculated as follows: number of reads mapping to the specific reference (Mapped_Raw) in comparison to the total number of reads (Total_Reads)

$PercentOnTarget = { MappedRaw \over TotalReads \* 100 }$

Where $Mapped_Raw$ is the number of reads mapping to the reference during the mapping step, and $Total_Reads$ are the total number of reads that went into the mapping step.

Percent on target modified (%): percent of mapping reads after filtering the raw bam based on quality, read length or any other filtering performed that mapped to a specific reference:

$PercentOnTargetModified = { MappedPostFiltering \over TotalReads \* 100 }$

Where $Mapped_Post_Filtering$ are the number of reads mapping to the reference after applying the filters.

Percent on target postdedup (%): percent of deduplicated reads (either raw mapped reads or filtered mapped reads) that mapped to a specific reference:

$PercentOnTargetPostdedup = { MappedPostDedup \over TotalReads \* 100 }$

Additionally, this script provides two different ways of estimating library complexity:

Clonality (Cluster Factor in eager1): ratio of reads that have duplicated reads

$Clonality = { TotalReadsPreDedup \over mappedDedup }$

Percent Duplicates (%): percent of mapping reads that have at least 1 duplicate. It is calculated as follows:

$PercentDuplicates = {(TotalReadsPreDedup - MappedReadsDedup) \over TotalReadsPreDedup \* 100}$

The combination of statistics calculated would vary depending on the steps taken: Mapping/bam input + filtering + deduplication (all): Percent on target (%) Percent on target modified (%) Percent on target postdedup (%) Clonality Percent Duplicates (%)

Mapping/bam input: Percent on target (%)

Mapping/bam input + filtering: Percent on target (%) Percent on target modified (%)

Mapping/bam input + deduplication: Percent on target (%) Percent on target postdedup (%) Clonality Percent Duplicates (%)

⚠️ Warning: When bam input is provided, please keep in mind it is assumed that this is an unfiltered bam. If you provide an already filtered bam, the percent on target calculations will be wrong since the original total number of reads in the initial fastq/bam cannot be inferred.

BAM Filtering

nf-core/eager can perform a range of different length, quality, and/or mapped-unmapped filtering of BAM files, and generate converted FASTQ files for manual analysis outside the pipeline.

If bam filtering is turned on, by default only mapped reads are retained for downstream genomic analysis. Additionally, please note that removal of PCR duplicates in nf-core/eager will additionally filter for mapped reads only.

Please be aware, that intermediate length and mapping quality filtered genomic BAM files are not saved in the results directory automatically, with the assumption that downstream deduplicated BAM files are preferred for further analyses - see the parameters page for more information if you wish to save these.

You may also receive the files above if metagenomic screening is turned on.

Metagenomics Screening

Bbduk

The entropy filter of BBDuk is used to remove reads from the fastq-files for metagenomics screening that don’t pass the --metagenomics_complexity_entropy threshold. When ommiting the flag, a default value of 0.3 is applied. To cite the docs:

A homopolymer such as AAAAAAAAAAAAAA would have entropy of zero; completely random sequence would have entropy approaching 1.

Using complexity-filtered fastq-files as input for metagenomic classifiers can reduce the number of false positive classifications, resulting in more precise taxonomic assignments of the sample through removal of reads that can align equally well to multiple reference genomes. Save the complexity-filtered fastq-files to the output directory to perform additional downstream analyses, such as testing multiple metagenomic profilers (see the nf-core/taxprofiler pipeline).

Note: To save output files, set the --metagenomics_complexity_savefastq flag

PRINSEQ++

PRINSEQ++ is an alternative to BBDuk for filtering the fastq files before metagenomics classification. From PRINSEQ++ we implemented filtering by the dust algorithm or by entropy, as explained above in the BBDuk section.

Using complexity-filtered fastq-files as input for metagenomic classifiers can reduce the number of false positive classifications, resulting in more precise taxonomic assignments of the sample. Save the complexity-filtered fastq-files to the output directory to perform additional downstream analyses, such as testing multiple metagenomic profilers (see the nf-core/taxprofiler pipeline).

The saved files are the good files, passing the dust or entropy filter treshold specified. The logs contain information about the amount of reads filtered.

Note: To save output files, set the --metagenomics_complexity_savefastq flag

Deduplication

Deduplication is carried by two possible tools, as described below. However the expected output files are should be the same.

picard MarkDuplicates

Picard is a toolkit for general BAM file manipulation with many different functions. nf-core/eager most visibly uses the ‘markduplicates’ tool, for the removal of exact PCR duplicates that can occur during library amplification and results in false inflated coverages (and overly-confident genotyping calls).

DeDup

DeDup is a merged-read deduplication tool capable of performing merged-read deduplication on paired-end sequencing data in BAM files.

Preseq

Preseq Preseq is a collection of tools that allow assessment of the complexity of the library, where complexity means the number of unique molecules in your library (i.e. not molecules with the exact same length and sequence). There are two algorithms from the tools we use: c_curve and lc_extrap. The former gives you the expected number of unique reads if you were to repeated sequencing but with fewer reads than your first sequencing run. The latter tries to extrapolate the decay in the number of unique reads you would get with re-sequencing but with more reads than your initial sequencing run.

The resulting histogram file will contain estimated deduplication statistics at different theoretical sequencing depths, and can be used to generate a complexity curve for estimating the amount unique reads that will be yield if the library is re-sequenced.

These curves will be displayed in the pipeline run’s MultiQC report, however you can also use this file for plotting yourself for further exploration e.g. in R to find your sequencing target depth.

Mapping Statistics

QualiMap

QualiMap is a tool which provides statistics on the quality of the mapping of your reads to your reference genome. It allows you to assess how well covered your reference genome is by your data, both in ‘fold’ depth (average number of times a given base on the reference is covered by a read) and ‘percentage’ (the percentage of all bases on the reference genome that is covered at a given fold depth). These outputs allow you to make decision if you have enough quality data for downstream applications like genotyping, and how to adjust the parameters for those tools accordingly.

NB: Neither fold coverage nor percent coverage on its own is sufficient to assess whether you have a high quality mapping. Abnormally high fold coverages of a smaller region such as highly conserved genes or un-removed-adapter-containing reference genomes can artificially inflate the mean coverage, yet a high percent coverage is not useful if all bases of the genome are covered at just 1x coverage.

Note that many of the statistics from this module are displayed in the General Stats table, as they represent single values that are not plottable.

You will receive output for each sample. This means you will statistics of deduplicated values of all types of libraries combined in a single value (i.e. non-UDG treated, full-UDG, paired-end, single-end all together).

⚠️ Warning: If your library has no reads mapping to the reference, this will result in an empty BAM file. Qualimap will therefore not produce any output even if a BAM exists!

Bedtools

bedtools utilities are a swiss-army knife of tools for a wide-range of genomics analysis tasks. Bedtools allows one to intersect, merge, count, complement, and shuffle genomic intervals from multiple files in widely-used genomic file formats such as BAM, BED, GFF/GTF, VCF.

The bedtools coverage tool computes both the depth and breadth of coverage of features in file B (alignment file) on the features in file A (provied by --mapstats_bedtools_featurefile when running the eager workflow). One advantage that bedtools coverage offers is that it not only counts the number of features that overlap an interval in file A, it also computes the fraction of bases in the interval in A that were overlapped by one or more features. Thus, bedtools coverage also computes the breadth of coverage observed for each interval in A.

The output from this module can be useful for things such as checking for the presence/absence of virulence factors in ancient pathogen genomes, or getting statistics on SNP capture positions.

Damage Manipulation

There are three different options for manipulation of ancient DNA damage.

Damage Rescaling

mapDamage2 can apply a probabilistic bayesian model to rescale quality scores of likely-damaged positions in the reads. A new BAM file is constructed by downscaling quality values for misincorporations likely due to ancient DNA damage according to their initial qualities, position in reads, and damage patterns. These BAM files have damage probabilistically removed via a bayesian model. Rescaling a BAM file in this way can help reduce/remove the effects of DNA damage from downstream analysis.

Be advised that this process introduces reference bias in the resulting rescaled BAMs, because only mismatches to the reference get corrected, while true mismatches that become reference alleles due to damage are not rescaled.

This functionality does not produce any MultiQC output.

⚠️ Warning: rescaled libraries will not be merged with non-scaled libraries of the same sample for downstream genotyping, as the underlying damage model may be different for each library. If you wish to merge these, please do this manually and re-run nf-core/eager using the merged BAMs as input.

Post-Mortem-Damage (PMD) Filtering

pmdtools implements a likelihood framework incorporating a postmortem damage (PMD) score, base quality scores and biological polymorphism to identify degraded DNA sequences that are unlikely to originate from modern contamination. Using the model, each sequence is assigned a PMD score, for which positive values indicate support for the sequence being genuinely ancient. By filtering a BAM file for damaged reads in this way, it is possible to ameliorate the effects of present-day contamination, and isolate sequences of likely ancient origin to be used downstream. By default, all positions are assessed for damage, but it is possible to provide a FASTA file masked for specific references (--damage_manipulation_pmdtools_masked_reference) or a BED file to mask the reference FASTA within nf-core/eager (--damage_manipulation_pmdtools_reference_mask). This can alleviate reference bias when running PMD filtering on capture data, where you might not want the allele of a SNP to be counted as damage when it is a transition.

Base Trimming

bamUtil trimBam trims the end of reads in a SAM/BAM file, changing read ends to N and quality to !, or by soft clipping (if command-line option, --clip is specified). By default reverse strands are reversed and then the left & right are trimmed. This means that --left actually trims from the right of the read in the SAM/BAM for reverse reads.

Within nf-core/eager, when BAM trimming is activated alongside PMD-filtering, trimming is performed on the PMD-filtered reads only.

Damage Estimation

DamageProfiler

DamageProfiler is a tool which calculates a variety of standard ‘aDNA’ metrics from a BAM file. The primary plots here are the misincorporation and length distribution plots. Ancient DNA undergoes depurination and hydrolysis, causing fragmentation of molecules into gradually shorter fragments, and cytosine to thymine deamination damage, that occur on the subsequent single-stranded overhangs at the ends of molecules.

mapDamage2

mapDamage2 is a computational framework written in Python and R, which tracks and quantifies DNA damage patterns among ancient DNA sequencing reads generated by Next-Generation Sequencing platforms. Based on a sample bam file and a fasta reference file, the software calculates the frequency of 5’ C to T and 3’ G to A misincorporations, as well as the read length distribution per strand.

Contamination estimation for Nuclear DNA

ANGSD nuclear contamination

ANGSD is a software for analyzing next generation sequencing data. Among other functions, ANGSD can estimate contamination for chromosomes for which one copy exists, i.e. X-chromosome for humans with karyotype XY. To do this, we first generate a binary count file for the X-chromosome (angsd) and then perform a Fisher’s exact test for finding a p-value and jackknife to get an estimate of contamination (contamination). Contamination is estimated with Method of Moments (MOM) and Maximum Likelihood (ML) for both Method1 and Method2. Method1 compares the total number of minor and major reads at SNP sites with the number of minor and major reads at adjacent sites, assuming independent errors between reads and sites, while Method2 only samples one read at each site to remove the previous assumption. The results of all methods for each library, as well as respective standard errors are summarised in nuclear_contamination.txt and nuclear_contamination_mqc.json.

Sex Determination

Sex.DetERRmine

Sex.DetERRmine calculates the coverage of your mapped reads on the X and Y chromosomes relative to the coverage on the autosomes (X-/Y-rate). This metric can be thought of as the number of copies of chromosomes X and Y that is found within each cell, relative to the autosomal copies. The number of autosomal copies is assumed to be two, meaning that an X-rate of 1.0 means there are two X chromosomes in each cell, while 0.5 means there is a single copy of the X chromosome per cell. Human females have two copies of the X chromosome and no Y chromosome (XX), while human males have one copy of each of the X and Y chromosomes (XY).

When a bedfile of specific sites is provided, Sex.DetERRmine runs much faster and additionally calculates error bars around each relative coverage estimate. For this estimate to be trustworthy, the sites included in the bedfile should be spaced apart enough that a single sequencing read cannot overlap multiple sites. Hence, when a bedfile has not been provided, this error should be ignored. When a suitable bedfile is provided, each observation of a covered site is independent, and the error around the coverage is equal to the binomial error estimate. This error is then propagated during the calculation of relative coverage for the X and Y chromosomes.

Note that in nf-core/eager this will be run on single- and double-stranded variants of the same library separately. This can also help assess for differential contamination between libraries.

Genotyping

pileupCaller

sequencetools pileupCaller is a tool to create genotype calls from bam files using read-sampling methods. It can call genotypes from samtools mpileup output, and is often used for ancient DNA. The output files are in Eigenstrat format, in which the gentype dataset is split across 3 different files (similar to PLINK). The .geno file contains the genotype table. Each number represents the genotype call of an individual at a specific position, with each column represents an individual, while each row represents a SNP. The .snp file contains the SNP annotation of the genotype table, and the .ind file contains the individual annotation of the genotype table. We provide an additional file named *_coverage.tsv which contains the number of covered SNPs and total number of SNPs for each individual in the eigenstrat dataset. One dataset is generated per reference genome, which includes genotypes from both single- and double-stranded libraries.

When using pileupCaller for genotyping, single-stranded and double-stranded libraries are genotyped separately. Single-stranded libraries are genotyped with the additional option --singeStrandMode, which ensure that deamination damage artefactts cannot affect the genotype calls, by only using the forward- or reverse-mapping reads when genotyping on transitions (depending on the alleles of the transition).

GATK UnifiedGenotyper

GATK’s UnifiedGenotyper uses a Bayesian genotype likelihood model to estimate simultaneously the most likely genotypes and allele frequency in a population of N samples, emitting a genotype for each sample. The system can either emit just the variant sites or complete genotypes (which includes homozygous reference calls) satisfying some phred-scaled confidence value. This tool has been deprecated by the GATK developers in favour of HaplotypeCaller, but is still cosidered a preferable genotyper for ancient DNA, given its ability to handle low coverage data. The output provided is a bgzipped VCF file for containing the genotype calls of each sample, it’s index file, as well as the statistics of the VCF file generated by bcftools stats.

GATK HaplotypeCaller

GATK’s HaplotypeCaller is capable of calling SNPs and indels simultaneously via local de-novo assembly of haplotypes in an active region. In addition, HaplotypeCaller is able to handle non-diploid organisms as well as pooled experiment data. This is the preferred genotyper for modern DNA. The output provided is a bgzipped VCF file for containing the genotype calls of each sample, it’s index file, as well as the statistics of the VCF file generated by bcftools stats.

FreeBayes

FreeBayes is a Bayesian genetic variant detector designed to find small polymorphisms, specifically SNPs (single-nucleotide polymorphisms), indels (insertions and deletions), MNPs (multi-nucleotide polymorphisms), and complex events (composite insertion and substitution events) smaller than the length of a short-read sequencing alignment. It calls variants based on the literal sequences of reads aligned to a particular target, not their precise alignment. This model is a straightforward generalization of previous ones (e.g. PolyBayes, samtools, GATK) which detect or report variants based on alignments. This method avoids one of the core problems with alignment-based variant detection - that identical sequences may have multiple possible alignments. The output provided is a bgzipped VCF file for containing the genotype calls of each sample, it’s index file, as well as the statistics of the VCF file generated by bcftools stats.

ANGSD

ANGSD is a software for analyzing next generation sequencing data. It can estimate genotype likelihoods and allele frequencies from next-generation sequencing data. The output provided is a bgzipped genotype likelihood file, containing likelihoods across all samples per reference. Users can specify the model used for genotype likelihood estimation, as well as the output format. For more information on the available options, see the ANGSD.