Estimating genotype probabilities in a mapping population
Data import example with UNEAK
In this example, we’ll import some data from an F1 mapping population of Miscanthus sinensis that were output by the UNEAK pipeline. These data are from a study by Liu et al. (2015; doi:10.1111/gcbb.12275; data available at http://hdl.handle.net/2142/79522), and can be found in the “extdata” folder of the polyRAD installation. Miscanthus is an ancient tetraploid that has undergone diploidization. Given the ability of the UNEAK pipeline to filter paralogs, we expect most loci to behave in a diploid fashion, but some may behave in an allotetraploid fashion.
We’ll start by loading polyRAD and importing the data into a RADdata object. The possiblePloidies argument indicates the expected inheritance modes: diploid (2) and allotetraploid (2 2).
With your own dataset, you will not need to use system.file. Instead, directly create a text string indicating the name of your file (and its location if it is not in the working directory.)
library(polyRAD)
maphmcfile <- system.file("extdata", "ClareMap_HapMap.hmc.txt",
package = "polyRAD")
maphmcfile## [1] "C:/Users/lvclark/AppData/Local/Temp/RtmpeS3WAw/Rinst3568731d4424/polyRAD/extdata/ClareMap_HapMap.hmc.txt"## ## RADdata object ##
## 299 taxa and 50 loci
## 766014 total reads
## Assumed sample cross-contamination rate of 0.001
##
## Possible ploidies:
## Autodiploid (2)
## Allotetraploid (2 2)We can view the imported taxa names (subsetted here for space).
## [1] "IGR-2011-001" "Kaskade-Justin" "Map1-001" "Map1-002"
## [5] "Map1-003" "Map1-005" "Map1-008" "Map1-011"
## [9] "Map1-016" "Map1-018" "Map1-488" "Map1-489"
## [13] "Map1-490" "Map1-491" "Zebrinus-Justin" "p196-150A-c"
## [17] "p877-348-b"All names starting with “Map” are progeny. “Kaskade-Justin” and “Zebrinus-Justin” are the parents. “IGR-2011-001”, “p196-150A-c”“, and”p877-348-b" aren’t part of the population, but were doubled haploid lines that were used to screen for paralogous markers. We can tell polyRAD which taxa are the parents; since this is an F1 population it doesn’t matter which is “donor” and which is “recurrent”.
mydata <- SetDonorParent(mydata, "Kaskade-Justin")
mydata <- SetRecurrentParent(mydata, "Zebrinus-Justin")The next thing we’ll want to do is add our genomic alignment data. For this dataset, we have alignment data stored in a CSV file, also in the “extdata” directory of the polyRAD installation. We’ll add it to the locTable slot of our RADdata object. Be sure to name the new columns “Chr” and “Pos”.
alignfile <- system.file("extdata", "ClareMap_alignments.csv",
package = "polyRAD")
aligndata <- read.csv(alignfile, row.names = 1)
head(aligndata)## Sorghum.LG Position.on.Sorghum.LG..bp.
## TP5489 1 4560204
## TP13305 1 4584260
## TP18261 1 2911329
## TP18674 1 387849
## TP19030 1 7576879
## TP26698 1 6972841mydata$locTable$Chr <- aligndata[GetLoci(mydata), 1]
mydata$locTable$Pos <- aligndata[GetLoci(mydata), 2]
head(mydata$locTable)## Chr Pos
## TP5489 1 4560204
## TP13305 1 4584260
## TP18261 1 2911329
## TP18674 1 387849
## TP19030 1 7576879
## TP26698 1 6972841If you don’t have alignment data in your own dataset, you can still use the pipeline described here. Just set useLinkage = FALSE in the code below. The advantage of including alignment data is that gentoypes at linked markers are used for imputing missing or correcting erroneous genotypes.
Quality control and genotype calling
It is important that the only individuals included in the analysis are those that are truly part of the population. Allele frequencies are used for inferring segregation pattern, and could be incorrect if many individuals are included that are not part of the population. Additionally, the genotype priors will be incorrect for individuals that are not part of the population, leading to incorrect genotypes.
At this point we would normally do
However, because a very small number of markers was used in this example dataset, the PCA does not accurately reflect the relatedness of individuals. Here I will load a PCA that was done with the full set of markers.
Now a plot can be used for visualizing the relationship among taxa.
Now we’ll extract a subset of taxa that we actually want to analyze. We can see from the plot that a fair number of them were the product of self-fertilization of “Zebrinus-Justin” and should be eliminated.
realprogeny <- GetTaxa(mydata)[mydata$PCA[,"PC1"] > -10 &
mydata$PCA[,"PC1"] < 10]
# eliminate the one doubled haploid line in this group
realprogeny <- realprogeny[!realprogeny %in% c("IGR-2011-001", "p196-150A-c",
"p877-348-b")]
# also retain parents
keeptaxa <- c(realprogeny, GetDonorParent(mydata), GetRecurrentParent(mydata))
mydata <- SubsetByTaxon(mydata, taxa = keeptaxa)
plot(mydata)
Now we can perform a preliminary run of the pipeline. The allowedDeviation argument indicates how different the apparent allele frequency (based on read depth ratios) can be from an expected allele frequency (determined based on ploidy and mapping population type) and still be classified as that allele frequency. The default settings assume an F1 population, but the population type can be adjusted using the n.gen.backcrossing, n.gen.intermating, and n.gen.selfing arguments. We’ll also lower minLikelihoodRatio from the default because one of the parents has many uncertain genotypes under the tetraploid model (which was determined by exploration of the dataset outside of this tutorial; many NA values were observed in priorProb under the default). Since this first run is for a rough estimate of genotypes, we’ll set useLinkage = FALSE to save a little computational time.
mydata2 <- PipelineMapping2Parents(mydata,
freqAllowedDeviation = 0.06,
useLinkage = FALSE,
minLikelihoodRatio = 2)## Making initial parameter estimates...## Generating sampling permutations for allele depth.## Done.We can use these preliminary estimates to determine whether we need to adjust the overdispersion parameter. Exactly how much does read depth distribution deviate from what would be expected under binomial distibution? The TestOverdispersion function will help us here.
## Optimal value is 12.sapply(overdispersionP[names(overdispersionP) != "optimal"],
quantile, probs = c(0.01, 0.25, 0.5, 0.75, 0.99))## 8 9 10 11 12 13
## 1% 0.005533597 0.004096143 0.002848244 0.001954561 0.001322713 0.0009932275
## 25% 0.304361182 0.279742132 0.258070046 0.239189269 0.223981900 0.2093677395
## 50% 0.583948430 0.562335702 0.545448898 0.528202715 0.512373571 0.4977739384
## 75% 0.818619207 0.808643509 0.799476463 0.793000540 0.785103002 0.7775481475
## 99% 1.000000000 1.000000000 1.000000000 1.000000000 1.000000000 1.0000000000
## 14 15
## 1% 0.0007527206 0.0005746429
## 25% 0.1957778953 0.1837745108
## 50% 0.4842515536 0.4706852334
## 75% 0.7704398101 0.7645633921
## 99% 1.0000000000 1.0000000000It looks like 12 follows the expected distribution of p-values most closely, both from the automated output of TestOverdispersion and from examining the quantiles of p-values. 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.
Next we can check for markers that are behaving in a non-Mendelian fashion. If we are expecting diploid segregation, all markers should show a \(H_{ind}/H_E\) value of 0.5 or less. (For an autopolyploid, the expected value is \(\frac{ploidy - 1}{ploidy}\).)
myhindhe <- HindHeMapping(mydata, ploidy = 2L)
hist(colMeans(myhindhe, na.rm = TRUE), col = "lightgrey",
xlab = "Hind/He", main = "Histogram of Hind/He by locus")
How does this compare to the distribution that we might expect from this dataset if all loci were Mendelian? We will use ExpectedHindHeMapping with the overdispersion parameter that we estimated above. One your own dataset, I recommend using the default value of reps.
set.seed(720)
ExpectedHindHeMapping(mydata, ploidy = 2, overdispersion = my_ovdisp, reps = 2,
contamRate = 0.001, errorRate = 0.001)## Generating sampling permutations for allele depth.## Simulating rep 1## Completed 2 simulation reps.## Mean Hind/He: 0.468
## Standard deviation: 0.027
## 95% of observations are between 0.426 and 0.529
Comparing the observed distribution to the expected distribution, we may want to filter some markers.
goodMarkers <- colnames(myhindhe)[which(colMeans(myhindhe, na.rm = TRUE) < 0.53 &
colMeans(myhindhe, na.rm = TRUE) > 0.43)]
mydata <- SubsetByLocus(mydata, goodMarkers)Now we can re-run the pipeline to properly call the genotypes.
mydata <- PipelineMapping2Parents(mydata,
freqAllowedDeviation = 0.06,
useLinkage = TRUE, overdispersion = my_ovdisp,
minLikelihoodRatio = 2)## Making initial parameter estimates...## Generating sampling permutations for allele depth.## Updating priors using linkage...## Done.Examining the output
We can examine the allele frequencies. Allele frequencies that fall outside of the expected ranges will be recorded as they were estimated from read depth. In this case most are within the expected ranges.
##
## 0.25 0.5 0.75
## 11 2 11Genotype likelihood is also stored in the object for each possible genotype at each locus, taxon, and ploidy. This is the probability of seeing the observed distribution of reads.
## TP19030_0 TP19030_1 TP28986_0 TP28986_1 TP31810_0 TP31810_1 TP34939_0 TP34939_1
## 2 13 0 6 0 16 8 9## TP19030_0 TP19030_1 TP28986_0 TP28986_1 TP31810_0 TP31810_1
## 0 9.706321e-04 4.742367e-09 9.987416e-01 1.237115e-07 9.973885e-01 4.444851e-11
## 1 2.328815e-02 2.328815e-02 3.741288e-02 3.741288e-02 1.567148e-03 1.567148e-03
## 2 4.742367e-09 9.706321e-04 1.237115e-07 9.987416e-01 4.444851e-11 9.973885e-01
## TP34939_0 TP34939_1
## 0 8.959640e-06 1.180948e-06
## 1 1.199589e-01 1.199589e-01
## 2 1.180948e-06 8.959640e-06## TP19030_0 TP19030_1 TP28986_0 TP28986_1 TP31810_0 TP31810_1
## 0 9.706321e-04 4.742367e-09 9.987416e-01 1.237115e-07 9.973885e-01 4.444851e-11
## 1 1.578635e-01 6.145133e-04 2.426471e-01 2.279025e-03 5.641026e-02 1.187627e-05
## 2 2.328815e-02 2.328815e-02 3.741288e-02 3.741288e-02 1.567148e-03 1.567148e-03
## 3 6.145133e-04 1.578635e-01 2.279025e-03 2.426471e-01 1.187627e-05 5.641026e-02
## 4 4.742367e-09 9.706321e-04 1.237115e-07 9.987416e-01 4.444851e-11 9.973885e-01
## TP34939_0 TP34939_1
## 0 8.959640e-06 1.180948e-06
## 1 5.115888e-02 3.296284e-02
## 2 1.199589e-01 1.199589e-01
## 3 3.296284e-02 5.115888e-02
## 4 1.180948e-06 8.959640e-06Above, for one individal (Map1-089), we see its read depth at eight alleles (four loci), followed by the genotype likelihoods under diploid and tetraploid models. For example, at locus TP19030, heterozygosity is the most likely state, although there is a chance that this individual is homozygous for allele 1 and the two reads of allele 0 were due to contamination. If this locus is allotetraploid, it is most likely that there is one copy of allele 0 and three copies of allele 1. Other individuals have higher depth and as a result there is less uncertainty in the genotype, particularly for the diploid model.
The prior genotype probabilities (expected genotype distributions) are also stored in the object for each possible ploidy. These distributions are estimated based on the most likely parent genotypes. Low confidence parent genotypes can be ignored by increasing the minLikelihoodRatio argument to PipelineMapping2Parents.
## TP19030_0 TP19030_1 TP28986_0 TP28986_1 TP31810_0 TP31810_1 TP34939_0
## 0 0.25 0.25 0.5 0.0 0.5 0.0 0.0
## 1 0.50 0.50 0.5 0.5 0.5 0.5 0.5
## 2 0.25 0.25 0.0 0.5 0.0 0.5 0.5
## TP34939_1
## 0 0.5
## 1 0.5
## 2 0.0## TP19030_0 TP19030_1 TP28986_0 TP28986_1 TP31810_0 TP31810_1 TP34939_0
## 0 0 0 NA NA 0 0 NA
## 1 0 0 NA NA 1 0 NA
## 2 1 1 NA NA 0 0 NA
## 3 0 0 NA NA 0 1 NA
## 4 0 0 NA NA 0 0 NA
## TP34939_1
## 0 NA
## 1 NA
## 2 NA
## 3 NA
## 4 NAHere we see some pretty big differences under the diploid and allotetraploid models. For example, if TP19030 is behaving in a diploid fashion we expect F2-like segregation since both parents were heterozygous. However, if TP19030 is behaving in an allotetraploid fashion, a 1:1 segregation ratio is expected due to one parent being heterozygous at one isolocus and the other being homozygous at both isoloci.
Now we want to determine which ploidy is the best fit for each locus. This is done by comparing genotype prior probabilities to genotype likelihoods and estimating a \(\chi^2\) statistic. Lower values indicate a better fit.
## TP19030_0 TP19030_1 TP28986_0 TP28986_1 TP31810_0 TP31810_1
## [1,] 0.3258861 0.3258861 0.07787878 0.07787878 0.1627961 0.1627961
## [2,] 162.8754654 162.8754654 NA NA 202.4987169 202.4987169
## TP34939_0 TP34939_1
## [1,] 1.021263 1.021263
## [2,] NA NAWe can make a plot to get an overall sense of how well the markers fit the diploid versus tetraploid model.
plot(mydata$ploidyChiSq[1,], mydata$ploidyChiSq[2,],
xlab = "Chi-squared for diploid model",
ylab = "Chi-squared for tetraploid model")
For each allele, whichever model gives the lower Chi-squared value is the one with the best fit. In this case it looks like everything is diploid with fairly high confidence, in agreement with our \(H_{ind}/H_E\) results.
Now we’ll examine the posterior genotype probabilities. These are still estimated separately for each ploidy.
## TP19030_0 TP19030_1 TP28986_0 TP28986_1 TP31810_0 TP31810_1
## 0 1.969255e-02 9.339240e-11 0.92839995 0.00000000 0.998431214 0.000000000
## 1 9.803075e-01 9.803075e-01 0.07160005 0.07160005 0.001568786 0.001568786
## 2 9.339240e-11 1.969255e-02 0.00000000 0.92839995 0.000000000 0.998431214
## TP34939_0 TP34939_1
## 0 0.000000e+00 1.778922e-07
## 1 9.999998e-01 9.999998e-01
## 2 1.778922e-07 0.000000e+00## TP19030_0 TP19030_1 TP28986_0 TP28986_1 TP31810_0 TP31810_1 TP34939_0
## 0 0 0 NA NA 0 0 NaN
## 1 0 0 NA NA 1 0 NaN
## 2 1 1 NA NA 0 0 NaN
## 3 0 0 NA NA 0 1 NaN
## 4 0 0 NA NA 0 0 NaN
## TP34939_1
## 0 NaN
## 1 NaN
## 2 NaN
## 3 NaN
## 4 NaNResults cleanup and export
Since we decided from the Chi-squared results that the markers were only segregating in a diploid manner, we can remove allotetraploidy from the dataset.
Typically in a mapping population, due to noisy data polyRAD will not be able to determine the segregation patterns of some markers, which end up having NA values for their prior and posterior probabilities. There may also be some cases where both parents were homozygous and as a result there is no segregation in an F1 population. In this example dataset, these issues are not present (as long as diploidy is assumed) because the markers were curated from a set that had already been filtered for mapping. Generally, however, you would want to find and remove such markers using RemoveUngenotypedLoci:
We can export the results for use in downstream analysis. The function below weights possible ploidies for each allele based on the results in mydata$ploidyChiSq, and for each taxon outputs a continuous, numerical genotype that is the mean of all possible genotypes weighted by genotype posterior probabilities (i.e. the posterior mean genotype). By default, one allele per locus is discarded in order to avoid mathematical singularities in downstream analysis. The continuous genotypes also range from zero to one by default, which can be changed with the minval and maxval arguments.
## TP53071_0 TP57018_0 TP110401_1 TP115159_0
## Kaskade-Justin 0.5 0.5 0.5 0.5
## Zebrinus-Justin 0.0 0.0 0.0 0.0
## Map1-001 0.5 0.5 0.0 0.0
## Map1-002 0.0 0.0 0.5 0.5
## Map1-003 0.0 0.0 0.0 0.5
## Map1-005 0.5 0.5 0.5 0.5
## Map1-008 0.0 0.0 0.5 0.0Note that the parent posterior mean genotypes were estimated using gentoype likelihood only, ignoring the priors set for the progeny. In some places they may not match the progeny genotypes, indicating a likely error in parental genotype calling. We can see the parental genotypes that were used for estimating progeny priors using $likelyGeno_donor and $likelyGeno_recurrent.
## TP19030_0 TP19030_1 TP28986_0 TP28986_1 TP31810_0 TP31810_1 TP34939_0 TP34939_1
## 1 1 1 1 0 2 1 1## TP19030_0 TP19030_1 TP28986_0 TP28986_1 TP31810_0 TP31810_1 TP34939_0 TP34939_1
## 1 1 0 2 1 1 2 0
