stuart.Rmd

---
title: "stuart"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{stuart}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---


```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

# stuart

Marie Bourdon

June 2021

## Goal

stuart is a R package which formats results of genotyping. It was developed to analyse data from MUGA arrays (Neogen) for Rqtl analysis, for backcross or F2 crosses, but can be used to analyze data of other laboratory animal strains with other arrays. It allows to filter the markers in arrays, from a genetic point of view. Indeed, markers will be selected depending on their proportion of each genotype, correspondance between F2 or N2 individuals alleles and parental strains alleles, etc.

The examples shown here require the use of dplyr package.


```{r setup}
library(dplyr)
library(stuart)
```


## Annotation files

In order to map the markers on the genome of the individuals, you need to load a table with the position of all markers in the array. The data frame must contain the following columns: `marker` with the markers names, `chr` with the chromosome of each marker, and a column with the position of the marker on the chromosome. For Rqtl analysis, you need to provide positions in cM. The data frame can contain other columns that you judge helpful.

The developer of Rqtl and Rqtl2 packages, Karl Broman, realised that the annotation of the MUGA arrays was not correct for some markers. Thus, he produced new annotation files for MUGA, miniMUGA, megaMUGA and gigaMUGA arrays. These files contain some informations about the markers including the chromosome and position where the probe of the marker matchs on the genome, wether the marker maps uniquely or not, etc. These files also contains the genetic position of the markers calculated with two methods : "cM_cox" and "cM_g2f1" (see https://kbroman.org/MUGAarrays/mini_revisited.html for more informations).

We recommand to use these annotation files to reconstruct the file use for Rqtl analysis. You can load the datasets with these annotations from GitHub (https://github.com/kbroman/MUGAarrays/tree/master/UWisc). Choose the file corresponding to the MUGA array that you used and use the URL to load the dataset in R.


```{r annot}
annot_mini <- read.csv(url("https://raw.githubusercontent.com/kbroman/MUGAarrays/master/UWisc/mini_uwisc_v2.csv"))
```


Here, we will present an example of the use of stuart with results of a F2 cross genotyped with miniMUGA. Examples of genotypes and phenotypes dataframe are available in stuart package. 

The genotype data frame must contain a first column with marker names, a second column with sample IDs, a third column with the first allele and a fourth column with the second allele. This format corresponds to the MUGA results. If your data differ, make sure to have these columns in this order.

You can load your genotype data with base R function `read.table()` or with `readr::read_tsv()` function. You have to skip the first lines with the header.

* `read.table("path/to/genotype_data.txt",sep = "\t",skip=9,header=TRUE)`

* `readr::read_tsv("path/to/genotype_data.txt",skip=9)`

The phenotype data produced by the lab must have the following format: one line per individual, one column per phenotype. The first column must be the identification number of each individual.

You can load your phenotype data depending on its format with `read.table()` or `readr::read_tsv()` for txt files, `read.table()` or `readr::read_csv()` for csv files, `read_excel::read_xls` or `read_excel::read_xlsx` for excel files. 

Here, we load the result of Neogen genotyping: `genos` (only useful columns with marker name, sample ID and alleles were kept) and the phenotype dataset produced by the lab: `phenos`, directly available in the package.




```{r load}
data(genos)
summary(genos)
data(phenos)
summary(phenos)
```

### Genotyping of parental strains

To use genotyping result for Rqtl analysis, we need to recode the genotypes of the individuals (originally encoded in A, T, G, C) depending on the genotype of the parental strains: homozygous for the first parental strain (0), heterozygous (1) or homozygous for the second parental strain (2).

We recommend to always genotype the parental strains of the cross. Here, their genotypes are in the `genos` file and correspond to the Sample.ID "StrainsA_1", "StrainsA_2", "StrainsB_1" and "StrainsB_2". Two individuals were genotyped for each parental strain. The first step will be to create a consensus genotype for each strain from the two genotyped individuals. The consensus genotype will be added to the annotation dataset in order to obtain a dataset with both annotation and reference genotype of the parental strains that will be used for recoding the genotypes or the F2 individuals. 

This is done with the `geno_strains()` function. If parental genotypes was in another dataset than the one with second generation individuals' geotypes, from a previous genotyping result for example, this dataset would have been used in this function, indicated with the `geno` argument. If you use reference genotypes data for your parental strain, you must create a table with genotypes of parental strains and marker informations by merging the `annot_mini` table with reference genotypes table. This can be done with `merge()` or `dplyr::full_join()` functions.

```{r strains}
parental_strains <- tibble::tibble(parent1 = c("StrainsA_1","StrainsA_2"),
                                   parent2 = c("StrainsB_1","StrainsB_2"))


strains <- geno_strains(annot=annot_mini,geno=genos,
                        strn=parental_strains,cols=c("chr","cM_cox"))
