{'fraction_non_alphanumeric': 0.03733998798465733, 'fraction_numerical': 0.02345302463145247, 'mean_word_length': 5.322717311906501, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 7, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
GIP: an open-source computational pipeline for mapping genomic instability from protists to cancer cells

Abstract Genome instability has been recognized as a key driver for microbial and cancer adaptation and thus plays a central role in many diseases. Genome instability encompasses different types of genomic alterations, yet most available genome analysis software are limited to just one type of mutation. To overcome this limitation and better understand the role of genetic changes in enhancing pathogenicity we established GIP, a novel, powerful bioinformatic pipeline for comparative genome analysis. Here, we show its application to whole genome sequencing datasets of Leishmania, Plasmodium, Candida and cancer. Applying GIP on available data sets validated our pipeline and demonstrated the power of our tool to drive biological discovery. Applied to Plasmodium vivax genomes, our pipeline uncovered the convergent amplification of erythrocyte binding proteins and identified a nullisomic strain. Re-analyzing genomes of drug adapted Candida albicans strains revealed correlated copy number variations of functionally related genes, strongly supporting a mechanism of epistatic adaptation through interacting gene-dosage changes. Our results illustrate how GIP can be used for the identification of aneuploidy, gene copy number variations, changes in nucleic acid sequences, and chromosomal rearrangements. Altogether, GIP can shed light on the genetic bases of cell adaptation and drive disease biomarker discovery.


INTRODUCTION
In recent years, the field of genomics has rapidly expanded with a fast increase in the number of newly sequenced genomes (1). This surge is a direct consequence of the development of new and ever more efficient, high-throughput capable sequencing technologies (2). On the one hand, the improvement in long reads technology allowed the generation of high-quality genome assemblies (3,4). On the other hand, the decreasing costs for short-reads sequencing and the parallel increase in sequencing throughput propelled the exponential increase of available whole genome sequencing (WGS) data (5). Thanks to these advances one can reasonably expect WGS to rapidly become a key component of personalized medicine and clinical applications.
Indeed, WGS technologies have revolutionized many applications in the fields of medicine and microbiology. WGS can be used for clinical sample screenings in the diagnosis, classification and surveillance of microbial pathogens (6). WGS allows strain identification with better resolution compared to genetic marker-based methods and can inform on the accessory genome of microbes undergoing horizontal gene transfer (6,7). Other important applications of WGS include antimicrobial resistance (AMR) profiling (8), the identification of candidate antigens in vaccine development (9,10), the tracking of outbreaks within hospitals and communities (6,(11)(12)(13)(14)(15) and microbial evolutionary adaptation (16)(17)(18), which can have profound effects on human infections. In this adaptation process, cycles of genetic mutation and environmental selection lead to microbial fitness gain, resulting in drug resistance and shifts in tropism or virulence. As such, genome instability often determines disease outcome (19)(20)(21)(22) and is a key driver of phenotypic and genetic variability of microbes and other pathogenetic systems that rely on genome instability for adaptation, such as cancer.
In this context, several consortia-based projects have been established with the goal to produce WGS for the study of different biological systems (5), and a number of publicly available databases have been compiled or updated (23). Parallel to data availability, many bioinformatics tools have been developed to perform specific genome analysis tasks (24,25). For instance, tools such as Freebayes (26), CN-Vnator (27) and DELLY (28) have been respectively used for the detection or characterization of DNA single nucleotide variants (SNVs), copy number variations (CNVs),

Genome sequencing and genome assembly data
WGS reads were downloaded from the Sequence Read Archive (SRA) (33), the European Nucleotide Archive (ENA) (34) repositories and the Encyclopedia of DNA Elements (ENCODE) dashboard (35) (Supplementary Table  S1). For L. infantum the GCA 900500625 genome reference and gene annotations available from the ENSEMBL protists server (release-48) were used (36). For C. albicans the assembly 21 of the SC5314 strain genome reference and gene annotations available from the Candida Genome Database (CGD) were used (37). For P. vivax the P01 reference genome and gene annotations available from Plas-moDB (release-50) were used (38). For the cancer cell lines, the human genome GRCh38 primary assembly and gene annotations available from ENSEMBL (release-102) were used (36).

GIP and giptools
All results presented in this study were generated using GIP and giptools version 1.0.9. GIP code is maintained and freely distributed at the github page: https://github.com/ giovannibussotti/GIP. Figure 1 and Supplementary Figure  S1 provide a schematic representation of the GIP workflow. The giptools container is accessible from the Singularity cloud at https://cloud.sylabs.io/library/giovannibussotti/ default/giptools. The GIP configuration files (Supplementary Data 1) and the giptools command options used to generate all results (Supplementary Data 2) are provided. The full documentation of GIP and giptools including a description of all options is available from https://gip.readthedocs. io/en/latest/.

Read alignment
The repetitive elements of reference genomes were softmasked by GIP using Red (39). WGS reads were mapped by GIP using BWA-mem (version 0.7.17) (40,41) run with option -M to label shorter split hits as secondary. Then the alignment files were sorted, indexed and reformatted by GIP using Samtools (version 1.8) (42). Finally, read duplicates were removed by GIP using Picard MarkDuplicates (http://broadinstitute.github. io/picard) (version 2.18.9) with the option 'VALIDA-TION STRINGENCY = LENIENT'. In the four considered datasets, WGS reads were aligned against full assemblies, including unsorted contigs if present. However just the canonical assembled chromosomes were considered for all downstream analyses ('chrs' option, Supplementary Data 1). A minimum read alignment MAPQ score was adopted to select genes for cluster analysis, and to call for SNVs and SVs ('MAPQ' option, Supplementary Data 1). Altogether a total of 6 306 951 266 reads were aligned to the respective reference genomes. The 'giptools overview' module was run to gather the alignment statistics as estimated by Picard CollectAlignmentSummaryMetrics (Supplementary Table S2).

Genomic bins and genes quantification
GIP was used to evaluate the mean sequencing coverage and the mean read MAPQ of genomic bins and genes. For genomic bins, GIP partitioned the input genomes into adjacent windows of user defined lengths ('binSize' option, Supplementary Data 1). The coverage GC-content score bias was corrected ('CGcorrect' option, Supplementary Data 1) fitting a LOESS regression with a 5-fold cross validation to optimize the model span parameter. A larger window length was utilized to bin the reference genomes for circos plot representations ('binSizeCircos' option, Supplementary Data 1). In Figure 1 Figure S7 bins and genes coverage scores were normalized by median genome coverage to account for sequenc- Figure 1. GIP workflow. The schema on the left recapitulates the GIP inputs, processes and outputs (see Materials and Methods). Blue, orange, green and purple boxes indicate genome reference pre-processing, read mapping, quantification and variants computation modules, respectively. The panels on the right demonstrate example plots included in the GIP report computed for individual samples. The 'Karyotype' plot shows the coverage distributions for each chromosome (y-axis). The 'Genomic bins' plot shows the genomic position (x-axis) and the normalized genomic bin sequencing coverage (y-axis). The 'Gene CNVs' plot shows the normalized gene sequencing coverage. The 'Structural variants' panel shows a Circos plot representing translocations (black links in the inmost part of the plot), and other possible structural variations in the outer tracks, including insertions, duplications, deletions and inversions. The outmost track shows the normalized sequencing coverage. The 'Single nucleotide variants' plot shows on the x and y axes respectively the genomic position and variant allele frequency of detected SNVs.
ing library size differences. GIP evaluated statistically significant copy number variant bins and genes ( Figure 1 'Genomic bins' and 'Gene CNVs' plots) using a P-value threshold of 0.001 ('covPerBinSigOPT' and 'covPerGeSigOPT' options, Supplementary Data 1). Estimated P-values for bins and genes CNVs were corrected for multiple testing using the Benjamini -Yekutieli ('-padjust BY') and the Benjamini -Hochberg ('-padjust BH') methods. The somy scores shown in Figure 1 ('Karyotype' plot) and Figure  5A were computed multiplying the median genome coverage normalized bin coverage by 2. GIP enabled the CNV analysis of genes sharing high sequence identity by clustering the nucleotide sequences of the genes with low mean MAPQ score into groups with cd-hit-est (version 4.8.1) (43) with options '-s 0.9 -c 0.9 -r 0 -d 0 -g 1'. Then for each gene cluster GIP computed the mean gene coverage normalized by median chromosome coverage ( Figure 2E, Supplementary Figure S3B). The predictions of CNV regions returned by GIP for the P. vivax dataset were compared to the ones available from previously published work (44). For the analysis we utilized the 178 available samples not affected by national export restrictions. The sample list is available from https://www.malariagen.net/sites/default/ files/PvGV May2016 sample data 2.xlsx. To enable the comparison, GIP was re-executed considering the same genome refence (Pvax Sal1) and CNV length cutoff (>3 kb) used in the published work.

Gene ontology and metabolic pathway enrichment
The FungiDB online tool (Release 52, 20 May 2021) (45) was used to evaluate the functional enrichment of network gene clusters. For the gene ontology analysis, the biological process (BP), molecular function (MF) and cellular compartment (CC) terms enrichments were tested, considering both computed and curated evidences and a P-value cutoff of 0.05. For the metabolic pathway enrichment, both KEGG (46) and MetaCyc (47) pathway sources were considered with a P-value cutoff of 0.05. Terms and pathways with Benjamini -Hochberg adjusted P-values <0.05 were considered statistically significant.

Sequencing coverage density estimates
GIP was used to convert the read alignment files (.bam format) in binary data files reflecting sequencing coverage (.bigWig format). The coverage file were produced using bamCoverage from the deepTools2 suite (48) (version 3.5.1) with options '-normalizeUsing RPKM -ignoreDuplicates -binSize 10 -smoothLength 30' ('bigWigOPT' option, Supplementary Data 1). The coverage track of sample PD0689 C was visualized with IGV using the 'Bar Chart', 'Autoscale' and windowing function 'Mean' options.

Single-nucleotide variant analysis
GIP was used to call SNVs using Freebayes (version 1.3.2) ('freebayesOPT' option, Supplementary Data 1) and filter its output ('filterFreebayesOPT' option, Supplementary Data 1). Filters included the minimum allele frequency ('-minFreq'), the minimum number of reads supporting the alternative alleles ('-minAO') and minimum mean mapping quality of the reads supporting the reference ('-minMQMR') or the alternative allele ('-minMQM'). A higher number of reads supporting the variants was requested for predictions positioned inside simple repeats of the same nucleotide (homopolymers) ('-minAOhomopolymer'). The homopolymers were defined as the DNA region spanning ±5 bases from the SNV ('-contextSpan 5'), with over 40% of identical nucleotides ('-homopolymerFreq 0.4'). Further, GIP discarded SNVs with sequencing coverage above or below 4 median absolute deviations (MADs) from the median chromosome coverage ('-MADrange'). SnpEff (version 4.3t) (49) was used to predict the impact of SNVs on coding sequence. The predicted effects that GIP considered synonymous mutations are: 'synonymous variant', 'stop retained variant' and 'start retained'. The predicted effects that GIP considered non-synonymous mutations are: 'missense variant', 'start lost', 'stop gained', 'stop lost' and 'coding sequence variant'. The phylogenetic tree was computed by the giptools module 'phylogeny' using IQtree2 (version 2.1.2) (50,51) with options '-seqtype DNA -alrt 1000 -B 1000'. The Venn-diagram comparison considered the strains QS0044 C, QS0001 C, QS0037 C, QS0016 C and SGH 358 that were sampled from different locations in Ethiopia, respectively Habala, Badowacho, Arbaminch, Hawassa and Jimma. The strains were selected to have comparable average genome coverage (52). To infer the tree GIP considered the set of filtered SNV and adopted the IUPAC ambiguous notation for the positions with allele frequency <70%. The tree was visualized by giptools using the ggtree R-package (53).

Analysis of structural variants
GIP was used to detect structural variants including insertions, tandem duplications, deletions, inversions and translocations with DELLY (version 0.8.7) (28). SVs predictions were performed on individual samples sepa-rately. To reduce incorrect predictions the DELLY output was additionally filtered ('filterDellyOPT' option, Supplementary Data 1). GIP discarded poor predictions with DELLY label 'LowQual' ('-rmLowQual') and low median MAPQ score of mapping reads ('-minMAPQ'). SVs positioned in proximity of chromosome ends were removed ('-chrEndFilter') to limit false predictions caused by potential misassembled regions close to the telomeric ends. To ease visualization and limit the analysis only to best supported SVs GIP restricted the output only to the top predictions ('-topHqPercentIns', '-topHqPercentDel', '-topHqPercentDup' and '-topHqPercentInv') based on the SV support score as in Formula (1), where DV, DR, RV and RR are respectively the number of high-quality variant pairs, reference pairs, variant junction reads and reference junction reads.

Comparison of CNV callers
GIP was compared to CNVnator (version 0.  Table S1). To run CN-Vnator the commands '-tree', '-his', '-stat', '-partition' and '-call' were used. Following the authors' recommendations, the 'binSize' parameter was benchmarked using the '-eval' command optimizing the ratio of average read depth signal to its standard deviation. A common optimal 'binSize' value of 150 was selected considering the read length and the read depth of the sequencing libraries. Cn.MOPS was run using a windows length (WL) of 150. The WL value was chosen accounting for read depth differences between samples to have approximately an average of 50-100 reads per segment. To account for ploidy differences between samples cn.MOPS was run applying separate normalizations for each chromosome. The estimated chromosome somy score rounded to its closest integer was used as normalization factor with the 'normalizeChromosomes' function.

GIP and giptools running time and computational resources
The '-with-timeline' Nextflow option was used to render HTML timelines for all GIP executed processes. The '-withreport' Nextflow option was used to summarize GIP computational resources. A detailed description of each of the two options is available from Nextflow documentation at https://www.nextflow.io/docs/latest/tracing.html#. The giptools computational requirements were estimated using the Linux 'time' command (URL https://man7.org/linux/manpages/man1/time.1.html) together with the '-v' option to obtain a verbose output.

The GIP workflow
GIP is a tool for scientific investigation compatible with Linux systems, requiring minimal configuration and dis-

PAGE 5 OF 16
Nucleic Acids Research, 2022, Vol. 50, No. 6 e36 tributed as a self-contained package. Our integrative approach provides a broad range of genomic data analyses and visualizations, and combines new and existing bioinformatic methods (see Materials and Methods). GIP consists of three files: the Nextflow pipeline code, the configuration and the Singularity container files. The pipeline container conveniently provides a working environment with 19 offthe-shelf applications for genomic analyses and 40 R packages listed at https://gip.readthedocs.io/en/latest/software/ index.html. The minimum required input is a paired-end WGS data set and a reference genome assembly in the standard fastq and FASTA formats, respectively. GIP analyses include (i) extracting genomic features such as assembly gaps or repetitive elements, (ii) mapping the reads, (iii) evaluating chromosome, gene and genomic bin copy numbers, (iv) identifying and visualizing copy number variation with respect to the reference genome, (v) identifying and quantifying gene clusters, (vi) detecting and annotating SNVs, (vii) measuring non-synonymous (N) and synonymous (S) mutations for all genes, (viii) detecting SVs including tandem duplications, deletions, inversions and breakends translocations using split-read and read-pair orientation information and (ix) producing a report file providing summary statistics, tables and visualizations ( Figure  1, Supplementary Figure S1). For SNVs, calling GIP relies on Freebayes, a popular tool that resolves the issue of having multiple potential ambiguous alignments between the read and the homologous genomic region by examining the whole haplotype of the read independently of the precise alignment positions. It was shown that compared to alternative variant calling software Freebayes demonstrates good performance across different aligners and Illumina platforms (56). For SV detection GIP runs DELLY, a method combining short insert paired-ends, long-range mate-pairs and split-read alignments to precisely define balanced and unbalanced forms of genomic rearrangements at single-nucleotide resolution. DELLY's prediction accuracy was previously validated by PCR (28) and its performance compared to other structural variant calling algorithms on simulated data under different sequencing parameter settings (28). The benchmark results indicate that DELLY possess very high-positive predictive value and a robust performance across the simulated sequencing parameter space (28). GIP allows to customize filtering and visualization options via the configuration file (see methods). The output of GIP can be used as input for giptools, a tool-suite to compare sample sub-sets and highlight chromosome copy number, gene copy number and SNV differences.

Applying giptools on a Leishmania infantum case study
GIP permits the batch analysis of a set of individual samples, where each sample is considered separately and compared only with respect to the provided reference genome assembly. As a consequence, all variants and copy number alterations detected in a sample merely reflect the differences between the sequenced and the reference genomes. While this application may be sufficient in some circumstances, research projects often involve downstream comparison between samples. Examples include the comparison of gene or chromosome copy variation number between  (Table 1). All modules in giptools are fully embedded in the Singularity container and they are provided with their own documentation.
To illustrate the type of exploratory data analyses and the biological questions that can be addressed, we tested giptools on a previously analyzed dataset of seven clinical L. infantum isolates from Tunisia (31). Leishmania is the etiological agent of leishmaniasis, a life-threatening human and veterinary disease affecting 12 million people worldwide (57). Parasites were derived from seven patients affected with visceral leishmaniasis, expanded in cell culture and their genomes were sequenced. This dataset includes four Glucantime drug susceptible isolates and three isolates from relapsed patients, and their comparison may inform on genetic factors resulting in treatment failure. Giptools allowed the detection and visualization of pervasive intra-chromosomal CNVs across the thirty-six Leishmania chromosomes ( Figure 2A). Additionally, giptools enables targeted comparison of normalized genomic bin sequencing coverage of sample pairs. We used giptools' 'binCNV' module to compute the ratio between corresponding genomic bins of the strains LIPA83 over ZK43, which correspond to a first-episode and a relapse leishmaniasis isolate, respectively. Giptools represents different chromosomes as separate panels (Figure 2B), as part of single genomewide overview (Supplementary Figure S2A) or as distinct plots (Supplementary Figure S2B). This analysis allowed the identification of 2,905 and 2,208 bins that were respectively amplified or depleted in LIPA83 with respect to ZK43. The results are returned by giptools as a Microsoft Excel table (.xlsx format) providing ratio scores at each genomic position (Supplementary Table S3). Likewise, giptools permits three-way comparisons of normalized genomic bin sequencing coverage with ternary plots  ( Figure 2C). We used this representation to display the genomic bin relative abundance in samples ZK43, LIPA83 and ZK28. The analysis shows important strain-specific differences in bin copy number that are visualized by shifts of the signals out of the center. Similar to genomic bin analysis, giptools makes it possible to compare the sequencing depth of annotated genes and thus the copy number in two or three samples ( Figure 2D and Supplementary Figure S3). However, the determination of gene copy number might be impeded by (i) short read length or fragment insert size, (ii) the complexity of the target genome, or (iii) the presence of repetitive elements. Together with the coverage, GIP also computes a mean read map quality (MAPQ) score for each gene and allows a different strategy to determine the copy number of low MAPQ genes (see methods). This gene MAPQ score is a measure that reflects how much each gene is supported by unambiguously mapped reads (high MAPQ) in contrast to multimapping reads (low MAPQ). The evaluation of LIPA83 and ZK43 gene coverage ratio scores revealed 13 gene CNVs with a stringent MAPQ cutoff of 50 (Supplementary Table S4). The maximum normalized coverage value of 6.5 was observed for a putative amastin surface glycoprotein (LINF 310009800). Other examples of gene CNVs in this set include the putative surface antigen protein 2 (LINF 120013500) and the heat shock protein HSP33 (LINF 300021600) (Supplementary Table S4). Genes falling below a user defined MAPQ score and sharing high level of sequence similarity are assigned to the same gene cluster, and their measured coverage scores are averaged across all members of the group. Low MAPQ scores can also be associated with single genes, e.g. in the case of internal repetitive elements that cause multiple ambiguous alignments inside the gene itself, or if mapping occurs in possibly misannotated intergenic regions. The LIPA83/ZK43 comparison showed 27 CNV gene clusters, including cluster cl303 (three genes annotated as 'amastin-like') and cluster cl16 (2 tb-292 membraneassociated protein-like proteins) ( Figure 2E, Supplementary Table S4). These results demonstrate the power of GIP and giptools to detect and compare intra-chromosomal CNVs in Leishmania at genomic bin level. Conveniently, analogous two-or three-ways comparisons can be applied to reveal copy number variations at individual gene or gene cluster levels.

Benchmarking of prediction tools
GIP uses Freebayes and DELLY to perform SNVs and SVs calls while providing additional options to filter and reformat their output for downstream analyses. Freebayes and DELLY predictions were already extensively benchmarked against other tools and validated in previous reports (28,56,58). GIP implements its own approach to quantify, normalize, compare and visualize genomic bins and genes. One gene CNV predicted using GIP was experimentally tested in a separate study where we confirmed by PCR the deletion of a NIMA-like kinase gene (59). In the following we used the L. infantum dataset to compare the gene CNV predictions returned by GIP with the ones of two popular CNV detection tools, CNVnator (27) Figure S4A). Overall CNVnator and cn.MOPS predict 13.4 and 3.2 times more CNV genes than GIP, respectively. While this result could reflect a superior sensitivity, it suggests a high rate of false positive predictions caused by confounding chromosome aneuploidy. The chromosomes with the greatest number of CNVnator and cn.MOPS predictions are the ones which are amplified in at least one of the considered isolates (Supplementary Figure S4B). CNVnator can be run on individual samples only, therefore it is not suited to account for potential between-samples chromosome copy number differences. cn.MOPS can be run on multiple samples together and allows the possibility to apply adapted normalizations on each chromosome separately. This approach alleviates the problem of aneuploidy-induced false gene CNV predictions, demonstrating to be effective for some chromosomes (e.g. LinJ12 or LinJ13, Supplementary Figure S4B) but still failing for chromosomes LinJ06, LinJ16, and partially for LinJ23 and LinJ24. In cn.MOPS the normalization factor is an integer value representing the chromosome ploidy. Conceivably, the normalization in cn.MOPS suffers from the within-sample population mosaicism, in which individual cells present chromosome copy number differences. On the contrary, to call bin or gene CNVs GIP normalizes by the actual measured median chromosome sequencing coverage, thus accounting for 'partial' chromosome copy number shifts. Overall, GIP uses a maximum of 1.057 Gigabytes of random access memory (RAM) in the 'covPerGe' step to calculate normalized gene sequencing coverage values, plus 0.39236 Gigabytes in the 'geInteraction' giptools step to evaluate the copy number variant genes. CNVnator and cn.MOPS require 2.727132 and 0.791924 Gb respectively. The cumulative running time to execute the GIP 'covPerGe' step on each individual sample and to execute the 'geInteraction' giptools step is of 18 min and 59.06 s. For comparison CNVnator employs 12 min and 46.53 s, while and cn.MOPS takes 6 min and 21.04 s. Altogether, these results show that to compute gene CNVs on the L. infantum dataset our pipeline has similar memory requirements but longer execution times compared to CNVnator and cn.MOPS. The analysis confirms that the CNV gene predictions produced by our pipeline largely recapitulate the ones produced by CNVnator or cn.MOPS. Furthermore, GIP provides the added benefit of a higher robustness when predicting intrachromosomal CNVs in the event of chromosomal aneuploidy.

Comparative genomics of a Plasmodium vivax WGS dataset
We next applied GIP and giptools on other biological systems to demonstrate its broad applicability outside the Leishmania field, including the human apicomplexan parasite P. vivax. Plasmodium vivax is a protist parasite and a human pathogen causing malaria. Plasmodium vivax gives rise every year to 130 million clinical cases (60), and it is estimated that 2.5 billion people are at risk of infection worldwide (61-63). We applied GIP and giptools to inves-  Table S1). The GIP quantification analysis was able to recover 5 out of the 11 previously reported CNV regions and predicted 73 new large (>3 kb) copy number variant areas (Supplementary Figure S5, Supplementary Table S5). Unmatched CNVs could be explained by sample sets differences. The published study considered a set of 228 samples from which two unspecified samples were removed showing excessive variation in read coverage. Conversely, in our analysis we considered the set of 178 available samples not affected by national export restrictions. The phylogenetic tree reconstruction and PCA analyses ( Figure 3A and B) showed a high correlation between genotypes and the geographic origin of the samples. However, we detected substantial genomic variability between isolates collected at smaller geographical scale, with 14 555 SNVs (∼42% of the total) uniquely characterizing representative samples from five Ethiopian study sites (52) (Figure 3C and D). This result may reflect diverging evolutionary trajectories radiating from few founder strains. At gene level we profiled the copy number variations of two gene panels. The first panel accounts for 43 previously described genes encoding for potential erythrocyte binding proteins suggested to operate at the interface of the parasitehost invasion process (52). The second panel includes two drug resistance markers comprising the chloroquine resistance transporter PVP01 0109300 and the multidrug resistance protein 1 PVP01 1010900, and four proteins implicated in red blood cell invasion, such as the merozoite surface protein gene MSP7 PVP01 1219700, the reticulocyte binding protein gene 2c PVP01 0534300, the serine-repeat antigen 3 PVP01 0417000, and the reticulocyte binding protein 2b PVP01 0800700) (64)(65)(66)(67)(68)(69)(70)(71)(72)(73). Read depth analysis indicated that four genes in the first panel (PVP01 0623800, PVP01 1031400, PVP01 1031200, PVP01 1031300) show a high degree of variability, with amplifications observed in samples from distinct geographic ( Figure 3E). This convergence is sign of strong natural selection, which further sustains the functional importance of these genes in the infection process. Furthermore, six genes positioned on chromosome 14 are absent in the Thai strain PD0689 C as a result of the loss of this chromosome (nullisomy) (Figure 3E and F). Finally, the comparison of synonymous and nonsynonymous SNVs in the panel of genes revealed important differences between sample groups. Our analysis indicates an overall higher number of non-synonymous mutations in Ethiopian compared to Cambodian isolates, therefore suggesting a stronger evolutionary pressure acting on the African strains ( Figure 3G). Taken together these analyses well illustrate how GIP and giptools can be readily applied for bulk analysis of P. vivax genomes to assess genome diversity, extract evolutionary information and identify potential disease biomarkers.

