DepthKopy v1.5.0

For a better rendering and navigation of this document, please download and open ./docs/depthkopy.docs.html, or visit https://slimsuite.github.io/depthkopy/. Documentation can also be generated by running DepthKopy with the dochtml=T option. (R and pandoc must be installed - see below.)

DepthKopy is an updated version of the regcnv methods of Diploidocus, using the same single-copy read depth estimation from BUSCO Complete genes as DepthSizer.

DepthKopy needs a genome assembly (fasta format, seqin=FILE), a set of long read (ONT, PacBio or HiFi) data for the assembly (reads=FILELIST and readtype=LIST) or mapped BAM file (bam=FILE), and a BUSCO/BUSCOMP full table of results (busco=TSVFILE) or pre-calculated single-copy read depth (scdepth=NUM). For kmer analysis, reads also need to be provided with kmerreads=FILELIST (and 10xtrim=T if barcoded 10x linked reads).

Optionally, it can then take one or more additional sources of assembly regions for which stats will be calculated. Delimited text files or GFF files can be provided using regfile=FILE for a single file, or regfile=CDICT for multiple files. Multiple files are given in the form Name1:File1,Name2:File2,...,NameN:FileN, from which each unique Name will be plotted as a separate violin (see output). GFF files will be parsed to extract any features matching the types given by gfftype=LIST (default = gene). Delimited files will need to have headers that match those provided by checkfields=LIST (default = SeqName,Start,End). By default, 100 kb non-overlapping windows will also be output (winsize=INT winstep=NUM). Unless seqstats=F, statistics will also be calculated per assembly scaffold. Subsets of scaffolds can be given separate window output plots using chromcheck=LIST (scaffold names) or chromcheck=INT (min size).

DepthKopy works on the principle that Complete BUSCO genes should represent predominantly single copy (diploid read depth) regions along with some poor quality and/or repeat regions. Assembly artefacts and collapsed repeats etc. are predicted to deviate from diploid read depth in an inconsistent manner. Therefore, even if less than half the region is actually diploid coverage, the modal read depth is expected to represent the actual single copy read depth.

DepthKopy uses samtools mpileup (or samtools depth if quickdepth=T) to calculate the per-base read depth. This is converted into an estimated single copy read depth using a smoothed density plot of BUSCO single copy genes. BUSCO single-copy genes are parsed from a BUSCO full results table, given by busco=TSVFILE (default full_table_$BASEFILE.busco.tsv). This can be replaced with any table matching the BUSCO fields: ['BuscoID','Status','Contig','Start','End','Score','Length']. Entries are reduced to those with Status = Complete and the Contig, Start and End fields are used to define the regions that should be predominantly single copy. Output from BUSCOMP is also compatible with DepthKopy. DepthKopy has been tested with outputs from BUSCO v3 and v5.

Output is a table of depth statistics and predicted copy number for each input dataset (BUSCO Complete, BUSCO Duplicated, regfile Regions, assembly scaffolds (seqstat=T) and sliding windows). Data visualisations are also output for each region set using ggstatsplot.

DepthKopy has been published as part of the Waratah genome paper:

Chen SH, Rossetto M, van der Merwe M, Lu-Irving P, Yap JS, Sauquet H, Bourke G, Amos TG, Bragg JG & Edwards RJ (2022). Chromosome-level de novo genome assembly of Telopea speciosissima (New South Wales waratah) using long-reads, linked-reads and Hi-C. Molecular Ecology Resources doi: 10.1111/1755-0998.13574

Please contact the author if you have trouble getting the full text version, or read the bioRxiv preprint version:

Chromosome-level de novo genome assembly of Telopea speciosissima (New South Wales waratah) using long-reads, linked-reads and Hi-C. bioRxiv 2021.06.02.444084; doi: 10.1101/2021.06.02.444084.

DepthKopy is written in Python 2.x and can be run directly from the commandline:

python $CODEPATH/depthkopy.py [OPTIONS]

If running as part of SLiMSuite, $CODEPATH will be the SLiMSuite tools/ directory. If running from the standalone DepthKopy git repo, $CODEPATH will be the path the to code/ directory.

Unless bam=FILE is given, minimap2 must be installed and either added to the environment $PATH or given to DepthSizer with the minimap2=PROG setting, and samtools needs to be installed. R and tidyverse will also need be installed, ideally with the ggstatsplot and writexl libraries also installed. To generate documentation with dochtml, a pandoc environment variable must be set, e.g.

export RSTUDIO_PANDOC=/Applications/RStudio.app/Contents/MacOS/pandoc

For DepthKopy documentation, run with dochtml=T and read the *.docs.html file generated.

A list of commandline options can be generated at run-time using the -h or help flags. Please see the general SLiMSuite documentation for details of how to use commandline options, including setting default values with INI files.

