Computational protocol: Beaver Fever: Whole-Genome Characterization of Waterborne Outbreak and Sporadic Isolates To Study the Zoonotic Transmission of Giardiasis

Similar protocols

Protocol publication

[…] The quality of the reads was assessed by Fastqc (, and the adapter sequences and the sequences of poor quality were trimmed by Trim Galore ( Reads were assembled using assembler SPAdes v.3.1.1 () with default parameters, and the resulting contigs were filtered to be >500 bp and have a coverage of >4. We assessed the assemblage type of each isolate by mapping reads from each isolate to the reference genomes of G. duodenalis (assemblages A, B, and E) that were retrieved from the GiardiaDB ( (). To remove any potential contaminant contigs (non-Giardia reads), contigs were mapped to the NCBI nt database using BLASTn megablast (-best_hit_overhang 0.25 -best_hit_score_edge 0.05). If 80% of hits within two times the minimum E value mapped to G. duodenalis genomes, the contigs were kept. Additionally, any contigs that did not fit this criterion were mapped to a combined reference of assemblage A (WB) and assemblage B (GS_B) using BLASTn megablast with the same parameters, and those for which >500 bp matched the reference were added to the final contig set.Even after removing contaminating reads, five samples had a total contig length greater than 14 Mb. This is significantly larger than the expected genome size of 10 to 11 Mb, suggesting further contamination from within Giardia species. Reads of these five samples were each mapped to a combined reference of assemblages A and B (WB, DH_A, and GS_B) and using BLASTn megablast, with the same parameters. In mixed samples, reads which hit only the minority assemblages (3.5 to 28.7% of total reads) were removed. SPAdes assembly and subsequent filtering were then performed on these filtered reads as described above. This resulted in two of the five samples having a combined contig length of <12,000,000 bp; however, three samples were still above this length and were excluded from analysis, as we could not separately the contamination from the major assemblage. Final contigs from each sample, along with the references WB and GS_B, were used to construct two assemblage specific pan-genomes using Panseq () with default parameters. The assembled contigs of each sample were arranged and compared using MUMmer (). [...] Reads were trimmed using cutadapt v1.8 () to remove Illumina adapters and reads with quality of <20 or length of <200. Assemblage assignment was then performed by mapping trimmed reads to a combined reference consisting of assemblage A (WB) and assemblage B (GS_B) using bowtie2 v2.1.0 () (with parameters --phred33 --local --dovetail --maxins 850) and determining the assemblage based on what the majority (70 to 99%) of reads mapped to. Trimmed reads were then remapped to the single majority assemblage using bowtie2 (--local).Four samples had >5% of reads mapped to the minor assemblage after the initial bowtie2 alignment, so these were included in the analysis for both assemblages. For analysis in their minor assemblage, reads were mapped to a combined reference of WB, DH_A, and GS_B, using BLASTn megablast, with the same parameters as described above, and reads which hit only the majority assemblages were removed. The remaining reads were then included in the discoSNP analysis as below.Single nucleotide variant (SNV) calling was performed by entering trimmed reads into discoSNP v2.1.2 () with parameters -c 4 –b 1 –k 31 –P 3, and a separate run for each assemblage. SNVs generated using discoSNP are output in “bubbles,” where each bubble consists of two distinct sequences (or branches) of k + 2 nodes, with the start and the end nodes in common. Each bubble was aligned to the assemblage-specific pan-genome using BLASTn and filtered so that only bubbles that hit a single unique position in the pan-genome were considered for further analysis. Additionally, bubbles were removed if any sample had a combined coverage over both branches of the bubble of <4 or any sample had <75% coverage in the majority branch. Filtered discoSNP output was then converted into a FASTA file using the SNV base from the branch with the highest coverage from each bubble. Phylogenetic trees were generated using FastTree v2.1.8 (default parameters) () and drawn using FigTree v1.4.0 ( The phylogenetic relationships among isolates were also inferred using SplitsTree, and a Φ-test was conducted to estimate the level of recombination (). The correlation between genetic and geographic distances was established using MEGA () and GenAlEx (). […]

Pipeline specifications