---
title: "Variant and Genotype Calling in Highly Duplicated Genomes"
author: "Lindsay V. Clark, University of Illinois, Urbana-Champaign"
date: "`r format(Sys.time(), '%d %B %Y')`"
output:
  rmarkdown::html_vignette:
    toc: true
vignette: >
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteIndexEntry{Variant and Genotype Calling in Highly Duplicated Genomes}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## When and why to use this pipeline

The variant and genotype calling pipeline described in this vignette is
intended for recent or ancient allopolyploid species, in which the reference
genome sequence contains many paralogous regions that do not recombine with
each other at meiosis.  In highly duplicated genomes such as these, conventional
alignment software, especially when restricted to only returning one alignment,
will often align sequences to the incorrect paralog, causing issues with variant
and genotype calling downstream.  The pipeline described in this vignette uses
read depth distribution across a population (the $H_{ind}/H_E$ statistic) to
determine whether a group of tags is behaving like a Mendelian locus. Attempting
to optimize this statistic, it rearranges tags among paralogous alignment locations
to correct alignments and generate well-behaved markers.

This pipeline is currently only designed for natural populations, diversity
panels, or mapping populations in which the most recent generation was produced
by random intermating among all progeny.  F1, self-fertilized, and backcrossed
populations will not be processed correctly because expected heterozygosity
cannot be estimated directly from allele frequencies in these populations.
However, the $H_{ind}/H_E$ statistic can still be used from within polyRAD
to filter non-Mendelian markers from these populations using the `HindHeMapping`
function.

If you use this pipeline, please cite:

Clark LV, Mays W, Lipka AE, and Sacks EJ (2020) A population-level statistic
for assessing Mendelian behavior of genotyping-by-sequencing data from highly
duplicated genomes.  *bioRxiv*
doi: [10.1101/2020.01.11.902890](https://doi.org/10.1101/2020.01.11.902890)

## Pipeline overview

Below is an overview of the pipeline.  It uses a combination of existing
bioinformatics tools, custom Python scripts included with polyRAD, and R
functions within polyRAD.

1. Use TASSEL-GBSv2 to find all unique sequence tags in a dataset, and the
depths of those tags in all taxa.
2. Align tags to the reference genome using Bowtie2, allowing multiple
alignments to be returned.
3. Use `process_sam_multi.py` to group tags based on sets of alignment
locations.
4. Import a small subset of the resulting data to polyRAD using
`readProcessSamMulti`.  Run `HindHe` on the imported data to identify
problematic samples and estimate inbreeding.
5. Use `process_isoloci.py` to sort tags into their correct alignment locations
within the groups of paralogs identified by `process_sam_multi.py`, and filter
out tags that cannot be adequately sorted.
6. Import the output of the previous step into polyRAD using
`readProcessIsoloci`, then perform genotype calling.

Everything except for the R portion of this pipeline will need to be run on your
operating system's Terminal/Shell/Command Prompt.  Bowtie2 requires Mac or
Linux.  Python 3 is needed for the Python portions of the pipeline.  If any of
this presents a barrier to you, I recommend finding someone experienced with
bioinformatics at your institution who can give you some advice.

Line breaks have been inserted in some terminal commands for readability in this
tutorial.  (Those same line breaks should not be typed at the command line.)

## TASSEL-GBSv2

An overview of the TASSEL-GBSv2 pipeline is available at
https://bitbucket.org/tasseladmin/tassel-5-source/wiki/Tassel5GBSv2Pipeline.
TASSEL can be downloaded from https://www.maizegenetics.net/tassel.

We will only run the first couple steps of the TASSEL-GBS pipeline.  The files
that you will need to start with are the
[key file](https://bitbucket.org/tasseladmin/tassel-5-source/wiki/Tassel5GBSv2Pipeline/KeyFileExample)
and your FASTQ files.

1. Run the [GBSSeqToTagDBPlugin](https://bitbucket.org/tasseladmin/tassel-5-source/wiki/Tassel5GBSv2Pipeline/GBSSeqToTagDBPlugin).
This requires your key file and FASTQ files as input, and will generate a
database file.
2. Run the [TagExportToFastqPlugin](https://bitbucket.org/tasseladmin/tassel-5-source/wiki/Tassel5GBSv2Pipeline/TagExportToFastqPlugin).
This requires your database file as input, and outputs a FASTQ with each tag
present once.  Set the `-c` parameter to something higher than the default
(10 is still pretty conservative but will discard most sequencing errors) to
save some processing downstream.
3. Run the [GetTagTaxaDistFromDBPlugin](https://bitbucket.org/tasseladmin/tassel-5-source/wiki/Tassel5GBSv2Pipeline/GetTagTaxaDistFromDBPlugin).  This requires the database as input, and outputs a tab-delimited
text file showing the depth of each tag in each sample.

The commands might look like (on Windows):

```
run_pipeline.bat -fork1 -GBSSeqToTagDBPlugin -e PstI-MspI -i D:\Msa\raw_data\
-db D:\Msa\Msa.db -k D:\Msa\key.txt -kmerLength 80 -endPlugin -runfork1

run_pipeline.bat -fork1 -TagExportToFastqPlugin -c 10 -db D:\Msa\Msa.db
-o D:\Msa\Msa_tags.fq -endPlugin -runfork1

run_pipeline.bat -fork1 -GetTagTaxaDistFromDBPlugin -db D:\Msa\Msa.db
-o D:\Msa\Msa_ttd.txt -endPlugin -runfork1
```

You may need to adjust the memory parameters to accomodate a large dataset.

### If TASSEL-GBS won't work for your dataset

If for some reason you can't use TASSEL-GBS (for example, your FASTQ files don't
have inline barcodes, or your protocol doesn't involve restriction enzymes) you
might be able to craft a custom alternative for yourself. What you will need out
of this step is a file formatted identically to a TagTaxaDist file from TASSEL,
and a FASTA or FASTQ file containing the same sequence tags as the
TagTaxaDist file.

A TagTaxaDist file is a tab-delimited text file formatted as below, and contains
the read depth of every tag in every sample (taxon).  I have used a simple
convention for sample names here, but the names can be anything you want.
Depending on the complexity of your dataset, it can be quite a large file.

| Tag | Sam1 | Sam2 | Sam3| Sam4 |
|-----|------|------|-----|------|
| TGCAGAAATCATAGATTAAGGATAT| 0 | 24 | 1 | 0 |
| TGCAGAACAGGATACGATACCCCTT| 2 | 0 | 0 | 0 |
| TGCAGAACAGTATACGATACCCCTT| 0 | 5 | 0 | 0 |

## Alignment with Bowtie2

Next, align the FASTQ file generated in step 2 of the previous section to your
reference genome.  If you have never used
[Bowtie2](http://bowtie-bio.sourceforge.net/bowtie2/index.shtml)
on your reference genome before, you will first need to
[make an index](http://bowtie-bio.sourceforge.net/bowtie2/manual.shtml#the-bowtie2-build-indexer).

```
bowtie2-build Msinensis_497_v7.0.hardmasked.fa Msi_DH1_7
```

Importantly, you will need to set the `-k` argument to something higher than the
number of subgenomes within your reference genome.  For example, in an
allotetraploid do at least `-k 3`, and in an allohexaploid at least `-k 4`.

For my *Miscanthus* dataset, I ran:

```
bowtie2 -k 3 --very-sensitive -x Msi_DH1_7 -U Msa_tags.fq -S Msa_align.sam
```

## Grouping tags by alignment sets

Next you will run `process_sam_multi.py`.  To locate this file, in R run

```{r}
system.file("python", "process_sam_multi.py", package = "polyRAD")
```

For convenience, you might copy this file to somewhere closer to your working
directory where you have the files generated in the previous step.  While you
are at it, also copy `isoloci_fun.py` and `process_isoloci.py` to somewhere
convenient.

As input, you will need the SAM file generated by Bowtie2 and the TagTaxaDist
file from TASSEL.  Be sure to set the `-g` argument to the number of subgenomes
if it is different from the default of 2 (allotetraploids); for example, in
an allohexaploid you should set `-g 3`.

For descriptions of other arguments, see

```
python process_sam_multi.py --help
```

The `-c` argument is useful if you want to divide your data into smaller chunks
for downstream processing (if you are limited on RAM or have many processors at
your disposal).  If there were some samples in your TASSEL database that you
already know you want to get rid of, you can list the samples that you want to
keep in a text file (one sample name per line) and pass that to the `-s`
argument.

In my *Miscanthus* data, the command looked like:

```
python process_sam_multi.py Msa_align.sam Msa_ttd.txt split_depths/Msa_split -c 5
```

The output will be a pair of CSVs for each chunk. One is named "align",
containing the set of alignment locations for each tag, with their respective
number of mutations and CIGAR strings.  The other, named "depth", is a subsetted
version of the TagTaxaDist file, containing the same tags in the same order as
"align".

## Filtering samples and estimating inbreeding

Before proceeding, we will want to get a look at the data to identify any
outlier samples that may be a different ploidy from the rest, interspecific
hybrids, or highly contaminated.  I recommend removing them to avoid biasing
$H_{ind}/H_E$ estimates, but keep in mind that this means they will be excluded
from all downstream analysis.  If a sample is important to your study and
you know of a reason why it would have a different heterozygosity from most
samples, you should keep it.

We will import just 1000 loci using `readProcessSamMulti`.  Sequence tags will
be assigned preliminary alignment locations based on where they had highest
sequence similarity, or to a random location if there was a tie.  This is
essentially a preview of what the data would look like if we went ahead with
variant calling using the top alignments returned by Bowtie2.

```{r}
library(polyRAD)
```
```{r eval = FALSE}
myRADprelim <- readProcessSamMulti("Msa_split_1_align.csv")
```

We will then estimate a $H_{ind}/H_E$ matrix for this data, as well as getting
a sum of the read depth for each individual.

```{r eval = FALSE, echo = FALSE}
# subset the object to have diploids and a few tetras
diploids <- readLines("diploids.txt")
myRADprelim <- SubsetByTaxon(myRADprelim, c(diploids, "KMS397", "KMS444", "UI11-00032"))
```
```{r eval = FALSE}
hh <- HindHe(myRADprelim)
TotDepthT <- rowSums(myRADprelim$locDepth)
```
```{r echo = FALSE}
load(system.file("extdata", "MsaHindHe.RData", package = "polyRAD"))
```

Now we will look at the distribution of $H_{ind}/H_E$ across samples, and
remove outlier samples.

```{r}
hhByInd <- rowMeans(hh, na.rm = TRUE)

plot(TotDepthT, hhByInd, xlog = TRUE,
     xlab = "Depth", ylab = "Hind/He", main = "Samples")
abline(h = 0.5, lty = 2)
```

For a diploid population, we aren't expecting values above 0.5.  (In general,
you shouldn't see values above $\frac{ploidy - 1}{ploidy}$.)  So, let's check
which are the outlier samples, then remove them from the dataset.

```{r}
threshold <- mean(hhByInd) + 3 * sd(hhByInd)
threshold

hhByInd[hhByInd > threshold]
hh <- hh[hhByInd <= threshold,]
```
```{r eval = FALSE}
myRADprelim <- SubsetByTaxon(myRADprelim, rownames(hh))
```

If you remove any samples like this, you should also make a text file indicating
samples to retain in the analysis, which we will use in the next Python step.

```{r eval = FALSE}
writeLines(rownames(hh), con = "samples.txt")
```

Now we can get a sense of what $H_{ind}/H_E$ looks like across loci.  From this
we can estimate inbreeding and get a rough sense of what proportion of loci
will need to be adjusted or filtered by `process_isoloci.py`.

```{r}
hhByLoc <- colMeans(hh, na.rm = TRUE)

hist(hhByLoc, breaks = 50, xlab = "Hind/He", main = "Loci", col = "lightgrey")
```

Again, we expect values below 0.5 (*i.e.* $\frac{ploidy - 1}{ploidy}$).  The
tail that we observe above that value likely represents groups of tags that
aligned to the same location but in fact represent paralogous loci.  The values
near zero likely represent loci with very high amplification bias or
overdispersion.  We will take the peak at 0.3 to represent typical well-behaved
markers.

To estimate inbreeding, we'll also need to estimate overdispersion.  You may
want to check a broader range of values than 9 to 14.
In the matrix of quantiles, we are
looking for columns where the 25th percentile is about 0.25, the 50th
percentile is about 0.5, etc.

```{r eval = FALSE}
overdispersionP <- TestOverdispersion(myRADprelim, to_test = 9:14)

sapply(overdispersionP[names(overdispersionP) != "optimal"],
       quantile, probs = c(0.01, 0.25, 0.5, 0.75, 0.99))
```

```{r echo = FALSE}
cat("Optimal value is 12.\n", sep = "\n")

print(structure(c(0.025793146160144, 0.253340971566736, 0.552888349753057, 
0.84390807760186, 1, 0.019449092419997, 0.241258741258741, 0.541687845916098, 
0.836394612538573, 1, 0.0147776703743338, 0.226466237547336, 
0.523304087800422, 0.829407211200752, 1, 0.0114170501660437, 
0.211650502292279, 0.505828679120272, 0.82350760103819, 1, 0.00893782083692289, 
0.197336617073809, 0.490038344650588, 0.818544664007547, 1, 0.00697809321403532, 
0.185307346326836, 0.4777134652488, 0.813849077168056, 1), .Dim = 5:6, .Dimnames = list(
    c("1%", "25%", "50%", "75%", "99%"), c("9", "10", "11", "12", 
    "13", "14"))))
```

If $H_{ind}/H_E$ is 0.3 and overdispersion is 12, we can estimate inbreeding
to be about 0.4 using the figure below, which was generated using
simulated data.

![Estimation of inbreeding](overdispersion_inbreeding.png)

If all markers were Mendelian, we can see what the expected distribution of
$H_{ind}/H_E$ would look like, given the sample size and read depth distribution
observed in the dataset.

```{r eval = FALSE}
ExpectedHindHe(myRADprelim, inbreeding = 0.4, ploidy = 2, overdispersion = 12)
```
```{r echo = FALSE}
message("Simulating rep 1")
message("Completed 5 simulation reps")
cat(c("Mean Hind/He: 0.286",
      "Standard deviation: 0.0891",
      "95% of observations are between 0.144 and 0.481"), sep = "\n")
load(system.file("extdata", "MsaHindHe3.RData", package = "polyRAD"))
hist(testhhdist, xlab = "Hind/He", main = "Expected distribution of Hind/He",
     breaks = 30)
```

## Sorting tags into isoloci and filtering isoloci

We are now ready to run the script `process_isoloci.py`.  It will be run
individually on each pair of "align" and "depth" files generated by
`process_sam_multi.py`.  Note that `isoloci_fun.py` must be present in the same
directory as `process_isoloci.py`.  We'll use an expected $H_{ind}/H_E$ of
0.3 and a maximum of 0.5, rounded from the output of `ExpectedHindHe` above.

```
python process_isoloci.py split_depths/Msa_split_1_align.csv
split_depths/Msa_split_1_depth.csv -e 0.3 -m 0.5 -s samples.txt
```

This will take some time to process, and will output a file called
`Msa_split_1_sorted.csv` in the `split_depths` folder.  This file will contain
corrected alignment results, tag sequences, position and sequence for the
variable portion of the tag with respect to the reference genome, and depth
data for all retained samples.  Tags are arranged into isoloci (*i.e.*
assigned to alignment locations) to optimize $H_{ind}/H_E$ first and number of
mutations from the genome second.  Isoloci are removed from the dataset if
they exceed the expected $H_{ind}/H_E$ by a certain amount (essentially, if
it is closer to what would be expected for double the ploidy than for the
input ploidy).

## Genotype calling

We can now import the sorted dataset into polyRAD.  Since there is a small
number of individuals in this dataset, we will lower the filtering thresholds
for loci we import.  The `possiblePloidies` argument is left at the default of
`2` because although the species is allotetraploid, the sorted loci should now
behave in a diploid fashion.

```{r eval = FALSE}
myRAD <- readProcessIsoloci("Msa_split_1_sorted.csv", min.ind.with.reads = 80,
                            min.ind.with.minor.allele = 5)
```

We can see what the distribution of $H_{ind}/H_E$ looks like after we have
sorted and filtered isoloci.

```{r eval = FALSE}
hh2 <- HindHe(myRAD)
hh2ByInd <- rowMeans(hh2, na.rm = TRUE)
hh2ByLoc <- colMeans(hh2, na.rm = TRUE)
```
```{r echo = FALSE}
load(system.file("extdata", "MsaHindHe2.RData", package = "polyRAD"))
```

```{r}
hist(hh2ByInd, xlab = "Hind/He", main = "Samples", breaks = 20, col = "lightgrey")
hist(hh2ByLoc, xlab = "Hind/He", main = "Loci", breaks = 50, col = "lightgrey")
```

By default, `readProcessIsoloci` runs `MergeRareHaplotypes` internally.  (This
behavior can be turned off with `mergeRareHap = FALSE`.)  Those
loci that now appear to have a very inflated $H_{ind}/H_E$ are ones where all
true minor alleles were rare enough to be merged into the common allele, and
any remaining minor alleles only represent sequencing error or alleles that were
amplified very poorly.  We can get rid of these.

```{r}
mean(hh2ByLoc <= 0.5) # proportion of loci retained
keeploci <- names(hh2ByLoc)[hh2ByLoc <= 0.5]
head(keeploci)
hist(hh2ByLoc[keeploci], xlab = "Hind/He", main = "Loci", breaks = 50, col = "lightgrey")
```
```{r eval = FALSE}
myRAD <- SubsetByLocus(myRAD, keeploci)
```

Genotype calling can then be performed as normal.

```{r eval = FALSE}
myRAD <- IteratePopStruct(myRAD, overdispersion = 12)
```

If you want genotypes output as traditional SNP markers, you can use the
`RADdata2VCF` function.

```{r eval = FALSE}
RADdata2VCF(myRAD, file = "Msa_test.vcf")
```