### ~ Main DepthKopy run options ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ###
seqin=FILE      : Input sequence assembly. (Must be *.fa or *.fasta for kmerself=T.) [None]
basefile=FILE   : Root of output file names [diploidocus or $SEQIN basefile]
scdepth=NUM     : Single copy ("diploid") read depth. If zero, will use SC BUSCO mode [0]
bam=FILE        : BAM file of long reads mapped onto assembly [$BASEFILE.bam]
bamcsi=T/F      : Use CSI indexing for BAM files, not BAI (needed for v long scaffolds) [False]
reads=FILELIST  : List of fasta/fastq files containing reads. Wildcard allowed. Can be gzipped. []
readtype=LIST   : List of ont/pb/hifi file types matching reads for minimap2 mapping [ont]
dochtml=T/F     : Generate HTML DepthKopy documentation (*.docs.html) instead of main run [False]
### ~ Depth and Copy Number options ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ###
busco=TSVFILE   : BUSCO full table [full_table_$BASEFILE.busco.tsv]
quickdepth=T/F  : Whether to use samtools depth in place of mpileup (quicker but underestimates?) [False]
depchunk=INT    : Chunk input into minimum of INT bp chunks for temp depth calculation [1e6]
deponly=T/F     : Cease execution following checking/creating BAM and fastdep/fastmp files [False]
depfile=FILE    : Precomputed depth file (*.fastdep or *.fastmp) to use [None]
homfile=FILE    : Precomputed homology depth file (*.fasthom) to use [None]
regfile=CDICT   : List of Name:Files (or single FILE) of SeqName, Start, End positions (or GFF) for CN checking [None]
checkfields=LIST: Fields in checkpos file to give Locus, Start and End for checking [SeqName,Start,End]
gfftype=LIST    : Optional feature types to use if performing regcheck on GFF file (e.g. gene) ['gene']
winsize=INT     : Generate additional window-based depth and CN profiles of INT bp (0 to switch off) [100000]
winstep=NUM     : Generate window every NUM bp (or fraction of winsize=INT if <=1) [1]
chromcheck=LIST : Output separate window analysis violin plots for listed sequences (or min size) + 'Other' []
seqstats=T/F    : Whether to output CN and depth data for full sequences as well as BUSCO genes [True]
### ~ Rscript options ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ###
outdir=PATH     : Redirect the outputs of the depthcopy.R script into outdir [./]
pointsize=INT   : Rescale the font size for the DepthKopy plots [24]
### ~ KAT kmer options ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ###
kmerself=T/F        : Whether to perform additional assembly kmer analysis [True]
kmerreads=FILELIST  : File of high quality reads for KAT kmer analysis []
10xtrim=T/F         : Whether to trim 16bp 10x barcodes from Read 1 of Kmer Reads data for KAT analysis [False]
### ~ System options ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ###
forks=X         : Number of parallel sequences to process at once [0]
killforks=X     : Number of seconds of no activity before killing all remaining forks. [36000]
forksleep=X     : Sleep time (seconds) between cycles of forking out more process [0]
tmpdir=PATH     : Path for temporary output files during forking [./tmpdir/]
### ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ###

The main inputs for DepthKopy copy number prediction are:

• seqin=FILE : Input sequence assembly [Required].
• reads=FILELIST: List of fasta/fastq files containing long reads. Wildcard allowed. Can be gzipped. A BAM file can be supplied instead with bam=FILE.
• readtype=LIST : List of ont/pb/hifi file types matching reads for minimap2 mapping [ont]
• busco=TSVFILE : BUSCO(MP) full table [full_table_$BASEFILE.busco.tsv] used for calculating single copy ("diploid") read depth.

The first step is to generate a BAM file by mapping reads on to seqin using minimap2. A pre-generated BAM file can be given instead using bam=FILE. There should be no secondary mapping of reads, as these will inflate read depths, so filter these out if they were allowed during mapping. Similarly, the BAM file should not contain unmapped reads. (These should be filtered during processing if present.) If no BAM file setting is given, the BAM file will be named $BASEFILE.bam, where $BASEFILE is set by basefile=X.

DepthKopy works on the principle that Complete BUSCO genes should represent predominantly single copy (diploid read depth) regions along with some poor-quality and/or repeat regions. Assembly artefacts and collapsed repeats etc. are predicted to deviate from diploid read depth in an inconsistent manner. Therefore, even if less than half the region is actually diploid coverage, the modal read depth is expected to represent the actual single copy read depth. This is estimated using a smoothed density distribution calculated using R density(). BUSCO single-copy genes are parsed from a BUSCO full results table, given by busco=TSVFILE (default full_table_$BASEFILE.busco.tsv). This can be replaced with any table with the fields: ['BuscoID','Status','Contig','Start','End','Score','Length']. Entries are reduced to those with Status = Complete and the Contig, Start and End fields are used to define the regions that should be predominantly single copy. BUSCOMP v0.10.0 and above will generate a *.complete.tsv file that can be used in place of BUSCO results. This can enable rapid re-annotation of BUSCO genes following, for example, vector trimming with Diploidocus.

