Overview: Using samples from amplicon deep sequencing from gene targets, either call haplotypes or call variants.

Method: Targeted amplicon deep sequencing which produces forward and reverse fastq files for each sample.

Snakemake is a workflow manager that enables massively parallel and reproducible analyses. Snakemake is a suitable tool to use when you can break a workflow down into discrete steps, with each step having input and output files.

We provide this example workflow as a template to get started running the pipeline with Snakemake. To adjust to your specific data, you can customize the config.yaml file.

For more details on Snakemake, see the Snakemake tutorial .

The Workflow

The Snakefile contains rules which define the output files we want and how to make them. Snakemake automatically builds a directed acyclic graph (DAG) of jobs to figure out the dependencies of each of the rules and what order to run them in.

The DAG shows how the pipeline can call haplotypes and call variants simultaneously, and how calls can be run in parallel if Snakemake is allowed to run more than one job at a time. Here, the pipeline is parallelized by being called with two targets each for haplotype and variant calling (AMA and CSP for haplotype calling, and PFCRT and PFMDR1 for variant calling) for and six quality scores (2, 5, 10, 15, 20, 25) at the same time.

call_fastqc uses FastQC to perform quality control on the raw reads and allow for visualization of the quality of the reads.

call_trimmomatic cleans and filters the raw reads. It uses CutAdapt to trim the primers and adapter sequences from sequencing reads and Trimmomatic to quality filter reads if the average of every 4 nucleotides had a Phred Quality Score < 15 or was less than 80 nucleotides long.

generate_multiqc_report uses the MultiQC reporting tool to parse summary statistics from results and log files generated by fastqc and cutadapt.

These first three rules are called for the given input data, regardless of whether the pipeline is being used to call haplotypes or call variants. Following this, the pipeline splits up into haplotype and/or variant calling.

variant_calling uses a package called Burrows-Wheeler Aligner, or BWA, which maps low-divergent sequences against a large reference genome. This is used for variant calling for drug resistance genes.

analyze_vcf takes in the vcf file(s) produced by the previous rule and returns a table of drug-resistance-associated mutations among the amplicon vcf file(s).

combine_vcf_analysis uses each target's table returned by analyze_vcf and combines them into a summary table that allows for easy analysis of the variant calling results.

get_marker_lengths calculates the length of each marker using the lengths of the reference and primer sequences.

synchronize_reads cleans, filters, and maps the raw reads. It uses BBmap to map all reads from the reference sequences, CutAdapt to trim the primers and adapter sequences from sequencing reads, and Trimmomatic to quality filter reads if the average of every 4 nucleotides had a Phred Quality Score < 15 or was less than 80 nucleotides long. This is the first step in read processing on the cluster.

trim_and_filter filters and trims the forward and reverse reads using the DADA2 program and outputs them into a filtered folder.

call_haplotypes calls haplotypes for your target(s) and writes the results out into a reads table.

optimize_reads finds the quality score that produces the highest number of read counts and uses only the filtered reads created using that quality score.

merge_read_tracking combines the read count files from the previous three rules into one table that shows the number of reads at different points of the workflow.

censor_haplotypes censors falsely detected haplotypes. Censoring criteria is applied in this order:

1. Haplotypes that occur in < 250 of the sample’s reads are removed. You can change this criteria in the config file ("read_depth").

2. Haplotypes that occur in < 3% of the sample’s reads are removed. You can change this criteria in the config file ("proportion").