Gene CNV analysis of Candida albicans evolutionary adapted strains
We next applied GIP and giptools to the human fungal pathogen C. albicans, an opportunistic yeast exhibiting major genome plasticity (74-83) and causing hun-dreds of thousands of severe infections each year (84). Candidemia, a bloodstream infection with Candida, are often associated with high rates of morbidity and mortality (15-50%) notwithstanding existing antifungal treatments (85,86). We applied GIP and giptools to a C. albicans WGS dataset described in a recent study that covers five different progenitor strains (P75063, P75016, P78042, SC5314, AMS3050) and investigates CNVs driving tolerance and resistance to anti-fungal azole drugs (87) (Supplementary Table S6). We analyzed nineteen samples, including (i) four clinical isolates (P75063, P75016, P78042, SC5314), (ii) seven strains selected in vitro against the anti-fungal drug fluconazole (FLC) (AMS4104, AMS4105, AMS4106, AMS4107, AMS4397, AMS4444, AMS4702), (iii) four isogenic colonies adapted to the drug miconazole (AMS3051, AMS3052, AMS3053 and AMS3054) together with their progenitor (AMS3050) and (iv) three colonies derived from a miconazole-adapted population and isolated on a rich medium (AMS3092, AMS3093 and AMS3094) (87-89). GIP and giptools were able to reproduce previous observations of the amplification of the genes for the drug efflux pumps TAC1 (orf19.3188) and ERG11 (orf19.922), for the stress response proteins HSP70 (orf19.4980), CGR1 (orf19.2722), ERO1 (orf19.4871), TPK1 (orf19.4892), ASR1 (orf19.2344), PBS2 (orf19.7388) and CRZ1 (orf19.7359), and for proteins involved in membrane and cell wall integrity, including CDR3 (orf19.1313), NCP1 (orf19.2672), ECM21 (orf19.4887), MNN23 (orf19.4874), RHB1 (orf19.5994) and KRE6 (orf19.7363) (Supplementary Table S6). Furthermore, the powerful comparative approach of our pipeline permitted the discovery of 1505 genes showing correlating or anticorrelating copy number variations ( Figure 4A and B, Supplementary Table S6), which could be assigned to nine distinct correlation clusters (CC) (Supplementary Figure S6A Figure S6D) (87). Eventually, by representing genes and absolute correlation respectively as nodes and edges of a network, we identified nine highly interconnected network clusters (NC) ( Figure  4C, Supplementary Table S8). NC7, NC8 and NC9 embody genes from individual chromosomes, respectively chromosomes 1, 3 and 4. The most parsimonious explanation for the high levels of correlation observed in these NCs (Figure 4C) is the occurrence of sub-chromosomal amplifications affecting several adjacent genes. A different scenario is pictured for each of the remaining NCs (NC1-6) where the genes are located on different chromosomes thus suggesting genetic interactions that causes coordinated changes in gene copy number. The gene ontology (GO) and metabolic pathway analyses revealed a significant functional enrichment of genes expressed on the cell surface and involved in the interaction with the host (NC2), gibberellin biosynthesis (NC3), transmembrane nucleobase transporters (NC4) and gluco- Comparison of non-synonymous (N) and synonymous (S) mutations between Ethiopia and Cambodia sample groups. Dots represent genes. The x-axis represents the difference between the mean non-synonymous mutation count in the two sample groups. The y-axis represents the difference between the mean synonymous mutation count in the two sample groups. The dot size demonstrates the ratio of the mean normalized sequencing coverage between the two sample groups for each gene. Red and blue dot colors indicate genes belonging to the 43 genes panel (52) and the custom 6 genes panel respectively.  Table S9). Altogether, GIP and giptools are validated by reproducing previously published results, and beyond that can drive new biological findings as documented by the discovery of a network of epistatic CNV interactions supporting genomic adaptation in C. albicans populations under drug selection.