head(strains) %>% print.data.frame()
```

After this step, we need to remove the genotyping result for these individuals from the `genos` dataset.

```{r no_parent}
genos <- genos %>% filter(!Sample.ID %in% c("StrainsA_1", "StrainsA_2", 
                                            "StrainsB_1","StrainsB_2"))
```


## Markers sorting

### Marker tab

The first step of the markers sorting is to create the marker dataframe with the tab_mark() function. This dataframe contains for each marker its chromosome (`chr`) and position (`cM_cox` was chosen from annot_mini, but you can choose the column that you want provided that it is in the annotation file), the two alleles that can be found in the F2/N2 population (`Allele_1` and `Allele_2`), the number of individuals for each genotype (homozygous for each allele (`n_HM1` and `n_HM2`) and heterozygous (`n_HT`)), and the number of non genotyped individuals (`n_NA`) This step can take several minutes, so you can also load the example output of this function.


```{r tab_mark}
# how to use the function:
# stuart_tab <- tab_mark(genos,annot_mini,"cM_cox")
data(stuart_tab)
summary(stuart_tab)
```

Then we will use the different mark_* functions in order to filter the markers. 

### mark_na

First, we can use `mark_na()` function in order to remove markers with high proportion of missing genotypes in our F2 individuals.

```{r mark_na}
tab2 <- mark_na(stuart_tab)
```



### mark_match

Then, we can use `mark_match()` function. Here, the parental strains were genotyped with the F2 individuals, but it can happen that you use previous genotyping results for the parental strains. `mark_match()` function excludes markers that are in your genotype file but not in the reference genotype dataset. We recomend using this function as the chip used for genotyping may change.

```{r mark_match}
tab2 <- mark_match(tab2,ref=strains)
tab2 %>% filter(exclude_match==1) %>% print.data.frame()
```

Here the reference strains were genotyped with the same version of the chip as the F2 individuals so no marker was excluded.

### mark_poly

Then, we can use the `mark_poly()` function, which will exclude the markers that are not polymorphic.

```{r mark_poly ex}
tab2 <- mark_poly(tab2)
head(tab2) %>% print.data.frame()
```

### mark_prop

The `mark_prop()` function can be used to filter markers depending on the proportion of each genotype. Here, we have a F2 and we use `homo=0.1, hetero=0.1`  so the function will exclude all markers with less than 10% of each genotype. For chromosome X, we use the `homo1X`, `homo2X` and `heteroX` arguments. The expected proportion of each genotype for markers on X chromosome depend on the directions of the cross performed in order to obtain F1 individuals.

For example, in an F2 between a strain A (a/a for all loci) and B (b/b for all loci), if we cross a female A (i.e. a/a for markers on X chromosome) with a male B (i.e. b/Y for markers on X chromosome), we obtain F1 females with a/b genotype and F1 males with a/Y for markers on X chromosome. Crossed together, they give birth to F2 females that can be either be a/a or a/b for X chromosome markers.

On the other hand, if we cross a female B (i.e. b/b for markers on X chromosome) with a male A (i.e. a/Y for markers on X chromosome), we still obtain F1 females with a/b genotype but F1 males are b/Y for these markers. The obtained F2 females can be either a/b or b/b for X chromosome markers. 

This example shows that depending on the directions of the crosses, genotypic proportions for markers on X chromosome can vary. In some cases, it is possible that a genotype cannot be present in the second generation individuals. You must determine these expected genotypic proportions. Moreover, it is possible that all second generation individuals were not obtained thanks to crosses with the same direction; in such cases, it can be quite complex to estimate the expected genotypic proportions.

However, it is still important to clean markers on X chromosome. The `homo1X`, `homo2X` and `heteroX` arguments are vectors of 2 numbers which are the lower and upper limits of the estimated genotypic proportion (`homo1X` for the homozygosity for the allele of parent 1, `homo2X` for the homozygosity for the allele of parent 2 and `heteroX` for the heterozygosity). As we saw previously, estimate these limits can be quite tough. Thus, if all genotypes are possible, we recommand to use the minimal limits `c(0.1,1)` for the 3 argumetns in order to eliminate most of the unwanted markers.

```{r mark_prop_ex_homo}
tab2 <- mark_prop(tab2,cross="F2",homo=0.1,hetero=0.1,homo1X=c(0.1,1),homo2X=c(0.1,1),heteroX=c(0.1,1))
head(tab2) %>% print.data.frame()
```
We could also use the `pval` argument which allows to exclude autosomal markers by performing a Chi2 test comparing observed distribution with Mendelian proportions. By using `pval=0.5` we would exclude all markers with a Chi2 p-value inferior to 0.05. However, for some markers, Chi2 approximation may be incorrect.

```{r mark_prop_ex_pval,warning=FALSE}
mark_prop(tab2,cross="F2",pval=0.05,homo1X=c(0.1,1),homo2X=c(0.1,1),heteroX=c(0.1,1)) %>% head() %>% print.data.frame()
```

### mark_allele

Last, we can use the `mark_allele()` function. This function excludes markers for which the alleles found in the F2/N2 individuals do not correspond to the alleles found in the parental strains. For example, if a marker is not polymorphic in the parental strains but we found two alleles in the F2/N2 individuals, it will be excluded. Moreover, for a backcross, this function will verify that the homozygous N2 individuals are homozygous for the allele of the strain used to backcross F1 individuals. This strain must be indicated as par1.

```{r mark_allele}
tab2 <- mark_allele(tab=tab2,ref=strains,cross="F2",par1="parent1",par2="parent2")
tab2 %>% arrange(desc(exclude_allele)) %>% head() %>% print.data.frame()
```

Indeed, we can see that the markers excluded with `mark_allele()` have different alleles in the parental strains.

```{r mark_allele-strains}
strains %>% filter(marker %in% c("gJAX00001099","gUNC83675","gUNC66010","mUNC010010509","gUNC30322","mUNC010515443")) %>% 
  select(marker,parent1,parent2) %>% print.data.frame()
