--- title: "An Introduction to microhaplot data preparation" author: "Thomas Ng and Eric C. Anderson" date: "`r Sys.Date()`" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{microhaplot data prep} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- This vignette concerns itself with an overview of microhaplot and on preparation of data for input into microhaplot Once your data has been transferred to the appropriate Haplot folder, see the `microhaplot walkthrough` vignette to learn more about microhaplot's graphical features. ## How does microhaplot work? `microhaplot` is a tool that makes it easy to extract, visualize, and interpret short-read haplotypes---or "microhaplotypes" in the parlance of some---from aligned sequences. Before we go further, it will be good to step back and confirm of what we are calling as a short-read haplotype, or a "microhaplotype" within the context of microhaplot. > *What is a microhaplotype?* To define a microhaplotype you need to first define a reference point or position that a read might align to. Then, `microhaplot` can analyze the reads aligning at those positions, and extracts the variant sequence of the reads at those positions. For example, let's say that we know SNPs occur at positions 19, 27, 68, and 101 of a particular reference sequence that we have aligned our reads to. If a read aligns to the ref and shows an `A` at position 19, a `G` at position 27, a `C` at position 68, and another `C` at position 101, then we say the microhaplotype on that read is `AGCC`. `microhaplot` is a tool that 1. automates the extraction of these haplotypes from aligned data (typically amplicon sequencing data) 2. provides useful filtering and visualizing capabilities for understanding your microhaplotype data and for calling genotypes at microhaplotypes. To understand what is involved in getting your data into `microhaplot`, it will be useful to have a brief introduction to what it does and how it operates. We will start by listing what it is **not**. - `microhaplot` **is not a SNP-caller or variant detector**. It is not its job to determine if a variant might exist in a particular location or not. We assume you already know what variant position you want to look for by supplying `microhaplot` a variant caller file (VCF). `microhaplot` just extracts the bases sequenced on different reads at those positions. - `microhaplot` **is not an alignment program**. `microhaplot` does not do alignments. You have to provide aligned data to it. From the aligned data, `microhaplot` can extract the bases sequenced at the desired locations. - `microhaplot` **is not a haplotype phasing program**. There is no special algorithm going on in `microhaplot` to try to assemble haplotypes (across reads or loci, for example.) It merely uses the simple fact that a single sequencing read will come from a single chromosome. So, if you can draw bases present at 4 variant positions along a 150-bp read, you know that those 4 bases constitute a haplotype. - `microhaplot` **is not a de-novo assembly algorithm**. This package will not infer or sort out any physical linkage between reads of different reference loci. ## What kind of data can I use with `microhaplot`? 'microhaplot' is designed to work with short read sequencing data where there is a relatively small number sequencing start sites. It is emphatically *not* geared to whole-genome shotgun sequencing data in which you may have overlapping reads with different starting place. We use 'microhaplot' primarily with amplicon sequencing data in which a few hundred fragments (of about 100 - 150 bp each) are amplified by PCR in a few hundred individuals. The DNA from different individuals is bar-coded so that it may be demultiplexed later, and then these amplified fragments are sequenced (usually to fairly high read depth) on an Illumina platform. Though we use 'microhaplot' for amplicon sequencing, it might also be useful for analyzing haplotype data from experiments that sequence captured or baited fragments, (e.g., baited DD-Rad, RAPTURE, myBaits, ExonCapture, etc.), and might even be applicable to standard ddRAD or RAD-PE experiments, however, as the number of loci increases and average read depths per individual-locus combination decreases, the unevenness of coverage can become a bigger problem, (and we have more reservations about allelic dropout/null alleles when using digestion-based protocols (ddRAD, Rapture, etc.) than we do with amplicon sequencing). ## Why did we write 'microhaplot'? We tried to use existing tools to give us what we needed for our amplicon-based microhaplotypes. - We were hopeful that [STACKS](http://catchenlab.life.illinois.edu/stacks/) would provide us useful microhaplotypes. What we found was that the presence of indels caused STACKS to fail, and even if it didn't we never found a reliable way to export inferred read haplotypes. (We could see what looked like haplotypes in the SQL-based viewer, but couldn't figure out how to get at those in an exported fashion!) - We spent countless months trying to get [freebayes](https://github.com/ekg/freebayes) to properly call microhaplotypes. Despite being an excellent variant detector, however, when presented with situations that were clearly simple microhaplotypes, `freebayes` would fail to call them, or would return results in a very difficult to manage "multinucleotide polymorphism" in the VCF. - We thought that `GATK's Unified Haplotype Caller` might deliver us from the bioinformatic hell we were experiencing, but alas, with the super-high read depths at some of the fragments in our amplicon sequencing data, GATK was brought to its knees and didn't offer a viable solution. It appears that the HaplotypeCaller is trying to solve a much harder problem than just extracting short-read haplotypes. After a few months of this, we realized that there is no software tailored for the relatively straightforward task of extracting and analyzing microhaplotypes. The other software programs have all been designed for different purposes, and, as a consequence, while they work great for identifying variants or performing de-novo assemblies, one might not expect them to work for the specialized task of extracting and analyzing microhaplotypes. ## How, specifically, should I prepare my data? The starting point for `microhaplot` requires: 1. An indexed, SAM file of reads aligned to the reference sequences for each amplicon. 2. A VCF file whose contents define the positions in each reference sequence from which you want to extract variants into microhaplotypes. Note that the *individuals* that appear in the VCF file need not be the same ones whose reads appear in the assemblies. The VCF is merely used as a way of defining positions to extract. Still, since the content of microhaplotype is based on the variant positions listed in the VCF file, it is crucial for the VCF file to only contain highly reliable variant sites. We describe here the workflow that we use to get these two necessary files from the short-read amplicon sequencing data that come off our Illumina Mi-seq. This workflow can be tweaked and tailored for your own data. Each step is described in one of the following subsections. ### Starting state Out description starts at the point where we have copied all of our fastq.gz files from the sequencer into a directory called `rawdata`. Next to that directory we have made two more directories, `flash` and `map`, in which we will be creating files and doing work. We have named all the Illumina files so that a directory listing of `rawdata` looks like this: ```{sh dir-list, eval = FALSE} /rawdata/--% ll | head total 3517864 drwxrwxr-x 2 biopipe biopipe 122880 Aug 18 11:57 ./ drwxrwxr-x 5 biopipe biopipe 4096 Aug 18 11:31 ../ -rw-r--r-- 1 biopipe biopipe 672293 Aug 18 11:53 satro_S100_L001_I1_001.fastq.gz -rw-r--r-- 1 biopipe biopipe 725385 Aug 18 11:53 satro_S100_L001_I2_001.fastq.gz -rw-r--r-- 1 biopipe biopipe 4656145 Aug 18 11:53 satro_S100_L001_R1_001.fastq.gz -rw-r--r-- 1 biopipe biopipe 4749637 Aug 18 11:53 satro_S100_L001_R2_001.fastq.gz -rw-r--r-- 1 biopipe biopipe 635863 Aug 18 11:53 satro_S101_L001_I1_001.fastq.gz -rw-r--r-- 1 biopipe biopipe 711534 Aug 18 11:53 satro_S101_L001_I2_001.fastq.gz -rw-r--r-- 1 biopipe biopipe 4410826 Aug 18 11:53 satro_S101_L001_R1_001.fastq.gz ``` In this case, there are 384 such sets of files. The number after the capital S in the filename gives the ID number of each individual fish. ### Flashing the reads together Our data are paired end reads of fragments that are short enough that the primary reads and the paired-end reads overlap. Rather than treat them as two separate reads, we combine them together into a single read using [flash](https://ccb.jhu.edu/software/FLASH/). Here is the command we use to cycle over all the files and flash them, run from the `flash` directory which is at the same level as the `rawdata` directory. (Note that the `/flash/--%` in the following code fragment is just the command prompt---it tells you which directory we are in...) ```{sh flash, eval = FALSE} /flash/--% for i in {1..384}; do flash -m 10 -M 100 -z --output-prefix=S${i} ../rawdata/satro_S${i}_L001_R1_001.fastq.gz ../rawdata/satro_S${i}_L001_R2_001.fastq.gz; done ``` The output that `flash` spits out to stdout looks like this: ```{sh flash-output, eval = FALSE} [FLASH] Starting FLASH v1.2.11 [FLASH] Fast Length Adjustment of SHort reads [FLASH] [FLASH] Input files: [FLASH] ../rawdata/satro_S1_L001_R1_001.fastq.gz [FLASH] ../rawdata/satro_S1_L001_R2_001.fastq.gz [FLASH] [FLASH] Output files: [FLASH] ./S1.extendedFrags.fastq.gz [FLASH] ./S1.notCombined_1.fastq.gz [FLASH] ./S1.notCombined_2.fastq.gz [FLASH] ./S1.hist [FLASH] ./S1.histogram [FLASH] [FLASH] Parameters: [FLASH] Min overlap: 10 [FLASH] Max overlap: 100 [FLASH] Max mismatch density: 0.250000 [FLASH] Allow "outie" pairs: false [FLASH] Cap mismatch quals: false [FLASH] Combiner threads: 12 [FLASH] Input format: FASTQ, phred_offset=33 [FLASH] Output format: FASTQ, phred_offset=33, gzip [FLASH] ``` This has created a series of files named like `S${i}.extendedFrags.fastq.gz` in which the `${i}` is replaced by a number between 1 and 384. These files are gzipped fastq files that hold the flashed reads. ### Mapping/Aligning the Reads The reference sequences we will use to align to are in a fasta file called `gtseq15_loci.fasta`. We put that in the `map` directory and then commence mapping the flashed reads as follows: ```{sh mapping, eval=FALSE} # index the reference. Output will be prefixed with "gtseq15_loci" /map/--% bwa index -p gtseq15_loci -a is gtseq15_loci.fasta # here is what comes by on stdout [bwa_index] Pack FASTA... 0.00 sec [bwa_index] Construct BWT for the packed sequence... [bwa_index] 0.00 seconds elapse. [bwa_index] Update BWT... 0.00 sec [bwa_index] Pack forward-only FASTA... 0.00 sec [bwa_index] Construct SA from BWT and Occ... 0.00 sec [main] Version: 0.7.5a-r405 [main] CMD: bwa index -p gtseq15_loci -a is gtseq15_loci.fasta [main] Real time: 0.061 sec; CPU: 0.007 sec # map reads. for each input file, the output is named satro_flashed_s${i}_aln.sam # Note the -R to insert appropriate read group header lines /map/--% for i in {1..384}; do bwa mem -aM -v 3 -t 10 -R "@RG\tID:s${i}\tLB:amplicon\tPL:ILLUMINA\tSM:satro${i}" ./gtseq15_loci ../flash/S${i}.extendedFrags.fastq.gz > ./satro_flashed_s${i}_aln.sam; done # convert the SAMs to BAMs /map/--% for i in {1..384}; do samtools view -bhS satro_flashed_s${i}_aln.sam > satro_flashed_s${i}.bam; done # Sort the BAMs /map/--% for i in {1..384}; do samtools sort satro_flashed_s${i}.bam -o satro_flashed_s${i}_sorted.bam; done # Index the BAMs /map/--% for i in {1..384}; do samtools index satro_flashed_s${i}_sorted.bam; done # Make a text-file list of the BAMs /map/--% ls -1 *sorted.bam > bamlist.txt # generate idx stats and put them in a directory called gtseq17_idx /map/--% mkdir gtseq17_idx /map/--% for i in {1..384}; do samtools idxstats satro_flashed_s${i}_sorted.bam > ./gtseq17_idx/s${i}_idxstats.txt; done ``` ### Variant detection We are now set up to detect variants using freebayes. We show how we do that here, but emphasize that if you have already chosen a set of variants that form your microhaplotypes, you could just specify those with a pre-existing VCF file. You *don't* need to do new variant detection if you don't want to. (Of course, if it doesn't take too long, you may as well do variant detection with every individual you have sequenced at the amplicons, so as to get as many variants as possible.) ```{sh freebayes, eval=FALSE} # index fasta reference file /map/--% samtools faidx gtseq15_loci.fasta # call variants, instructing freebayes to NOT return MNPs or other complex variants nohup freebayes-parallel <(fasta_generate_regions.py gtseq15_loci.fasta.fai 150) 10 -f gtseq15_loci.fasta -L bamlist.txt --haplotype-length 0 -kwVa --no-mnps --no-complex > satro384_noMNP_noComplex_noPriors.vcf ``` ### runHaplot After the above is done, the file `satro384_noMNP_noComplex_noPriors.vcf` can be filtered, if desired, to make sure that everything in it has a solid SNP. Then that file is used with the function `prepHaplotFiles` to extract haplotypes from the aligned reads in the `map` directory. An example of an vcf file and SAM files extracted from an actual GT-seq rockfish data is available in the `inst/extdata`. ```{r runHaplot, eval=FALSE} library(microhaplot) run.label <- "sebastes" sam.path <- tempdir() untar(system.file("extdata", "sebastes_sam.tar.gz", package="microhaplot"), exdir = sam.path) label.path <- file.path(sam.path, "label.txt") vcf.path <- file.path(sam.path, "sebastes.vcf") mvShinyHaplot(tempdir()) app.path <- file.path(tempdir(), "microhaplot") haplo.read.tbl <- prepHaplotFiles(run.label = run.label, sam.path = sam.path, out.path = tempdir(), label.path = label.path, vcf.path = vcf.path, app.path = app.path) ```