Breaking Barriers in Candida auris Genomics: Analysis Tools for Whole Genome Sequencing Amid Database Scarcity

Quality Check

The first inevitable step after NGS is to perform QC check (Table 1) and to analyze the quality of obtained NGS data. Standard quality control of reads obtained after short-read sequencing (typically from Illumina to Ion torrent) is performed using FastQC or Trimmomatic. FastQC provides a visual illustration of distribution of GC content, nucleotide bias, and base quality. Graphs and tables generated by FastQC helps quickly to assess data quality. If poor quality reads are obtained then Trimmomatic software, which is available freely in GALAXY, can be used to remove or trim those reads to obtained better quality data.⁴² GALAXY is an open-source, web-based platform for computational analysis. It allows scientists to analyze, visualize, share, and efficiently manage their data. After removing spurious and contaminated reads, data should be reanalyzed using FastQC to check the final quality of reads.

Similarly, PRINSEQ tool can be used to filter, reformat, or trim genomic data obtained from long-read sequencing technologies (like Pacific Biosciences and Oxford Nanopore). PRINSEQ tool is best used for metagenomic data but has also been used for C. auris to assess the quality of NGS data and to filter out low-quality sequences.^40,42 Another tool for quality control is QUALIMAP for long reads and can be used as graphical user interface (GUI) or command line interface.⁴¹ It analyzes alignment sequence data file for the biases and helps in taking decision for further analyzes. The advantage with QUALIMAP is that it can be used to assess the data quality irrespective of the NGS platform being used. Another option of tools to trim reads are FASTX-Toolkit, long-read error correction (LoRDEC), and Porechop.³⁷ FASTX-Toolkit is a collection of command line tools, which can be used to filter data based on different quality parameters and trims fastq or fasta reads.³⁸ LoRDEC is mostly used for PacBio long reads since error rate is high and it uses graph based approach to fetch corrective sequence from erroneous sequence data. On similar lines, Porechop is used for data generated on Oxford Nanopore platform to find and trim adapter reads from the sequences.⁴¹ Porechop also provides demultiplexing support for reads generated by Oxford Nanopore.

Genome Assembly