DepthKopy uses samtools mpileup (or samtools depth if quickdepth=T) to calculate the per-base read depth and extracts the smoothed modal read depth for all single-copy (Complete BUSCO genes) using the density() function of R. To avoid a minority of extremely deep-coverage bases disrupting the density profile, the depth range is first limited to the range from zero to 1000, or four time the pure modal read depth if over 1000. If the pure mode is zero coverage, zero is returned. The number of bins for the density function is set to be greater than 5 times the max depth for the calculation.

By default, the density bandwidth smoothing parameter is set to adjust=12. This can be modified with depadjust=INT. The raw and smoothed profiles are output to *.plots/*.raw.png *.plots/*.scdepth.png to check smoothing if required. Additional checking plots are also output (see Outputs below).

The full output of depths per position is output to $BAM.fastmp (or $BAM.fastdep if quickdepth=T). The single-copy is also output to $BAM.fastmp.scdepth. By default, generation of the fastdep/fastmp data is performed by chunking up the assembly and creating temporary files in parallel (tmpdir=PATH). Sequences are batched in order such that each batch meets the minimum size criterion set by depchunk=INT (default 1Mbp). If depchunk=0 then each sequence will be processed individually. This is not recommended for large, highly fragmented genomes. Unless dev=T or debug=T, the temporary files will be deleted once the final file is made. If DepthSizer crashed during the generation of the file, it should be possible to re-run and it will re-use existing temporary files.

For each region analysed, the same density profile calculation is used to predict the dominant read depth across the region, which is then converted into copy number (CN) by dividing by the single-copy read depth. Confidence intervals are calculated based on random sampling of the observed single copy read depth. (Details to follow: available on request.) By default, the R script will parallelise this using the number of threads set with forks=INT. If this causes memory issues, it can be forced to run with a single thread using memsaver=T.

Finally, CN and depth statistics are integrated with kmer and within-assembly homology data, and output as a series of plots and tables (below).

The main DepthKopy outputs are:

• *.regcnv.tsv = CN and depth statistics for different input region datasets.
• *.xlsx = compiled Excel file of each region file.
• *.log = DepthKopy log file with key steps and details of any errors or warnings generated.
• *.plots/ = Directory of PNG plots (see below)

The BUSCO input file will get *.regcnv.tsv and *.dupcnv.tsv copy number prediction tables generated. Other region datasets will have a *.regcnv.tsv file. These will have the following fields added:

• MeanX = Mean sequencing depth across region.
• MedX = Median sequencing depth across region.
• ModeX = Pure modal sequencing depth across region.
• DensX = Density-plot adjusted modal sequencing depth across region.
• SelfK = Density-plot adjusted modal kmer frequency across region. (In development.)
• HomPC = Percentage of region with within-assembly homology identified. (In development.)

DepthKopy will also generate a number of plots of results, which are output in $BASE.plots/ as PDF and PNG files.

First, the raw and smoothed read depth profiles will be output to:

• $BASE.plots/$BASE.raw.png = raw depth profile
• $BASE.plots/$BASE.scdepth.png = smoothed depth profile with SC depth marked

In addition, violin plots will be generated for the following $STAT values:

• $BASE.plots/$BASE.CN.png = Estimated copy number.
• $BASE.plots/$BASE.MeanX.png = Mean depth of coverage
• $BASE.plots/$BASE.MedX.png = Median depth of coverage
• $BASE.plots/$BASE.ModeX.png = Pure Modal depth of coverage
• $BASE.plots/$BASE.DensX.png = Smoothed density modal depth of coverage

If seqstats=T then each assembly sequence will also be output in a Sequences violin plot for comparison. Each point in the plots is a separate gene or sequence.

Values for individual genes/sequences are also output as a density scatter plot named $BASE.plots/$BASE.$REGIONS.$STAT.png, where

• $REGIONS is the type of region plotted (BUSCO complete, Duplicated, or assembly Sequences).
• $STAT is the output statistic: MeanX, MedX, ModeX or DensX.

If DepthKopy fails to run to completion, or you wish to change axes limits and/or edit out some of the features, you can also generate plots from the *.xlsx output file using the depthcopyplot.R file:

Rscript $CODEPATH/depthcopyplot.R basefile=$BASE [scdepth=NUM] [xlsx=FILE] [xsheets=LIST] [reghead=LIST] [pngdir=PATH] [cnmax=INT] [sigdif=T/F] [rdir=PATH]