Scripts for grysbok (Raphicerus) mitogenomes announcement

Contains scripts used in the manuscript: The complete mitogenomes of the two grysbok species: Raphicerus melanotis Thunberg 1811 (Cape grysbok) and Raphicerus sharpei Thomas 1897 (Sharpe’s grysbok).

• The mitogenomes were assembled from the reads generated for high-coverage nuclear genome sequencing, hence the fastp read filtering step, and the subsampling step to 5Gbp of reads to reduce computational time and avoid excessive coverage.
• Numbers at the start of the file names correspond to the steps below and are not technically part of the file name (so they must be removed for scripts calling other scripts to work, or the names must be updated inside the scripts).
• Sample names (RmM001 and RsM001) are often replaced by the asterisk (*) wildcard/placeholder below.

1. Process raw reads

• The fastp_*.slurm scripts call the fastp_wrapper_modern.sh script, which contains the actual fastp command and parameters (from line 107). Type fastp_wrapper_modern.sh -h for help.

2. Subsample clean reads

• Contains the reformat.sh commands to subsample the clean reads to ~5Gbp, which were then used to assemble the mitogenomes. The reformat.sh script/command is part of the BBMap set of tools.

3. Assemble mitogenomes with NOVOPlasty

• The *_novoplasty.slurm scripts call NOVOPlasty and provide the path to the config_*_5Gbp.txt files, which contain the actual settings for the NOVOPlasty assembly.

4. Construct phylogenetic tree with IQ-TREE

• The data for maximum-likelihood phylogentic analysis in IQ-TREE was prepared in Geneious Prime to extract, concatenate, and align the coding regions of the mitogenomes (control region excluded) into three files:
  • n7_PCGs_MAFFT.fasta: Contains the 13 protein-coding genes (PCGs) MAFFT alignment with stop codons removed.
  • n7_rRNAs_MAFFT.fasta: Contains the 2 rRNA genes MAFFT alignment.
  • n7_tRNAs_MAFFT.fasta: Contains the 22 tRNA genes MAFFT alignment.
• The n7_charset.nex file is given to IQ-TREE. It defines the partitions and tells IQ-TREE in which files to find the alignments, which positions in the alignments are relevant for each partition (an asterisk (*) means use all positions in the alignment), and which type of substitution models to use. For example, for the PCGs partition, the flag CODON2 is added, because it is an in-frame alignment of codons, so IQ-TREE uses codon models with the vertebrate mitochondrial code (2) instead of standard DNA substitution models. For the RNA alignments, no flag is added, so IQ-TREE will detect the data type (DNA) and use standard DNA substitution models. The flag MFP+MERGE results in IQ-TREE testing all relevant models for each partition, choosing the best-fitting model (chosen with BIC), and then checking whether any partitions have the same best-fitting model and if so merges them into a single partition. This was the case for the rRNA and tRNA partitions. IQ-TREE then automatically proceeds with the maximum-likelihood tree estimation using the best partition scheme and models for each partition.
• The ultrafast bootstrap and SH-aLRT methods were used in IQ-TREE, which have slightly different interpretations than traditional bootstrap values. However, IQ-TREE still outputs a consensus tree with traditional bootstrap support values in the output files with the .iqtree and .contree extensions (not included here), and this is what was reported in the manuscript.
• The iqtree_Raphicerus_mitogenomes.slurm contains the commands to run the analysis.
• Note: The evolutionary divergence/p-distance analysis in the manuscript (done in MEGA11) used these same alignments, except that they were manually concatenated into a single file for use in MEGA.

5. Make mitogenome coverage plots (Figure S1)

The file GetCoverage_ForCircularMitogenome.txt explains how the plots in Figure S1 were generated. In short, the reads that were assembled by NOVOPlasty into the final mitogenome were saved (see the config files in Step 3 above) and these were mapped back to the relevant circular mitogenome reference to get the per-base read depth/coverage. Simple as it sounds, it is a bit more complicated than that:

• If you use BWA to map reads to the linearised fasta file of the assembled mitogenome (as I initially did), then you get an artificial decrease in coverage towards the ends of the scaffold/mitogenome, because reads that map to the ends and have overhangs are not mapped (depending on the size of the overhang, I assume), so fewer reads are mapped than expected, resulting in an artifical decrease in coverage at the ends.
• Thus, the reads were mapped using the Geneious mapper, because it is aware of the circular nature of the reference and can map reads across the start/end of the mitogenome.
• The *_NovoplastyAssembledReads_GeneiousMapped_Coverage.csv files contain the per-base read depth/coverage and were exported from Geneious after mapping. The Plot_coverage_mito_Revision1.R then uses these csv files to make the plot in R.
• Another way to avoid the artificial decrease in coverage at the ends of the mitogenome is to concatenate short pieces (e.g. 40 bases) from the beginning and end of the mitogenome to the opposite end, so that reads can map across the artificial breakpoint, as in Morin et al. (2015) (open access), but I did not try this approach here. I think the length of the concatenated pieces are important, because if they are too long, then reads will map equally well to multiple locations in the reference and might be discarded, or a random one is kept (depends on the aligner settings), which then doesn't really solve the problem. So don't make these pieces too long, and consider your read length.