## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = TRUE)

## ----eval=FALSE---------------------------------------------------------------
#  if (!require("BiocManager", quietly=TRUE))
#      install.packages("BiocManager")
#  BiocManager::install("RAREsim")

## -----------------------------------------------------------------------------
library(RAREsim)
library(ggplot2)

## -----------------------------------------------------------------------------
# load the target data
data("nvariant_afr")
print(nvariant_afr, row.names  =  FALSE)

## -----------------------------------------------------------------------------
nvar <- fit_nvariant(nvariant_afr)
nvar

## -----------------------------------------------------------------------------
plot(nvariant_afr$n, nvariant_afr$per_kb,
     xlab = 'Sample Size', ylab = 'Variants per Kb')
lines(nvariant_afr$n, (nvar$phi*(nvariant_afr$n^nvar$omega)))


## -----------------------------------------------------------------------------
nvariant(phi = nvar$phi, omega = nvar$omega, N = 8128)

## -----------------------------------------------------------------------------
nvariant(N=8128, pop = 'AFR') 

## -----------------------------------------------------------------------------
nvariant(phi = 0.1638108, omega = 0.6248848, N = 8128)

## -----------------------------------------------------------------------------
19.029*nvariant(nvar$phi, omega = nvar$omega, N = 8128)

## -----------------------------------------------------------------------------
# load the data
data("afs_afr")
print(afs_afr)

## -----------------------------------------------------------------------------
af <- fit_afs(Observed_bin_props = afs_afr)
print(af)

## -----------------------------------------------------------------------------
# Label the target and  fitted data
afs_afr$Data <- 'Target Data'
af$Fitted_results$Data <- 'Fitted Function'

# Combine into one data frame
af_all<-rbind(afs_afr, af$Fitted_results)
af_all$Data <- as.factor(af_all$Data)

#plot
ggplot(data = af_all)+
  geom_point( mapping = aes(x= Lower, y= Prop, col = Data))+
  labs(x = 'Lower Limit of MAC Bin',
             y = 'Propoortion of Variants')


## -----------------------------------------------------------------------------
mac <- afs_afr[,c(1:2)]

## -----------------------------------------------------------------------------
afs(mac_bins = mac, pop = 'AFR')

## -----------------------------------------------------------------------------
afs(alpha = 1.594622, beta =  -0.2846474, b  = 0.297495, mac_bins = mac)

## -----------------------------------------------------------------------------
bin_estimates <- expected_variants(Total_num_var = 865.0634, mac_bin_prop = af$Fitted_results)
print(bin_estimates)

## -----------------------------------------------------------------------------
bin_estimates <- expected_variants(Total_num_var =
                                     19.029*nvariant(phi = 0.1638108,
                                                     omega = 0.6248848, N = 8128),
                                   mac_bin_prop = afs(mac_bins = mac, pop = 'AFR'))
print(bin_estimates)

## -----------------------------------------------------------------------------
sessionInfo()