Nothing
## ---- echo = FALSE------------------------------------------------------------
library(knitr)
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", cache = FALSE,
fig.width = 5, fig.height = 4) # warning = FALSE
## ---- eval = FALSE------------------------------------------------------------
# install.packages("inbreedR")
## ---- eval = FALSE------------------------------------------------------------
# install.packages("devtools")
# devtools::install_github("mastoffel/inbreedR", build_vignettes = TRUE)
## -----------------------------------------------------------------------------
library("inbreedR")
## ---- eval = FALSE------------------------------------------------------------
# ?inbreedR
## ---- echo=TRUE, results='hide'-----------------------------------------------
library(inbreedR)
data("mouse_msats") # microsatellite data
data("mouse_snps") # snp data
data("bodyweight") # fitness data
## -----------------------------------------------------------------------------
data("mouse_snps")
mouse_snps[1:10, 1:10]
## -----------------------------------------------------------------------------
check_data(mouse_snps, num_ind = 36, num_loci = 13198)
## -----------------------------------------------------------------------------
data("mouse_msats")
mouse_msats[1:8, 1:8]
## -----------------------------------------------------------------------------
mouse_microsats <- convert_raw(mouse_msats)
mouse_microsats[1:8, 1:8]
## ---- eval = FALSE------------------------------------------------------------
# # install.packages("vcfR")
# # install.packages("reshape")
# library(vcfR)
# library(reshape2)
# vcf_file <- "yourvcffile.vcf"
# # read vcf
# vcf <- read.vcfR(vcf_file, verbose = FALSE )
# # extract genotypes
# gt <- extract.gt(vcf)
# # transpose and data.frame
# gt <- as.data.frame(t(gt), stringsAsFactors = FALSE)
# # NA handling
# gt[gt == "."] <- NA
# # split columns
# snp_geno <- do.call(cbind, apply(gt, 2, function(x) colsplit(x, "/", c("a","b"))))
# # convert
# mouse_snp_genotypes <- inbreedR::convert_raw(snp_geno)
# # check data
# check_data(mouse_snp_genotypes)
## ---- echo=TRUE, results='hide'-----------------------------------------------
g2_mouse_microsats <- g2_microsats(mouse_microsats, nperm = 100, nboot = 100, CI = 0.95)
g2_mouse_snps <- g2_snps(mouse_snps, nperm = 100, nboot = 10,
CI = 0.95, parallel = FALSE, ncores = NULL)
## -----------------------------------------------------------------------------
g2_mouse_microsats
## ----g2_bootstrapping, echo=TRUE, fig.cap="Distribution of g2 from bootstrapping with confidence interval", fig.width=7, fig.height=4----
par(mfrow=c(1,2))
plot(g2_mouse_microsats, main = "Microsatellites",
col = "cornflowerblue", cex.axis=0.85)
plot(g2_mouse_snps, main = "SNPs",
col = "darkgoldenrod1", cex.axis=0.85)
## ---- echo=TRUE, results='hide'-----------------------------------------------
HHC_mouse_microsats <- HHC(mouse_microsats , reps = 1000)
HHC_mouse_snps <- HHC(mouse_snps, reps = 100)
## -----------------------------------------------------------------------------
HHC_mouse_microsats
## ----HHC, echo=TRUE, fig.cap="Distribution of heterozygosity-heterozygosity correlations", fig.width=7, fig.height=4----
par(mfrow=c(1,2))
plot(HHC_mouse_microsats, main = "Microsatellites",
col = "cornflowerblue", cex.axis=0.85)
plot(HHC_mouse_snps, main = "SNPs",
col = "darkgoldenrod1", cex.axis=0.85)
## ---- echo=TRUE---------------------------------------------------------------
# r^2 between inbreeding and heterozygosity
hf <- r2_hf(genotypes = mouse_microsats, type = "msats")
# r^2 between inbreeding and fitness
Wf <- r2_Wf(genotypes = mouse_microsats, trait = bodyweight,
family = gaussian, type = "msats")
## ---- echo=TRUE, results = "hide"---------------------------------------------
# r^2 between inbreeding and heterozygosity with bootstrapping
hf <- r2_hf(genotypes = mouse_microsats, nboot = 100, type = "msats", parallel = FALSE)
## -----------------------------------------------------------------------------
plot(hf)
## ---- echo=TRUE, results='hide'-----------------------------------------------
# g2
g2 <- g2_microsats(mouse_microsats)
# calculate sMLH
het <- sMLH(mouse_microsats)
# variance in sMLH
het_var <- var(het)
# Linear model of fitness trait on heterozygosity
mod <- lm(bodyweight ~ het)
# regression slope
beta <- coef(mod)[2]
# r2 between fitness and heterozygosity
Wh <- cor(bodyweight,predict(mod))^2
# r2 between inbreeding and heterozygosity
hf <- r2_hf(genotypes = mouse_microsats, type = "msats")
# r2 between inbreeding and fitness
Wf <- r2_Wf(genotypes = mouse_microsats, trait = bodyweight,
family = gaussian, type = "msats")
## ---- echo=FALSE--------------------------------------------------------------
library(knitr)
df_msats <- data.frame(g2[["g2"]], het_var, beta, Wh, hf$r2_hf_full, Wf$r2_Wf_full)
row.names(df_msats) <- c("microsats")
kable(df_msats, digits=3, caption="Descriptors of HFCs",
col.names = c("$\\hat{g}_2$", "$\\hat{\\sigma}^2(h)$", "$\\hat{\\beta}_{Wh}$", "$\\hat{r}^2_{Wh}$", "$\\hat{r}^2_{hf}$","$\\hat{r}^2_{Wf}$"))
## ---- echo=TRUE, results='hide'-----------------------------------------------
sim_g2 <- simulate_g2(n_ind = 20, H_nonInb = 0.5, meanF = 0.2, varF = 0.05,
subsets = c(4,6,8,10,12), reps = 100,
type = "msats", CI = 0.95)
## ----sim_g2, echo=FALSE, fig.cap="Simulation: Sensitivity of g2 estimated from an increasing number of markers", fig.width=7, fig.height=4----
plot(sim_g2, main = "Microsatellites", cex.axis=0.85)
## ----sim_g2_true, echo=FALSE, fig.cap="Simulation: Sensitivity of g2 estimated from an increasing number of markers with true g2 value", fig.width=7, fig.height=4----
plot(sim_g2, true_g2 = TRUE, main = "Microsatellites", cex.axis=0.85)
## ---- results = "hide"--------------------------------------------------------
g2_seals <- g2_microsats(mouse_microsats, nperm = 100,
nboot = 100, CI = 0.95)
## -----------------------------------------------------------------------------
str(g2_seals)
## -----------------------------------------------------------------------------
g2_bootstrap_results <- g2_seals$g2_boot
str(g2_bootstrap_results)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.