Exploring instability of larger genomes using cancer cell lines as a benchmark
The larger genome size, and the higher number of genes and WGS reads can represent a challenge when working with higher eukaryotes. For the purpose of comparison, the human genome is ∼216 times larger than the one of C. albicans we analyse in this study. Therefore, we sought to evaluate the applicability of the GIP and giptools framework to human data and utilized a panel of genomes from cancer cell lines as a test set. In our analyses we considered publicly available WGS data of the cell lines T47D, NCI H460 and K562 (90), which respectively derive from human breast, lung and blood cancers. The karyotype analysis revealed aneuploidy for all chromosomes except chromosome 4 ( Figure 5A). The observed heterogeneity in read depth across chromosomes, illustrated by large interquartile range in the boxplot, suggests sub-chromosomal or episomal copy number variations, or the co-existence of karyotypically different sub-populations. Indeed, the coverage analysis confirmed the pervasive presence of CNVs both at chromosomal and sub-chromosomal levels (Supplementary Figure S7) with remarkable instability observed for specific chromosomes, e.g. chromosomes 6, 9, 10 and 16 ( Figure 5B). Overall, we detected 1 647 016 SNVs (Supplementary Data 3) and allele frequency shifts with respect to the reference genome, suggesting haplotype selection and the preferential expression of distinct alleles in different cell lines ( Figure 5C and Supplementary Figure S8). Furthermore, we identified repeated loss of heterozygosity events and uneven distribution of SNVs that form 'patches' of high frequency correlating with chromosomal and sub-chromosomal CNVs ( Figure  5D, E and Supplementary Figure S9). These results identify GIP and giptools as a powerful new platform to reveal loci, genes or alleles that are under natural selection in cancer cells, thus allowing important new insight into the genetic basis of tumor development, cancer cell evolution and drug resistance.