```

## Creation of the R/qtl file

After excluding the problematic markers, we can create the R/qtl file. The individuals must have the same ID in the geno and in the pheno file. If there is a prefix in the geno file that must be removed in order to acheive this, you can use the "prefix" argument. The "path" argument can be used in order to create a CSV file that you can laod with `qtl::read.cross`. 

```{r write_rqtl}
stuart_output <- write_rqtl(geno=genos,pheno=phenos,tab=tab2,
                        ref=strains,par1="parent1",par2="parent2",prefix="ind_",pos="cM_cox")

stuart_output[1:10,1:7] %>% print.data.frame()
```

## Check markers with est.map()

The best way to verify that the markers are correctly sorted is to calculate the estimated map with R/qtl. Indeed, if there is a problematic marker, it will present high recombination fraction with adjacent markers.

The CSV file produced with stuart must be loaded in R as a cross with `qtl::read.cross`. The dataframe produced by `write_rqtl` cannot be used with R/qtl as it does not have the "cross" class.

The cross object was saved in stuart. Here we can load it as well as the newmap that was produced with the function `est.map(cross=stuart_cross,error.prob=0.01)`.

```{r load_cross}
library(qtl)

data(stuart_cross)
summary(stuart_cross)

data(stuart_newmap)
```


Then we can plot the newmap to compare it with the positions found in the annotation file.

```{r plot_map}
plotMap(stuart_cross,stuart_newmap,shift=TRUE)
```

### mark_estmap

It seems to be still 2 problematic markers. We can identify them with the `mark_estmap()`function.

```{r mark_estmap}
tab2 <- mark_estmap(tab=tab2,map=stuart_newmap,annot=annot_mini)
tab2 %>% filter(exclude_estmap == 1) %>% print.data.frame()
```
We can see that these two markers have correct proportions of each genotypes so they were not excluded. Let's see the alleles of the parental strains for these markers.


```{r bad_markers_annot}
strains %>% filter(marker %in% c("S6J010381992","SS6071602326","S6J205609960"))
```

We can see that for 2 of these markers, one of the two parents has a missing genotype. In the cases where one parent has a missing of a heterozygous genotype, the `mark_allele()` function does not by default exclude the marker if the observed allele of the genotyped parent correspond to one of the alleles observed in the crossed individuals.

### Exclusion of last problematic markers

To avoid the issues with these situations, you can use the `parNH=FALSE` argument in the `mark_allele()` function. This will exclude all the markers with one parent with a missing or a heterozygous genotype. However, this will exclude a very large number of markers. More than 1000 are in this situation in this cross, but about 700 were excluded with other filters. If we exclude all these markers with `parNH=FALSE` we would lose about 300 markers for only 2 problematic ones.

We rather advise you to keep these markers and identify possible problematic markers with `mark_estmap()` and exlude them. 


If there are still some markers that need to be excluded after all these steps, you have two possibilities to exclude them. First, you have to identify the name of these markers by looking at the `newmap` object as follow : in RStudio, click on this objectin the Environment to view it, then click on the arrow on the left of the desired chromosome. Each chromosome is a named vector that is visible here in two columns: one with the name of the marker, and one with the calculated position of the marker. This is how you can detect markers that present a huge gap with the previous and the following marker.

You can see an example for the stuart_newmap object in the following screen shot:

![screenchot_newmap](screenshot_newmap.png)

You can drop these markers in the cross object with `qtl::drop.markers()`, or if you prefere to have a correct CSV object that can always be used for QTL analysis, you can delete these markers in the marker tab and create a new cross file with `write_rqtl()`.