Using better quality of reads, genome is assembled so as to reconstruct and recreate original chromosomal structure of organism. Tools like SPAdes, Velvet, Flye, PHRAP, NECAT, Canu, SMRT analysis, etc. are used for genome assembly (Table 2).^43-45 For de novo sequencing, SPAdes assembler is used, which an open access downloadable tool and is also available in Illumina’s basespace sequence hub platform (https://basespace.illumina.com/).⁴⁵ Velvet assembler uses algorithm of de Bruijn graph to assemble short reads and resolve the small repeats and produces contigs of reasonable length.⁴⁶ The mean number of contigs is lower and the mean N50 value is higher in the genomes assembled using SPAdes, which is more desirable in comparison to Velvet assembled genomes. The shortest contig’s sequence length at 50% of the assembly’s overall length is known as N50, which determines the assembly’s quality. Flye assembler is mostly used for long sequencing reads such as those produced by PacBio and Oxford Nanopore Technologies. For assembly of shotgun sequence data, PHRAP assembler is used. PHRAP assembly program follows match “words” criteria and can also change parameters for assembly using command line. For assembling large genome which uses long noisy reads generated by Nanopore, NECAT assembler is used. Canu provides graph-based assembly outputs of higher N50 in GFA format. Canu assembler can be used for both, short as well as long reads generated by PacBio and Oxford Nanopore Technologies. SmrtAnalysis suite v2.3 was used to assemble the single-molecule real-time (SMRT) sequencing reads of C. auris that were generated on PacBio RS II SMRT DNA sequencing system (Pacific Biosciences).⁴⁷

Since there are many assemblers available for use, choice of assembler will depend on several factors including characteristics of the sequencing data, such as the size and complexity of the genome and specific research question or application. It is important to select an assembler that is compatible with sequencing data generated for your experiment. When evaluating assemblers, it is important to consider performance metrics such as contiguity, completeness, and accuracy of the assembled genome.

Assembled genome needs to be polished for better representation of genome from the read data. Pilon and Medaka are the tools used for polishing of assembly.⁴⁹ Pilon is a free online automated tool and by correcting bases, fixing mis-assemblies, and filling gaps, it improves draft genome assemblies. Medaka uses signal-based and sequence-graph based logic to create consensus sequences and generate polished genomic data.

The genome of C. auris is relatively large (around 12.7 Mb) and complex, with multiple copies of some genes and regions of repetitive DNA.⁵⁰ Candida species exhibit significant genetic diversity, both within and between species.⁵¹ This can make it challenging to compare genomes and identify conserved and variable regions. Accurate gene predictions and functional annotations is a concern because C. auris is poorly characterized.⁵²

Sequence Alignment

To identify regions of similarity or nucleotide variations that may be a consequence of evolutionary relationships, alignment of sequences is done. Sequence alignment (Table 3) can be local or global. Three algorithms make up the BWA software package: BWA-backtrack, BWA-SW, and BWA-MEM. BWA-MEM is faster, most accurate, and highly recommended for long sequences.⁵⁴ BWA focuses on local alignment and hence need indexes for alignment. To perform multiple sequence alignment, ClustalW software is used, which is available online for free, but comes with limitation of file upload size. ClustalW uses progressive algorithm to perform global alignment of sequences. MUSCLE uses iterative algorithm and run combination of both global and local alignment. TOPHAT2 is alignment tool for transcriptome data of RNA seq. Bowtie 2 supports gapped, local, and paired-end alignment modes with most accurate and fast result.⁵⁵ MUMmer/nucmer tool is basically used for alignment of whole genome of organism and generates results in minimum amount of time.⁵⁶ MUMmer/nucmer is a tool to check “Maximal Unique Matches” and is best used for sequences with high similarity but have large rearrangements.

After alignment of sequences, either BAM, SAM, or CRAM file is generated. With the use of samtools, data can be sorted, merged, and converted between the formats. It can also create indexes of the alignment file as well as can mark duplicates.⁵⁷ Sequence alignment/map (Samtools) can be very useful in calculating coverage of sequences and read depth at each position or region. Similar to samtools, PICARD tool can be used for sorting, format conversion, to locate, and tag duplicate reads in a BAM or SAM file.

New SAM/BAM file obtained after edit can be visualized using Integrative Genomics Viewer (IGV).⁵⁸ In IGV, user can upload BAM file along with their reference genomes to have clear visualization of mutations at relevant positions of chromosomes. IGV can also be used to view alignment, depth of reads and regions of poor reads. IGV can either be downloaded or can be used on web server independently. Along with BAM file, index file is necessary to be kept under same folder in local machine for error free visualization in IGV. Previously published study, used IGV to visualize rearrangement sites based on alignment of reads of C. auris.⁴⁰

Gene Prediction

Process of identifying the coding region in genome is called gene prediction (Table 4). The output of gene prediction is the position of start, stop, and open reading frame of gene in chromosome or scaffold. Gene prediction can be done for both prokaryotic as well as eukaryotic organism. AUGUSTUS, based on a Generalized Hidden Markov Model, predicts genes for eukaryotic as well prokaryotic genome on web server.⁵⁹ AUGUSTUS not only predicts coding parts of the genes but also gives respective protein sequences. Glimmer3 software is used for gene prediction for prokaryotes and is available for free on the GALAXY tool package.

For all those predicted genes, functional annotation is required, which can be done using software like MAKER2, BRAKER1, BRAKER2, and GENEMARK.^60,61 For adding annotation to genes, trained dataset is required. Annotation quality can be improved by using genomic GFF file available in databases and mRNA-sequence data. The functional annotation of C. auris and other Candida species was performed using the Clusters of Orthologous Groups of proteins (COGs) database.⁶²

For C. auris, gene level data is not appropriately documented in the NCBI or “Candida Genome Database”.⁴⁴ For example, “FKS1” is annotated as “CJI97_001351” and “ERG11” as “CJI97_00156.” Similarly, description for annotation, “CJI97_000015,” is given as “hypothetical protein” and role of this gene is unclear. TAC1b, which is one of the important gene responsible for antifungal resistance, cannot be searched in the databases by its name. This major limitation hinders postsequencing data analysis process and needs to be addressed by the scientific community.

Orthologs Identification

Orthologs genes are common genes evolved due to speciation event and can be traced to ancestral genome if same gene exists. Using orthology relationships, the biologists can predict functional data like genome annotation, comparative genomics, etc., computationally. In case of C. auris, closest organism is Candida albicans for which data are available on all platforms or databases. Also, biofilm and antimicrobial mechanism are studied for C. albicans in depth, hence information available on C. albicans can give functional insights for C. auris genome.⁶³ OrthoMCL and BUSCO are two software used for identification of orthologs.⁶⁴ OrthoMCL approach is similar to INPARANOID algorithm with some exceptions of parameters used during analyses.⁶⁵ BUSCO is an open source software available in conda package for assessment of orthologs. In the study focusing on clade specific rearrangement of C. auris, orthologs were assigned using bidirectional blast and Markov index 1.5; maximum e-value 1e-5 were parameters used in OrthoMCL.⁴⁰ Chatterjee et al. used orthologues of C. albicans and reported that C. auris shared significant virulence attributes, including secreted proteases, genes involved in biofilm formation, mannosyl transferases, and oligopeptide transporters.⁵⁰

Variant Calling

Variant calling uses high-throughput sequencing technology to calculate genomic variations like single nucleotide polymorphisms (SNPs), insertions, deletions, structural variations, and copy number variations (Table 5). GATK (Genome Analysis Toolkit) Haplotype Caller is one of the promising tool used for variant calling.⁶⁶ Output of the Haplotype caller is variant calling file (VCF) and VCF index files. For variant calling, reference of Candida genome is necessary. Reference genomes can either be downloaded from NCBI or de novo assembled genome can be used. Freebayes is based on Bayesian algorithm, which checks whole haplotype of the read independently of alignment positions.⁶⁷ VarScan2 employs a heuristic approach to call variance and it relies on threshold set for various parameters to determine variants.⁶⁸ VarScan2 performs variant calling based on aligned read and by going through each read position by position.

SnpEff and SnpSift are tools to annotate genomic variants and predict its functional effects.^69,70 SnpEff can also provide noncoding and regulatory annotations. SnpSift helps in filtering and manipulation of VCF files like joining, extracting fields, interval indexing, and calculation of concordance of two VCF files etc.⁷¹ NASP is a comprehensive, open-source tool, which is easily executable and aids in SNP identification as well as offers differentiation between large numbers of microbial genome datasets.³² Another software for SNP identification is CLC Genomics workbench (Qiagen), which is a commercially available comprehensive analysis package that supports analysis and visualization of NGS data.⁵³ VCFtools is another open source tool that is designed to work on VCF/BCF format to filter, summarize, manipulate, and convert data into other useful file formats.⁴⁷ It is a very easy tool to filter out variants with less quality score and uncharacterized genes of C. auris.⁴⁷ Variants are filtered out based on several quality parameters such as QualByDepth (QD) or FisherStrand (FS) or mapping quality (MQ). QD is assessment of the variant’s quality relative to the sequencing depth. FS is Phred score to evaluate the variation being seen on only the forward or only the reverse strand. Root mean square of the MQ of the reads across all samples in a run is known as RMSMappingQuality. To further enhance data accuracy and to remove technical biases, hard filters such as MQRankSum, ReadPosRankSum, and SOR can be applied. Variants can also be filtered out based upon genotype quality and alternate percent allele.

Prediction of tRNAs and rRNAs

tRNA and rRNA are associated with each other and are involved in protein synthesis. tRNAscan is widely used tool to predict tRNA genes in genomic sequences by using covariance model.⁷² Structural information of tRNA is stored as training dataset, which is then used to search for complete gene details in query sequences. Output of tRNAscan tool is genomic coordinate, predicted function, and secondary structure of tRNA. Hidden Markov Model is used as foundation for RNAmmer software.⁷³ RNAmmer is a free webserver tool that uses ERRD to predict ribosomal genes by mapping aligned data. Very few researchers have actually worked on prediction of tRNAs and rRNAs in Candida and further utility of this information is unclear.⁵³

Phylogenetic Analysis

C. auris has five clades (I–V) based upon the number of SNPs identified after alignment. These clades are highly identical, sharing an average pairwise nucleotide of 98.7%. Clade II and Clade III appear most similar and share 99.3% identity.⁴⁷ Whole genome based large SNP data file is used to study relatedness among the isolates. Based on SNP density, average divergence for each chromosome and whole genome is calculated. Using aligned file, phylogeny of multiple samples can be generated together to depict evolutionary relationship among each other, hence very useful in comparative genomics. Maximum Parsimony, Neighbour Joining, UPGMA, Maximum Likelihood, and MCMC are some of the algorithm used for generating phylogeny by MEGA software (Table 6).⁷⁴ Each tree can be represented as cladogram or dendrogram. MEGA is a very popular tool used for computational molecular evolution analysis. It is an interactive software and has integrated feature of alignment, tree construction, and tree visualization. MEGA is a free tool that can be downloaded and installed on any computer system. RAxML is another popular software, which return trees with good likelihood score.⁷⁵ RAxML uses 1,000 bootstrap replicates with the GTR model of nucleotide substitution and γ-distributed rates to speed up the process and gives very fast result of large data. Based upon Bayesian phylogenetic inference methods, BEAST software uses MCMC algorithm to construct phylogeny.⁷⁶ IQ-Tree and FastTree are different tools used for tree visualization.⁷⁷ To better understand the spread and evolution of C. auris over time, it is important to perform phylogenetic analysis on strains from different geographic locations. Phylogeny also helps in study of transmission of C. auris within a particular healthcare facility.

Genomic Sequence Databases

There are limited databases, which hold information of C. auris and its genetic clades. One of the databases is CANDIDAGENOME (http://www.candidagenome.org/), which has information on genomic sequences, genes, and protein of C. albicans, Candida dubliniensis, Candida parapsilosis, Candida glabrata, and C. auris.^78-80 Candidagenome database has lot of information for C. albicans, but minimal data is available for C. auris.⁴⁴ Hence, for the researchers who are working on different clades of C. auris, Candidagenome database is of limited use.

NCBI is another database, which is a premium source of knowledge on Candida species. However, in case of C. auris, most of the genomic data are at scaffold level in NCBI (e.g., C. auris strain B11221) and very few at chromosome level (e.g., C. auris strain B11220). Since chromosome/scaffold level information is not consistent or uniform among clades, it hampers downstream analysis.

FungiDB provides a range of analysis tools and resources to facilitate data exploration and analysis.⁸¹ These tools include genome browsers, BLAST search functionality, sequence alignment tools, and data visualization tools. It provides gene sequence along with mRNA and protein sequences, which can be utilized for tertiary analysis such as identification of drug resistance mutations after sequencing. FungiDB has benefit over other databases in terms of gene information, which is not be present in NCBI explicitly.

Another database is MARDy, Mycology Antifungal Resistance Database version 1.1, which has limited information (mutations associated with drug resistance) for three genes, namely ERG11, FKS1, and FUR1 for C. auris and C. albicans.⁸² The information available in the database can be retrieved by applying relevant filters such as name of the organism, gene, or drug.

JGI (Joint Genome Institute) MycoCosm is a web-based platform that provides access to genomic and functional data for fungal organisms.⁸³ The platform allows users to search and browse the available fungal genomes, explore gene annotations, and access various tools for analyzing and visualizing the data. Researchers can investigate specific genes, metabolic pathways, or genetic features of interest, and compare them across different fungal species.⁸⁴ MycoCosm is linked to another database such as Ensembl fungi, which has sequence level information and focuses specifically on providing genomic data, annotations, and analysis tools for fungal species.⁸⁵ Both MycoCosm and Ensembl fungi do not have information for C. auris but does have information on other species of Candida.

Antimicrobial Resistance

Most of the researchers have determined Minimal Inhibitory Concentration (MIC₅₀) for antifungal susceptibility testing for drugs like fluconazole, caspofungin, and amphotericin B for C. auris.^32,38,47,53 At molecular level, the resistant mutations are identified by Sanger sequencing, allele-specific PCR, or real-time PCR.^48,84 Some of the researchers have managed to figure out antifungal resistant mutations by manually aligning the whole genome sequence to the reference genome by ClustalW alignment using MEGA.⁷⁴ In case of fungi, the utility of the WGS data is limited as far as antifungal resistance is concerned because there are no automated tools to retrieve mutations associated with AMR. Instead, different gene sequences of a particular isolate have to be manually aligned to the reference gene sequences to pick up nucleotide variations.