Running times and computational resources
We used the Nextflow built-in options to render the process execution time and the computational resources of each task executed by GIP for each sample considered in this study. As exemplified for the cancer cell line dataset, the GIP read mapping process accounts for the longest execution time for all three samples (Supplementary Figure S10A) and the highest I/O (Input/Output) data access, with peaks of read and written bytes of 773.4 and 598.8 Gb, respectively (Supplementary data 4). The process requiring the highest physical memory is 'bigWigGenomeCov' (mean usage 37.18 Gb) (Supplementary Figure S10B). This step runs deepTools2 and bedGraphToBigWig (91) to create wiggle (wig) type files representing continuous-valued sequencing coverage data in indexed binary format (bigWig). The time line and the computational resources usage including CPU, memory, job duration and I/O of all datasets are provided (Supplementary data 4). The computational resources and process execution statistics to run all giptools commands described in this study are provided (Supplementary Table  S10). The longest running time (474 578.57 and 4 417.61 s in user and system time, respectively) and the maximum required RAM size (9 465 100 Kb) were measured for the estimation of P. vivax phylogeny (Supplementary Table S10), and largely explained by the execution of IQtree2.

DISCUSSION
Genome instability is a key driver of evolution for microbial pathogens and cancer cells (92) and a major source of human morbidity. Here we introduce GIP and giptools, an integrated framework for the genotype profiling of biological systems exploiting genome instability for adaptation. While our pipeline relies on existing genome analysis tools (Supplementary Figure S1), it also implements its own algorithm that quantifies, corrects by GC content, normalizes, compares, estimates significance and visualizes coverage of genomic bins and genes. Furthermore, it defines and quantifies gene clusters, detects chromosome aneuploidies and discovers gene interactions. We document the power and versatility of GIP and giptools by performing genomic screenings of three major pathogenic eukaryotes and human cancer cell lines. While originally deployed for Leishmania genome analysis, in this study we validate the use of our pipeline on other organisms reproducing expected results. For example, in C. albicans we confirmed the CNVs correlating to drug resistance as well as a loss of heterozygosity event (Supplementary Figure S6B-D). Parallel to this we also show how GIP and giptools can be used for data mining and discovery of new signals that escaped existing tools. New findings include (i) the discovery of the convergent amplification of erythrocyte binding proteins in P. vivax strains sampled from distinct geographic areas ( Figure 3E), (ii) the detection of a nullisomic strain ( Figure 3F), (iii) the identification of correlated copy number variations between genes positioned on separate chromosomes of C. albicans adapting strains, and (iv) the functional association of such genes, strongly supporting a mechanism of epistatic interactions exerted through gene-dosage changes, and corroborating previous reports on adapting Leishmania populations (59).
Importantly, GIP and giptools overcome key limitations of current analysis tools, such as the breadth of analysis that is often limited to individual types of mutations, and the lack of genome-wide, comprehensive reports. To ease genome instability investigations our pipeline offers a single solution to karyotype, gene CNV, SNV and SV batch analyses, providing summary reports and high-quality, genomewide visualizations. Furthermore, some tools (27,(93)(94)(95) identify variations with respects to a reference assembly only, which leaves the between samples comparisons to external tools that need installing, may be incompatible in terms of file format, and may rely on different analytical assumptions. To address this limitation giptools enables custom sample comparisons, to explore differences and common features between genomes, and provides a vast choice of analytical tools with compatible features. Likewise, current tools are often restricted to the analysis of data from one or few species only (96-100), but are not generally applicable to different biological systems, which interferes with the investigation of genome variations across multiple species and the exploration potentially conserved genomic adaptation mechanisms. By contrast, our pipeline limits as much as possible the use of hardcoded parametrization, which could limit its use to a specific organism. Therefore, GIP's flexible design makes it adapted for the genome analysis of both model and non-model organisms, including Leishmania or human.
Many current tools are further limited in software portability and reproducibility across different computer environments, which can produce faulty results calling their clinical application into question. Conversely, thanks to the Singularity implementation all required software are embedded and provided within the software container. As a consequence, users can easily recreate the same work environment just by downloading the pipeline container, and reproduce exactly the same publication-quality plots and tables presented in this study. Lastly, one more common limitation is posed by software scalability. In the WGS domain, with the rapid increase new samples made available and the enormous amount of data generated in each sequencing run, the CPU and memory resources of local workstation risk to quickly become inadequate for data analysis. Therefore, it is paramount that WGS tools are implemented to run on high-performance computing (HPC) clusters and feature remote cloud computing solutions. Because of its Nextflow implementation GIP can be executed on a local machine, on cluster resource manager or the cloud. GIP can be applied on individual samples and without additional effort on large WGS data sets for batch computation as shown for the 222 P. vivax genomes.
These results well illustrate how GIP and giptools can be applied to perform extended genomic analyses in different biological systems and drive biomedical discovery. To conclude, we believe that GIP and giptools represent a step forward toward reproducible research in genomics, and provide a robust computational framework to study how microbes and tumor cells harness genome instability for environmental adaptation and fitness gain.

DATA AVAILABILITY
All data is available in the main text or the supplementary materials.

SUPPLEMENTARY DATA
Supplementary Data are available at NAR Online.