Skip to content

duke-malaria-collaboratory/BRAVA

Folders and files

NameName
Last commit message
Last commit date

Latest commit

 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 

Repository files navigation

BRAVA: Basic and Rigorous Amplicon Variant Analyzer

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

Input data: forward and reverse fastq files for each sample from targeted amplicon deep sequencing.

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.

dag

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/duke-malaria-collaboratory/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 config/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:

    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.

multiqc

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:

FIELD1 Target Sample POS REF ALT ALT_DEPTH DEPTH ALT_FREQ
0 PFCRT BF252 404 A . 0 11 0.0
1 PFCRT BF252 454 A . 0 11 0.0
2 PFCRT BF252 466 C . 0 12 0.0
3 PFCRT BF252 610 C . 0 7 0.0

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:

Target Sample POS BASE TOTAL_DEPTH DEPTH FREQ
PFMDR1 BF246 57 T 111 34 0.3063063063063063
PFMDR1 BF259 57 T 115 64 0.5565217391304348
PFMDR1 BF260 57 T 89 64 0.7191011235955056
PFMDR1 BF261 57 T 139 67 0.48201438848920863
PFMDR1 BF262 57 T 125 96 0.768
PFMDR1 BF263 57 T 137 127 0.927007299270073
PFMDR1 BF246 352 T 116 116 1
PFMDR1 BF247 352 T 68 68 1
PFMDR1 BF248 352 T 119 119 1
PFMDR1 BF249 352 T 88 88 1
PFMDR1 BF251 352 T 99 99 1
PFMDR1 BF252 352 T 84 67 0.7976190476190477
PFMDR1 BF254 352 T 80 80 1
PFMDR1 BF255 352 T 108 108 1
PFMDR1 BF257 352 T 79 78 0.9873417721518988
PFMDR1 BF259 352 T 118 49 0.4152542372881356
PFMDR1 BF261 352 T 143 43 0.3006993006993007
PFMDR1 BF246 57 A 111 77 0.6936936936936937
PFMDR1 BF259 57 A 115 51 0.4434782608695652
PFMDR1 BF260 57 A 89 25 0.2808988764044944
PFMDR1 BF261 57 A 139 72 0.5179856115107914
PFMDR1 BF262 57 A 125 29 0.23199999999999998
PFMDR1 BF263 57 A 137 10 0.07299270072992703
PFMDR1 BF246 352 A 116 0 0
PFMDR1 BF247 352 A 68 0 0
PFMDR1 BF248 352 A 119 0 0
PFMDR1 BF249 352 A 88 0 0
PFMDR1 BF251 352 A 99 0 0
PFMDR1 BF252 352 A 84 17 0.20238095238095233
PFMDR1 BF254 352 A 80 0 0
PFMDR1 BF255 352 A 108 0 0
PFMDR1 BF257 352 A 79 1 0.012658227848101222
PFMDR1 BF259 352 A 118 69 0.5847457627118644
PFMDR1 BF261 352 A 143 100 0.6993006993006994
PFMDR1 BF247 57 A 67 67 1
PFMDR1 BF248 57 A 118 118 1
PFMDR1 BF249 57 A 88 88 1
PFMDR1 BF250 57 A 81 81 1
PFMDR1 BF251 57 A 95 95 1
PFMDR1 BF252 57 A 85 85 1
PFMDR1 BF253 57 A 103 103 1
PFMDR1 BF254 57 A 79 79 1
PFMDR1 BF255 57 A 105 105 1
PFMDR1 BF256 57 A 121 121 1
PFMDR1 BF257 57 A 79 79 1
PFMDR1 BF258 57 A 94 94 1
PFMDR1 BF250 352 A 82 82 1
PFMDR1 BF253 352 A 102 102 1
PFMDR1 BF256 352 A 121 121 1
PFMDR1 BF258 352 A 95 95 1
PFMDR1 BF260 352 A 96 96 1
PFMDR1 BF262 352 A 126 126 1
PFMDR1 BF263 352 A 142 142 1
PFCRT BF252 404 A 11 11 1
PFCRT BF252 454 A 11 11 1
PFCRT BF252 466 C 12 12 1
PFCRT BF252 610 C 7 7 1

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:

reads.in reads.out
BF1_AMA_1.fastq.gz 13072 8297
BF10_AMA_1.fastq.gz 14875 9243
BF2_AMA_1.fastq.gz 15269 10247
BF3_AMA_1.fastq.gz 10116 5075
BF4_AMA_1.fastq.gz 13672 8852
BF5_AMA_1.fastq.gz 12650 7720
BF6_AMA_1.fastq.gz 14591 9754
BF7_AMA_1.fastq.gz 12760 8208
BF8_AMA_1.fastq.gz 11474 7382
BF9_AMA_1.fastq.gz 12432 6318

And the AMA read count table looked like this:

25 76565
15 103295
20 103224
2 81096
5 81096
10 116036

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:

merged tabled nonchim
BF1 8296 8296 8296
BF10 9241 9241 9241
BF2 10239 10239 10239
BF3 5056 5056 5056
BF4 8850 8850 8850
BF5 7614 7614 7210
BF6 9752 9752 9752
BF7 8125 8125 7795
BF8 7378 7378 7378
BF9 6317 6317 6317

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:

merged tabled nonchim
BF1 11768 11768 11768
BF10 13377 13377 13377
BF2 13681 13681 13681
BF3 8450 8450 8450
BF4 12211 12211 12211
BF5 10584 10584 10360
BF6 13078 13078 13078
BF7 11049 11049 10773
BF8 10314 10314 10314
BF9 10549 10549 10549

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:

sample reads.in reads.out merged tabled nonchim
1 BF1 13072 11777 11768 11768 11768
2 BF10 14875 13384 13377 13377 13377
3 BF2 15269 13718 13681 13681 13681
4 BF246 0 0 NA NA NA
5 BF247 0 0 NA NA NA
6 BF248 0 0 NA NA NA
7 BF249 0 0 NA NA NA
8 BF250 0 0 NA NA NA
9 BF251 0 0 NA NA NA
10 BF252 0 0 NA NA NA
11 BF253 0 0 NA NA NA
12 BF254 0 0 NA NA NA
13 BF255 0 0 NA NA NA
14 BF256 0 0 NA NA NA
15 BF257 0 0 NA NA NA
16 BF258 0 0 NA NA NA
17 BF259 0 0 NA NA NA
18 BF260 0 0 NA NA NA
19 BF261 0 0 NA NA NA
20 BF262 0 0 NA NA NA
21 BF263 0 0 NA NA NA
22 BF3 10116 8491 8450 8450 8450
23 BF4 13672 12226 12211 12211 12211
24 BF5 12650 11119 10584 10584 10360
25 BF6 14591 13085 13078 13078 13078
26 BF7 12760 11350 11049 11049 10773
27 BF8 11474 10326 10314 10314 10314
28 BF9 12432 10560 10549 10549 10549

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:

H1 H2 H3 H4 H5 H6 H7 MiSeq.ID
0 0 0 0 11768 0 0 BF1
0 0 0 13377 0 0 0 BF10
13665 0 0 16 0 0 0 BF2
0 8450 0 0 0 0 0 BF3
18 12186 7 0 0 0 0 BF4
268 0 9862 0 0 230 0 BF5
13078 0 0 0 0 0 0 BF6
0 603 10170 0 0 0 0 BF7
0 0 0 0 0 10314 0 BF8
51 0 0 0 0 0 10498 BF9

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:

H2 MiSeq.ID
8450 BF3
12186 BF4
603 BF7

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:

read_type read_ct sample target
Input Read Pairs 33577 BF10 NA
Both Surviving Reads 32298 BF10 NA
Forward Only Surviving Reads 454 BF10 NA
Reverse Only Surviving Reads 659 BF10 NA
Dropped Reads 166 BF10 NA
Input Read Pairs 26485 BF1 NA
Both Surviving Reads 25565 BF1 NA
Forward Only Surviving Reads 359 BF1 NA
Reverse Only Surviving Reads 461 BF1 NA
Dropped Reads 100 BF1 NA
Input Read Pairs 28008 BF2 NA
Both Surviving Reads 26932 BF2 NA
Forward Only Surviving Reads 412 BF2 NA
Reverse Only Surviving Reads 516 BF2 NA
Dropped Reads 148 BF2 NA
Input Read Pairs 25840 BF3 NA
Both Surviving Reads 24730 BF3 NA
Forward Only Surviving Reads 443 BF3 NA
Reverse Only Surviving Reads 504 BF3 NA
Dropped Reads 163 BF3 NA
Input Read Pairs 29520 BF4 NA
Both Surviving Reads 28500 BF4 NA
Forward Only Surviving Reads 392 BF4 NA
Reverse Only Surviving Reads 546 BF4 NA
Dropped Reads 82 BF4 NA
Input Read Pairs 30472 BF5 NA
Both Surviving Reads 29347 BF5 NA
Forward Only Surviving Reads 455 BF5 NA
Reverse Only Surviving Reads 533 BF5 NA
Dropped Reads 137 BF5 NA
Input Read Pairs 32025 BF6 NA
Both Surviving Reads 30863 BF6 NA
Forward Only Surviving Reads 444 BF6 NA
Reverse Only Surviving Reads 591 BF6 NA
Dropped Reads 127 BF6 NA
Input Read Pairs 29455 BF7 NA
Both Surviving Reads 28429 BF7 NA
Forward Only Surviving Reads 391 BF7 NA
Reverse Only Surviving Reads 524 BF7 NA
Dropped Reads 111 BF7 NA
Input Read Pairs 27293 BF8 NA
Both Surviving Reads 26279 BF8 NA
Forward Only Surviving Reads 303 BF8 NA
Reverse Only Surviving Reads 609 BF8 NA
Dropped Reads 102 BF8 NA
Input Read Pairs 26675 BF9 NA
Both Surviving Reads 25541 BF9 NA
Forward Only Surviving Reads 477 BF9 NA
Reverse Only Surviving Reads 514 BF9 NA
Dropped Reads 143 BF9 NA
filt 11777 BF1 ama
prefilt 13072 BF1 ama
filt 11352 BF1 csp
prefilt 12491 BF1 csp
filt 13384 BF10 ama
prefilt 14875 BF10 ama
filt 15886 BF10 csp
prefilt 17422 BF10 csp
filt 13718 BF2 ama
prefilt 15269 BF2 ama
filt 10659 BF2 csp
prefilt 11659 BF2 csp
filt 8491 BF3 ama
prefilt 10116 BF3 ama
filt 13099 BF3 csp
prefilt 14611 BF3 csp
filt 12226 BF4 ama
prefilt 13672 BF4 ama
filt 13584 BF4 csp
prefilt 14828 BF4 csp
filt 11119 BF5 ama
prefilt 12650 BF5 ama
filt 15006 BF5 csp
prefilt 16695 BF5 csp
filt 13085 BF6 ama
prefilt 14591 BF6 ama
filt 14725 BF6 csp
prefilt 16269 BF6 csp
filt 11350 BF7 ama
prefilt 12760 BF7 ama
filt 14179 BF7 csp
prefilt 15668 BF7 csp
filt 10326 BF8 ama
prefilt 11474 BF8 ama
filt 13359 BF8 csp
prefilt 14801 BF8 csp
filt 10560 BF9 ama
prefilt 12432 BF9 ama
filt 11500 BF9 csp
prefilt 13108 BF9 csp

And the wide summary looked like this:

sample NA_Input Read Pairs NA_Both Surviving Reads NA_Forward Only Surviving Reads NA_Reverse Only Surviving Reads NA_Dropped Reads ama_filt ama_prefilt csp_filt csp_prefilt
BF10 33577 32298 454 659 166 13384 14875 15886 17422
BF1 26485 25565 359 461 100 11777 13072 11352 12491
BF2 28008 26932 412 516 148 13718 15269 10659 11659
BF3 25840 24730 443 504 163 8491 10116 13099 14611
BF4 29520 28500 392 546 82 12226 13672 13584 14828
BF5 30472 29347 455 533 137 11119 12650 15006 16695
BF6 32025 30863 444 591 127 13085 14591 14725 16269
BF7 29455 28429 391 524 111 11350 12760 14179 15668
BF8 27293 26279 303 609 102 10326 11474 13359 14801
BF9 26675 25541 477 514 143 10560 12432 11500 13108

More resources

About

No description, website, or topics provided.

Resources

License

Stars

Watchers

Forks

Releases

No releases published

Packages

No packages published

Languages