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

## Packages for this tutorial

If you don't already have the following packages, install them:

```{r eval = FALSE}
CRAN_packages <- c("polyRAD", "BiocManager", "qqman")
Bioconductor_packages <- c("pcaMethods", "VariantAnnotation")

BiocManager::install(Bioconductor_packages)
install.packages(CRAN_packages)
```

Then load them into your R workspace:

```{r message = FALSE, warning = FALSE}
library(polyRAD)
library(VariantAnnotation)
library(qqman)
```

As another quick setup item, I'll set the seed for random number generation
to ensure that your results are identical to mine.

```{r}
set.seed(50819)
```

## Import the VCF

The VCF for this tutorial was generated by the
[GBSv2](https://bitbucket.org/tasseladmin/tassel-5-source/wiki/Tassel5GBSv2Pipeline)
pipeline in TASSEL 5.  A large diversity panel of *Miscanthus sacchariflorus*
was sequenced using a RAD-seq approach with the enzymes *Pst*I and *Msp*I.
To make a file small enough for this tutorial,
[VariantAnnotation](http://bioconductor.org/packages/release/bioc/html/VariantAnnotation.html)
was used to filter the dataset by missing data rate and minor allele frequency,
as well as subset it to only contain tetraploid individuals, and only SNPs on
chromosome 5.  The *Miscanthus sinensis* v7 reference from Phytozome was used.

We can get some basic information about the VCF, including the names of the
268 samples:

```{r}
vcf_file <- "Msa_tetraploids_Chr05.vcf.bgz"
hdr <- scanVcfHeader(vcf_file)
hdr
samples(hdr)
```

The VCF can be imported into a `RADdata` object for genotype calling in polyRAD.
There can be a lot of missing data using RAD-seq, so we'll set the filtering
threshold to only retain markers where at least 140 individuals (slightly more
than half) have non-missing data.  We'll also say that we want at least five
individuals to have the minor allele, in order to prevent false positives that
can result from individuals with extreme phenotypes being the only ones to
posess a rare allele.  Using the `possiblePloidies` argument, we will indicate
that the organism is autotetraploid.  Since I know that this is a small
dataset, I am also setting `expectedAlleles` and `expectedLoci` below their
default values to save computation time and memory.

```{r cache = TRUE}
myRAD <- VCF2RADdata(vcf_file,
                     min.ind.with.reads = 140,
                     min.ind.with.minor.allele = 5,
                     possiblePloidies = list(4),
                     expectedAlleles = 5e4,
                     expectedLoci = 2e4)
myRAD
```

We can see that we imported 5097 loci.  The default behavior of `VCF2RADdata` is
also to phase SNPs that appeared to originate from the same tag location,
assuming a reduced-representation sequencing technique.  This means that for
many loci, we will have multiallelic haplotype-based markers rather than biallelic
SNPs.  The following code gives us a sense of how many alleles per locus we tend
to have.

```{r}
table(table(myRAD$alleles2loc))
```

Most markers are biallelic or triallelic, but there are many with more alleles.

We can also get a quick picture of population structure.

```{r cache = TRUE}
myRAD <- AddPCA(myRAD)
plot(myRAD)
```

There appear to be three major groups, which makes sense given
[previous findings](https://doi.org/10.1093/aob/mcy161) indicating
mainland Asian tetraploids, S Japan tetraploids, and N Japan tetraploids
being distinct populations.  We can see which individuals are in which group.

```{r}
sort(GetTaxa(myRAD)[myRAD$PCA[,1] < -5]) # N Japan
sort(GetTaxa(myRAD)[myRAD$PCA[,1] > -5 & myRAD$PCA[,1] < 1]) # S Japan
sort(GetTaxa(myRAD)[myRAD$PCA[,1] > 1]) # mainland and a few from Japan
```

## Preliminary genotype calling and estimation of overdispersion

We will do a quick genotype estimation so that we can test how overdispersed
the distribution of sequencing reads is.  We may want to adjust the model
if it seems like read depth diverges very far from the expected ratios.

```{r cache = TRUE}
myRAD_pre <- IterateHWE(myRAD, tol = 0.001)
testod <- TestOverdispersion(myRAD_pre, to_test = 7:11)
```

The we can visualize the overdispersion results.

```{r cache = TRUE}
qq(testod[["7"]])
qq(testod[["8"]])
qq(testod[["9"]])
qq(testod[["10"]])
qq(testod[["11"]])
```

It seems that the dots follow the red line most closely for 8 and 9, so we
will keep the default of `overdispersion = 9`.  The lower the value that we
choose, the messier the data are assumed to be.

## Calling the genotypes

Now we will use the population structure pipeline to call the genotypes,
since we know there is population structure.  *Miscanthus* are wild, outcrossing
grasses with very little linkage disequilibrium, so we will use the model
without linkage disequilibrium.

```{r cache = TRUE}
myRAD <- IteratePopStruct(myRAD)
```

We can plot population structure again and see that it has tightened up now
that we have estimated allele dosage more accurately.

```{r}
plot(myRAD)
```

## Exploring the results

We can look at the allele depth matrix, and see that the first four alleles
belong to the same locus.

```{r}
myRAD$alleleDepth[1:10,1:5]
```

polyRAD not only estimated the overall frequency of these alleles in the
population, but also the frequency in the local subpopulation that each
individual originated from.

```{r}
myRAD$alleleFreq[1:4]
myRAD$alleleFreqByTaxa[1:10,1:4]
```

These local allele frequencies were used for estimating local genotype
frequencies under Hardy-Weinberg equilibrium (although we could have also
specified a self-fertilization rate).

```{r}
myRAD$priorProb[[1]][,1:5,1:4]
```

Genotype likelihoods were estimated from allelic read depth for each individual,
under a beta-binomial distribution using our selected overdispersion parameter.

```{r}
myRAD$genotypeLikelihood[[1]][,1:5,1:4]
```

Together, the genotype priors and likelihoods were used to estimate the genotype
posterior probabilities.

```{r}
myRAD$posteriorProb[[1]][,1:5,1:4]
```

From these, we can estimate posterior mean genotypes that give us a rough idea
of allele copy number for each individual and allele.

```{r}
posteriormean <- GetWeightedMeanGenotypes(myRAD, maxval = 4,
                                          omit1allelePerLocus = FALSE)
posteriormean[1:10,1:4]
```

## Exporting the genotypes

We can export to GAPIT format so that we may perform GWAS next.

```{r}
gapit_geno <- ExportGAPIT(myRAD)
```

Since many genotypes are uncertain, `ExportGAPIT` exports the
*posterior mean genotypes*, that is, numeric genotypes averaged across the
posterior probability distribution for allele copy number.  To match the
GAPIT documentation, homozygotes are coded as `0` and `2`, with heterozygotes
(or possible heterozygotes) as intermediate numbers.

```{r}
gapit_geno$GD[1:20,1:6]
gapit_geno$GM[1:5,]
```

We don't see many values near 2, since the common allele is omitted for each
marker.  A value around 0.5 indicates an AAAB genotype, 1 an AABB genotype,
and 1.5 an ABBB genotype.

We'll export this data to CSV.

```{r}
write.csv(gapit_geno$GD, file = "Msa_tetraploids_Chr05_GD.csv",
          row.names = FALSE)
write.csv(gapit_geno$GM, file = "Msa_tetraploids_Chr05_GM.csv",
          row.names = FALSE)
```

For genomic prediction, we will export to `rrBLUP`.  Similarly, posterior mean
genotypes will be exported and common alleles omitted.  Genotypes range from
-1 to 1.

```{r}
rrb_geno <- Export_rrBLUP_Amat(myRAD)
rrb_geno[1:20,1:5]
```

Since this is a matrix, we will save it to an R workspace to simplify import
later.

```{r}
save(rrb_geno, file = "Msa_tetraploids_Chr05_rrb.RData")
```