3. Haplotypes that are a different length than the majority of haplotypes (300 nucleotides for pfama1, 288 nucleotides for pfcsp. You can change the length in the config file ("haplotype_length").

4. For haplotypes that have 1 SNP difference, occur in the same sample, and have a > 8 times read depth difference between them within that sample, removed the hapltoype with the lower read depth from that sample. You can change the 8 to a different read depth ratio in the config file ("read_depth_ratio").

5. If a haplotype is defined by a single variant position that is only variable within that haplotype, then it is removed.

get_read_summaries runs a bash script that pulls statistics from the summary files produced during the trim_and_filter rule for later manipulation in R.

create_summaries outputs csv files that summarizes the trim and read counts.

Quick Start

1. Clone or download this repo.

   git clone https://github.com/kathiehuang/haplotype_calling_pipeline.git
   

   Alternatively, if you're viewing this on GitHub, you can click the green Use this template button to create your own version of the repo on GitHub, then clone it.

2. Install the dependencies.

   1. If you don't have conda yet, we recommend installing miniconda .

      If you use a Macbook that has an M1 Chip, you will need to install Rosetta2:

      softwareupdate --install-rosetta
      

      Download the x86-64 bit version of miniconda instead of the M1 version.

   2. Next, install mamba , a fast drop-in replacement for conda:

      conda install mamba -n base -c conda-forge
      

   3. Finally, create the environment and activate it:

      mamba env create -f environment.yaml # you only have to do this once
      conda activate haplotype_calling # you have to do this every time
      

   • Alternatively, you can install the dependencies listed in environment.yaml however you like.

3. Edit the configuration file config.yaml .

   • haplotype_calling_targets : the list of target(s) you want to perform haplotype calling on.

   • variant_calling_targets : the list of target(s) you want to perform variant calling on.

   • pair1 : the path to the folder containing the forward reads.

   • pair2 : the path to the folder containing the reverse reads.

   • refs : the path to the folder containing reference sequences for the polymorphic gene target that will be used to map the raw reads to the appropriate gene targets of interest.

   • forward : the path to the file with the list of forward primers.

   • rev : the path to the file with the list of reverse primers.

   • variant_table : the path to the table that holds positions of interest for each target ran through variant calling.

   • root : the path to your desired output directory.

   • truncQ_values : the list of values of truncQ to be used in the filterAndTrim function to find the optimal truncQ value.

   • cutoff : the cutoff for which samples with less than this number of reads after sampling should be removed.

   • seed : the seed of R's random number generator for the purpose of obtaining a reproducible random result.

   • read_depth : the cutoff for which haplotypes that occur in less than this amount of the sample reads should be removed.

   • proportion : the cutoff for which haplotypes that occur in less than this percentage of the sample reads should be removed.

   • read_depth_ratio : cutoff for which haplotypes that have 1 SNP difference, occur in the same sample, and have a greater than this amount times read depth difference between them within that sample with the lower read depth should be removed

   You can leave these options as-is if you'd like to first make sure the workflow runs without error on your machine before using your own dataset and custom parameters.

4. Modify your adapters files to notify CutAdapt that you want to use linked adapters.

   A linked adapter combines a 5’ and a 3’ adapter, so we use this if a sequence is surrounded by a 5’ and a 3’ adapter and we want to remove both adapters. Linked adapters aren't required for this pipeline to run, but we found that for our data, using linked adapters produced better results.

   Linked adapters are specified as two sequences separated by ... (three dots), with the second sequence being the reverse complement of the reverse primer for the forward primer and the reverse complement of the forward primer for the reverse primer. We use -g instead of -a (5' vs 3' adapters) when calling CutAdapt, which causes both adapters to be required, even if they are not anchored. However, we want the non-anchored adapters to be optional. You can mark each adapter explicitly as required or optional using the search parameters required and optional . As a result, we add the optional parameter to the end of the sequence after the three dots. This is what our forwardPrimers.fasta file looks like for CSP:

   >Pfcsp-f
   TTAAGGAACAAGAAGGATAATACCA...CATTTCGGTTTGGGTCATTT;optional
   

   And our reversePrimers.fasta file looks like this:

   >Pfcsp-r_rc
   AAATGACCCAAACCGAAATG...TGGTATTATCCTTCTTGTTCCTTAA;optional
   

   If you choose to forgo using linked adapters, this is what your adapters file would look like:

   >Pfcsp-f
   TTAAGGAACAAGAAGGATAATACCA
   

   >Pfcsp-r_rc
   AAATGACCCAAACCGAAATG
   

   For more information about the linked adapters parameter, visit the CutAdapt documentation.

5. Do a dry run to make sure the snakemake workflow is valid.

   snakemake -n
   

6. Run the workflow.

   Run it locally with:

   snakemake --cores {NCORES} --stats output/stats # NCORES = number of cores, ie. without parallelization use snakemake --cores 1. You can use snakemake -n to see how many jobs you'll need given the number of targets and quality scores you put in the config file. The stats command produces a JSON file where you can see statistics about your run, such as runtime.
   

   To run the workflow on an HPC with Slurm :

   1. Edit your email address ( YOUR_EMAIL_HERE ) in:

      • scripts/submit_slurm.sh

   2. Edit the number of cores you want to use for this pipeline run (right now it is defaulted to 12 cores).

   3. Submit the snakemake workflow with:

      sbatch scripts/submit_slurm.sh
      

      Slurm output files will be written to logs/ . You will receive an email when the job is finished.

Usage Examples

Below shows what output files looked like after running this pipeline with a small sample of haplotypes.

call_fastqc should produce an output directory with a fastqc_out folder that contains FastQC reports for each input fastq file. This allows you to visualize and assess quality collectively across all reads within a sample. To open the .zip files, download and extract the files. A copy of the config.yaml file will also be added to the output folder so that the parameters for each run can be saved.

Example of FastQC output files for a single sample: BF1_S1_L001_R1_001_fastqc.html , BF1_S1_L001_R1_001_fastqc.zip , BF1_S1_L001_R2_001_fastqc.html , BF1_S1_L001_R2_001_fastqc.zip

call_trimmomatic should add a trimmed_reads folder within the output directory containing the results from the Trimmomatic program, which filters poor quality reads and trims poor quality bases from your samples and produces statistical summaries.

The trimmed_reads folder contains four folders labeled 1 , 2 , singleton , and summaries . As an example, the trim/1 folder contains a BF1.1.fastq.gz file, the trim/2 folder contains a BF1.2.fastq.gz file, the singleton folder contains BF1.1_unpaired.fq.gz and BF1.2_unpaired.fq.gz files, and the summaries folder contains a BF1.summary file.

generate_multiqc_report should create a file multiqc_report.html in the output directory. If you download this file to your local computer, you can open it in a browser to see the aggregated results as an HTML report.

Variant Calling

variant_calling should create a variant_output directory, as well as a directory within it for each target, with each one containing four folders: bam , sam , sort_bam , and vcf .

analyze_vcf should temporarily create a file called DR_mutations.csv in the directory for each target. After running the pipeline on a small set of samples, the PFCRT table looked like this:

combine_vcf_analysis should create a file called dr_depths_freqs.csv in the variant_output directory. After running the pipeline on a small set of samples with targets PFCRT and PFMDR1, our depths frequency table looked like this:

Haplotype Calling

get_marker_lengths should add a csv file containing the names of each marker and their corresponding lengths called marker_lengths.csv to the output directory.

synchronize_reads should create a haplotype_output directory, as well as a directory within it for each target. For each target, the rule creates a bbsplit_out folder containing three folders: mapped_reads , ref , and map_info . mapped_reads contains the cleaned and mapped fastq files, ref contains information about the indexes produced and referenced by the BBSplit program, and map_info contains summaries of the BBSplit call on each sample.

As an example, for AMA, the mapped_reads folder contains BF1_AMA_1.fastq.gz and BF1_AMA_2.fastq.gz files, the ref folder contains a genome/1 folder that has merged_ref_64917.fa.gz , namelist.txt , and reflist.txt files, and map_info contains a BF1.txt file.

trim_and_filter should produce a trim_filter_out folder in each target's directory that contains a summary for read trimming and filtering for each q value (which can be changed in the config file), {target}_{q_values}_trim_and_filter_table . It also produces a {target}/read_count file that lists the read counts for each q value.

After running the pipeline on a small sample of haplotypes, the AMA trim and filter table for a q value of 2 looked like this:

And the AMA read count table looked like this:

call_haplotypes should add temporary files {target}_{q_values}_haplotypes.rds to the haplotype_output folder as well as a file {target}_{q_values}_track_reads_through_pipeline.csv to the trim_filter_out folder.

• {target}_{q_values}_haplotypes.rds : R file that stores the haplotype results data set for each q value for further manipulation in censor_haplotypes .

• {target}_{q_values}_track_reads_through_pipeline.csv : tracks the reads, looking at the number of reads that made it through each step of the pipeline for each q value.

With our small sample, the AMA track reads through pipeline table for a q value of 2 looked like this:

optimize_reads should produce an optimize_reads_out folder with 3 files: {target}_final_q_value , {target}_max_read_count , and {target}_final_track_reads_through_pipeline , and remove the temporary files from the last step.

• {target}_final_q_value : contains the quality score that produced the highest number of read counts.

• {target}_max_read_count : contains the value of the highest number of read counts—the read count associated with the optimal quality score.

• {target}_final_track_reads_through_pipeline : contains the reads table that was created with the optimal quality score.

With our small sample, the AMA final q value was: 10

The AMA max read count was: 114561

And the final AMA trim and filter table looked like this:

merge_read_tracking should add a file {target}_track_reads_through_dada2.csv to the optimize reads directory. This file combines the trim and filter table from the trim_and_filter rule that corresponds to the optimal q value and the {target}_final_track_reads_through_pipeline table that was produced in the last step.

With our small sample, the AMA track_reads_through_dada2 table looked like this:

censor_haplotypes should produce a haplotypes folder for each target with two files {target}_haplotype_table_precensored.csv and {target}_haplotype_table_censored_final_version.csv . It should also produce a sequences folder with four files: {target}_snps_between_haps_within_samples.fasta , {target}_unique_seqs.fasta , {target}_aligned_seqs.fasta , and {target}_unique_seqs_final_censored.fasta .

• {target}_haplotype_table_precensored.csv : outputs the haplotype data set prior to beginning the censoring process (essentially {target}_haplotypes.rds in a formatted csv file).

• {target}_snps_between_haps_within_samples.fasta : fasta file of the haplotypes after the first three steps of the censoring process are completed. This is used to tally up the number of SNPs between all haplotype pairings.

• {target}_unique_seqs.fasta : fasta file of the haplotypes after the fourth step of the censoring process is completed.

• {target}_aligned_seqs.fasta : fasta file of the sequences after alignment.

• {target}_unique_seqs_final_censored.fasta : fasta file of the haplotype results after all five steps of the censoring process are completed.

• {target}_haplotype_table_censored_final_version.csv : outputs the final censored haplotype data set in a formatted table.

With our small sample, the AMA haplotype_table_precensored file looked like this:

The AMA snps_between_haps_within_samples file looked like this:

>Seq1
GTAAAGGTATAATTATTGAGAATTCAAATACTACTTTTTTAACACCGGTAGCTACGGGAAATCAATATTTAAAAGATGGAGGTTTTGCTTTTCCTCCAACAGAACCTCATATGTCACCAATGACATTAGATGAAATGAGACATTTTTATAAAGATAATAAATATGTAAAAAATTTAGATGAATTGACTTTATGTTCAAGACATGCAGGAAATATGATTCCAGATAATGATAAAAATTCAAATTATAAATATCCAGCTGTTTATGATGACAAAGATAAAAAGTGTCATATATTATATATTG
>Seq2
GTAAAGGTATAATTATTGAGAATTCAAATACTACTTTTTTAAAACCGGTAGCTACGGGAAATCAAGATTTAAAAGATGGAGGTTTTGCTTTTCCTCCAACAAATCCTCTTATATCACCAATGACATTAAATGGTATGAGAGATTTTTATAAAAATAATGAATATGTAAAAAATTTAGATGAATTGACTTTATGTTCAAGACATGCAGGAAATATGAATCCAGATAATGATGAAAATTCAAATTATAAATATCCAGCTGTTTATGATGACAAAGATAAAAAGTGTCATATATTATATATTG
>Seq3
GTAAAGGTATAATTATTGAGAATTCAAATACTACTTTTTTAAAACCGGTAGCTACGGGAAATCAAGATTTAAAAGATGGAGGTTTTGCTTTTCCTCCAACAGAACCTCTTATATCACCAATGACATTAGATGATATGAGAGATTTTTATAAAAATAATGAATATGTAAAAAATTTAGATGAATTGACTTTATGTTCAAGACATGCAGGAAATATGAATCCAGATAATGATCAAAATTCAAATTATAAATATCCAGCTGTTTATGATTACGAAGATAAAAAGTGTCATATATTATATATTG
>Seq4
GTAAAGGTATAATTATTGAGAATTCAAATACTACTTTTTTAACACCGGTAGCTACGGGAAAACAAGATTTAAAAGATGGAGGTTTTGCTTTTCCTCCAACAAATCCTCTTATATCACCAATGACATTAGATCATATGAGAGATTTTTATAAAAAAAATGAATATGTAAAAAATTTAGATGAATTGACTTTATGTTCAAGACATGCAGGAAATATGAATCCAGATAATGATGAAAATTCAAATTATAAATATCCAGCTGTTTATGATTACAAAGATAAAAAGTGTCATATATTATATATTG
>Seq5
GTAAAGGTATAATTATTGAGAATTCAAATACTACTTTTTTAAAACCGGTAGCTACGGGAAATCAAGATTTAAAAGATGGAGGTTTTGCTTTTCCTCCAACAGAACCTCTTATATCACCAATGACATTAAATGGTATGAGAGATTTATATAAAAATAATGAAGATGTAAAAAATTTAGATGAATTGACTTTATGTTCAAGACATGCAGGAAATATGAATCCAGATAATGATAAAAATTCAAATTATAAATATCCAGCTGTTTATGATGACAAAAATAAAAAGTGTCATATATTATATATTG
>Seq6
GTAAAGGTATAATTATTGAGAATTCAAAAACTACTTTTTTAACACCGGTAGCTACGGAAAATCAAGATTTAAAAGATGGAGGTTTTGCTTTTCCTCCAACAAATCCTCCTATGTCACCAATGACATTAAATGGTATGAGAGATTTATATAAAAATAATGAATATGTAAAAAATTTAGATGAATTGACTTTATGTTCAAGACATGCAGGAAATATGAATCCAGATAATGATAAAAATTCAAATTATAAATATCCAGCTGTTTATGATTACAATGATAATAAGTGTCATATATTATATATTG
>Seq7
GTAAAGGTATAATTATTGAGAATTCAAATACTACTTTTTTAAAACCGGTAGCTACGGGAAATCAAGATTTAAAAGATGGAGGTTTTGCTTTTCCTCCAACAGAACCTCTTATATCACCAATGACATTAAATGGTATGAGAGATTTTTATAAAAATAATGAATATGTAAAAAATTTAGATGAATTGACTTTATGTTCAAGACATGCAGGAAATATGAATCCAGATAAGGATGAAAATTCAAATTATAAATATCCAGCTGTTTATGATGACAAAGATAAAAAGTGTCATATATTATATATTG

The AMA unique_seqs file looked like this:

>Seq1
GTAAAGGTATAATTATTGAGAATTCAAATACTACTTTTTTAACACCGGTAGCTACGGGAAATCAATATTTAAAAGATGGAGGTTTTGCTTTTCCTCCAACAGAACCTCATATGTCACCAATGACATTAGATGAAATGAGACATTTTTATAAAGATAATAAATATGTAAAAAATTTAGATGAATTGACTTTATGTTCAAGACATGCAGGAAATATGATTCCAGATAATGATAAAAATTCAAATTATAAATATCCAGCTGTTTATGATGACAAAGATAAAAAGTGTCATATATTATATATTG
>Seq2
GTAAAGGTATAATTATTGAGAATTCAAATACTACTTTTTTAAAACCGGTAGCTACGGGAAATCAAGATTTAAAAGATGGAGGTTTTGCTTTTCCTCCAACAAATCCTCTTATATCACCAATGACATTAAATGGTATGAGAGATTTTTATAAAAATAATGAATATGTAAAAAATTTAGATGAATTGACTTTATGTTCAAGACATGCAGGAAATATGAATCCAGATAATGATGAAAATTCAAATTATAAATATCCAGCTGTTTATGATGACAAAGATAAAAAGTGTCATATATTATATATTG
>Seq3
GTAAAGGTATAATTATTGAGAATTCAAATACTACTTTTTTAAAACCGGTAGCTACGGGAAATCAAGATTTAAAAGATGGAGGTTTTGCTTTTCCTCCAACAGAACCTCTTATATCACCAATGACATTAGATGATATGAGAGATTTTTATAAAAATAATGAATATGTAAAAAATTTAGATGAATTGACTTTATGTTCAAGACATGCAGGAAATATGAATCCAGATAATGATCAAAATTCAAATTATAAATATCCAGCTGTTTATGATTACGAAGATAAAAAGTGTCATATATTATATATTG
>Seq4
GTAAAGGTATAATTATTGAGAATTCAAATACTACTTTTTTAACACCGGTAGCTACGGGAAAACAAGATTTAAAAGATGGAGGTTTTGCTTTTCCTCCAACAAATCCTCTTATATCACCAATGACATTAGATCATATGAGAGATTTTTATAAAAAAAATGAATATGTAAAAAATTTAGATGAATTGACTTTATGTTCAAGACATGCAGGAAATATGAATCCAGATAATGATGAAAATTCAAATTATAAATATCCAGCTGTTTATGATTACAAAGATAAAAAGTGTCATATATTATATATTG
>Seq5
GTAAAGGTATAATTATTGAGAATTCAAATACTACTTTTTTAAAACCGGTAGCTACGGGAAATCAAGATTTAAAAGATGGAGGTTTTGCTTTTCCTCCAACAGAACCTCTTATATCACCAATGACATTAAATGGTATGAGAGATTTATATAAAAATAATGAAGATGTAAAAAATTTAGATGAATTGACTTTATGTTCAAGACATGCAGGAAATATGAATCCAGATAATGATAAAAATTCAAATTATAAATATCCAGCTGTTTATGATGACAAAAATAAAAAGTGTCATATATTATATATTG
>Seq6
GTAAAGGTATAATTATTGAGAATTCAAAAACTACTTTTTTAACACCGGTAGCTACGGAAAATCAAGATTTAAAAGATGGAGGTTTTGCTTTTCCTCCAACAAATCCTCCTATGTCACCAATGACATTAAATGGTATGAGAGATTTATATAAAAATAATGAATATGTAAAAAATTTAGATGAATTGACTTTATGTTCAAGACATGCAGGAAATATGAATCCAGATAATGATAAAAATTCAAATTATAAATATCCAGCTGTTTATGATTACAATGATAATAAGTGTCATATATTATATATTG
>Seq7
GTAAAGGTATAATTATTGAGAATTCAAATACTACTTTTTTAAAACCGGTAGCTACGGGAAATCAAGATTTAAAAGATGGAGGTTTTGCTTTTCCTCCAACAGAACCTCTTATATCACCAATGACATTAAATGGTATGAGAGATTTTTATAAAAATAATGAATATGTAAAAAATTTAGATGAATTGACTTTATGTTCAAGACATGCAGGAAATATGAATCCAGATAAGGATGAAAATTCAAATTATAAATATCCAGCTGTTTATGATGACAAAGATAAAAAGTGTCATATATTATATATTG

The AMA aligned_seqs file looked like this:

>Seq5
GTAAAGGTATAATTATTGAGAATTCAAATACTACTTTTTTAAAACCGGTAGCTACGGGAAATCAAGATTTAAAAGATGGA
GGTTTTGCTTTTCCTCCAACAGAACCTCTTATATCACCAATGACATTAAATGGTATGAGAGATTTATATAAAAATAATGA
AGATGTAAAAAATTTAGATGAATTGACTTTATGTTCAAGACATGCAGGAAATATGAATCCAGATAATGATAAAAATTCAA
ATTATAAATATCCAGCTGTTTATGATGACAAAAATAAAAAGTGTCATATATTATATATTG
>Seq6
GTAAAGGTATAATTATTGAGAATTCAAAAACTACTTTTTTAACACCGGTAGCTACGGAAAATCAAGATTTAAAAGATGGA
GGTTTTGCTTTTCCTCCAACAAATCCTCCTATGTCACCAATGACATTAAATGGTATGAGAGATTTATATAAAAATAATGA
ATATGTAAAAAATTTAGATGAATTGACTTTATGTTCAAGACATGCAGGAAATATGAATCCAGATAATGATAAAAATTCAA
ATTATAAATATCCAGCTGTTTATGATTACAATGATAATAAGTGTCATATATTATATATTG
>Seq4
GTAAAGGTATAATTATTGAGAATTCAAATACTACTTTTTTAACACCGGTAGCTACGGGAAAACAAGATTTAAAAGATGGA
GGTTTTGCTTTTCCTCCAACAAATCCTCTTATATCACCAATGACATTAGATCATATGAGAGATTTTTATAAAAAAAATGA
ATATGTAAAAAATTTAGATGAATTGACTTTATGTTCAAGACATGCAGGAAATATGAATCCAGATAATGATGAAAATTCAA
ATTATAAATATCCAGCTGTTTATGATTACAAAGATAAAAAGTGTCATATATTATATATTG
>Seq1
GTAAAGGTATAATTATTGAGAATTCAAATACTACTTTTTTAACACCGGTAGCTACGGGAAATCAATATTTAAAAGATGGA
GGTTTTGCTTTTCCTCCAACAGAACCTCATATGTCACCAATGACATTAGATGAAATGAGACATTTTTATAAAGATAATAA
ATATGTAAAAAATTTAGATGAATTGACTTTATGTTCAAGACATGCAGGAAATATGATTCCAGATAATGATAAAAATTCAA
ATTATAAATATCCAGCTGTTTATGATGACAAAGATAAAAAGTGTCATATATTATATATTG
>Seq2
GTAAAGGTATAATTATTGAGAATTCAAATACTACTTTTTTAAAACCGGTAGCTACGGGAAATCAAGATTTAAAAGATGGA
GGTTTTGCTTTTCCTCCAACAAATCCTCTTATATCACCAATGACATTAAATGGTATGAGAGATTTTTATAAAAATAATGA
ATATGTAAAAAATTTAGATGAATTGACTTTATGTTCAAGACATGCAGGAAATATGAATCCAGATAATGATGAAAATTCAA
ATTATAAATATCCAGCTGTTTATGATGACAAAGATAAAAAGTGTCATATATTATATATTG
>Seq3
GTAAAGGTATAATTATTGAGAATTCAAATACTACTTTTTTAAAACCGGTAGCTACGGGAAATCAAGATTTAAAAGATGGA
GGTTTTGCTTTTCCTCCAACAGAACCTCTTATATCACCAATGACATTAGATGATATGAGAGATTTTTATAAAAATAATGA
ATATGTAAAAAATTTAGATGAATTGACTTTATGTTCAAGACATGCAGGAAATATGAATCCAGATAATGATCAAAATTCAA
ATTATAAATATCCAGCTGTTTATGATTACGAAGATAAAAAGTGTCATATATTATATATTG
>Seq7
GTAAAGGTATAATTATTGAGAATTCAAATACTACTTTTTTAAAACCGGTAGCTACGGGAAATCAAGATTTAAAAGATGGA
GGTTTTGCTTTTCCTCCAACAGAACCTCTTATATCACCAATGACATTAAATGGTATGAGAGATTTTTATAAAAATAATGA
ATATGTAAAAAATTTAGATGAATTGACTTTATGTTCAAGACATGCAGGAAATATGAATCCAGATAAGGATGAAAATTCAA
ATTATAAATATCCAGCTGTTTATGATGACAAAGATAAAAAGTGTCATATATTATATATTG

The AMA unique_seqs_final_censored file looked like this:

>Seq1
GTAAAGGTATAATTATTGAGAATTCAAATACTACTTTTTTAAAACCGGTAGCTACGGGAAATCAAGATTTAAAAGATGGAGGTTTTGCTTTTCCTCCAACAAATCCTCTTATATCACCAATGACATTAAATGGTATGAGAGATTTTTATAAAAATAATGAATATGTAAAAAATTTAGATGAATTGACTTTATGTTCAAGACATGCAGGAAATATGAATCCAGATAATGATGAAAATTCAAATTATAAATATCCAGCTGTTTATGATGACAAAGATAAAAAGTGTCATATATTATATATTG

And the AMA haplotype_table_censored_final table looked like this:

get_read_summaries should produce two files, pre-filt_fastq_read_counts and filt_fastq_read_counts , in the haplotype_output folder, and two files, trim_summaries and trim_summary_names , in the output/trimmed_reads folder.

• trim_summaries consolidates all the output/trimmed_reads/summaries data into one file.

• trim_summary_names contains all the file names in output/trimmed_reads/summaries .

• pre-filt_fastq_read_counts contains the read counts prior to filtering.

• filt-fastq_read_counts contains the read counts after filtering.

With our small sample, one of the consolidated trim summaries looked like this:

Input Read Pairs: 33577
Both Surviving Reads: 32298
Both Surviving Read Percent: 96.19
Forward Only Surviving Reads: 454
Forward Only Surviving Read Percent: 1.35
Reverse Only Surviving Reads: 659
Reverse Only Surviving Read Percent: 1.96
Dropped Reads: 166
Dropped Read Percent: 0.49

The trim summary names file looked like this:

BF10.summary
BF1.summary
BF246.summary
BF247.summary
BF248.summary
BF249.summary
BF250.summary
BF251.summary
BF252.summary
BF253.summary
BF254.summary
BF255.summary
BF256.summary
BF257.summary
BF258.summary
BF259.summary
BF260.summary
BF261.summary
BF262.summary
BF263.summary
BF2.summary
BF3.summary
BF4.summary
BF5.summary
BF6.summary
BF7.summary
BF8.summary
BF9.summary

Some of the pre-filtered read counts looked like this:

AMA/out/fastq/all_samples/BF10_AMA_1.fastq.gz	14875
AMA/out/fastq/all_samples/BF10_AMA_2.fastq.gz	14875
AMA/out/fastq/all_samples/BF1_AMA_1.fastq.gz	13072
AMA/out/fastq/all_samples/BF1_AMA_2.fastq.gz	13072
AMA/out/fastq/all_samples/BF2_AMA_1.fastq.gz	15269
AMA/out/fastq/all_samples/BF2_AMA_2.fastq.gz	15269
AMA/out/fastq/all_samples/BF3_AMA_1.fastq.gz	10116
AMA/out/fastq/all_samples/BF3_AMA_2.fastq.gz	10116

And some of the filtered read counts looked like this:

AMA/out/fastq/all_samples/final_filtered/BF10_final_F_filt.fastq.gz	13384
AMA/out/fastq/all_samples/final_filtered/BF10_final_R_filt.fastq.gz	13384
AMA/out/fastq/all_samples/final_filtered/BF1_final_F_filt.fastq.gz	11777
AMA/out/fastq/all_samples/final_filtered/BF1_final_R_filt.fastq.gz	11777
AMA/out/fastq/all_samples/final_filtered/BF2_final_F_filt.fastq.gz	13718
AMA/out/fastq/all_samples/final_filtered/BF2_final_R_filt.fastq.gz	13718
AMA/out/fastq/all_samples/final_filtered/BF3_final_F_filt.fastq.gz	8491
AMA/out/fastq/all_samples/final_filtered/BF3_final_R_filt.fastq.gz	8491

create_summaries should produce 2 files in the haplotype_output folder: long_summary and wide_summary . long_summary is a csv file that summarizes the read counts with columns read_type, read_ct, sample, and target. This allows for easy manipulation of the dataframe if needed in the future. wide_summary is a csv file that summarizes the read counts with each row being a sample. This allows for easy visualization of the data analysis to see how each sample did throughout the workflow.

With our small samples of AMA and CSP, the long summary looked like this:

And the wide summary looked like this:

More resources

• Snakemake tutorial

• Conda user guide

• MultiQC report tutorial

• MultiQC documentation