Skip to content

Commit

Permalink
change name
Browse files Browse the repository at this point in the history
  • Loading branch information
kelsj committed Aug 3, 2022
1 parent ac31451 commit 4ad3755
Show file tree
Hide file tree
Showing 5 changed files with 17 additions and 17 deletions.
2 changes: 1 addition & 1 deletion LICENSE.txt
@@ -1,2 +1,2 @@
YEAR: 2020
YEAR: 2022
COPYRIGHT HOLDER: Kelsey E Johnson
8 changes: 4 additions & 4 deletions R/sampler.R
Expand Up @@ -27,7 +27,7 @@ sampler = function(dists,ac,alphaRec,betaRec,alphaIBD,alpha0,beta0,piPriors,nite
piPriors = piPriors/sum(piPriors)

#partitions from getPartitions
part = ibdibsR::getPartitions(ac)
part = EVICORD::getPartitions(ac)
nrec_k = part$nrec_k
nibd1_k = part$nibd1_k
nibd2_k = part$nibd2_k
Expand Down Expand Up @@ -59,19 +59,19 @@ sampler = function(dists,ac,alphaRec,betaRec,alphaIBD,alpha0,beta0,piPriors,nite

for (i in 1:niter){
#sample beta for ibd
beta_i[i+1,1] = ibdibsR::sample_beta(c(t1_i[i,],t2_i[i,]) %>% na.omit(), alphas[1],alpha0,beta0)
beta_i[i+1,1] = EVICORD::sample_beta(c(t1_i[i,],t2_i[i,]) %>% na.omit(), alphas[1],alpha0,beta0)
beta_i[i+1,2] = betaRec

#update z,v,t
zvt = ibdibsR::sample_t_z_v(alphas,beta_i[i+1,],pi_chain[i,],dists$distL_cM,dists$distR_cM,ac,nibd1_k,nibd2_k,nrec_k,n,nk)
zvt = EVICORD::sample_t_z_v(alphas,beta_i[i+1,],pi_chain[i,],dists$distL_cM,dists$distR_cM,ac,nibd1_k,nibd2_k,nrec_k,n,nk)
z_chain[i+1,] = zvt$z
v_chain[i+1,] = zvt$v
t1_i[i+1,] = zvt$t[,1]
t2_i[i+1,] = zvt$t[,2]
t3_i[i+1,] = zvt$t[,3]

#update pi
pi_chain[i+1,] = sort(ibdibsR::sample_pi(z_chain[i+1,],piPriors,nk),decreasing=T)
pi_chain[i+1,] = sort(EVICORD::sample_pi(z_chain[i+1,],piPriors,nk),decreasing=T)
}

outlist = list("nVars"=n,"beta"=beta_i,"pi"=pi_chain,"t1"=t1_i,"t2"=t2_i,"t3"=t3_i,"z"=z_chain,"v"=v_chain)
Expand Down
12 changes: 6 additions & 6 deletions README.md
Expand Up @@ -12,9 +12,9 @@ The purpose of this method is to classify variants as likely IBD or non-IBD usin
```R
#install from github
library(devtools)
github_install("kelsj/ibdibsR")
github_install("kelsj/EVICORD")
#load package
library(ibdibsR)
library(EVICORD)
```


Expand All @@ -24,13 +24,13 @@ library(ibdibsR)
1. Load pairwise recombination distances in centimorgans (file format should have 3 columns: varID, distL, distR). The distances should be ordered by variant and by allele pairs. Within each variant, the order fo allele pair recombination distances matters; e.g. for allele count 3, the pairwise distances should be in the order 1-2,1-3,2-3; for allele count 4, the order should be 1-2,1-3,1-4,2-3,2-4,3-4; etc.

```R
dists = ibdibsR::load_dists("sim_dists_ex.txt",header=T,ac=4)
dists = EVICORD::load_dists("sim_dists_ex.txt",header=T,ac=4)
```

2. Run the Gibbs sampler on the pairwise recombination distances. Several parameter values must be provided at this step: the allele count (**ac**), values for alpha and beta for recurrent/non-IBD variants (**alphaRec**, **betaRec**; can be calculated with nonIBD\_param\_values.R, see below), alpha for IBD variants (**alphaIBD**; recommended to try multiple values and see how this affects your results), priors for beta for IBD variants (**alpha0**, **beta0**), vector of priors for pi (**piPriors**, should have as many values as there are categories of k; i.e. 1+number of non-IBD partitions), and the number of iterations to run (**niter**).

```R
ibdibsR::sampler(dists,ac=4,alphaRec=40,betaRec=0.1,alphaIBD=10,alpha0=1,beta0=10,piPriors=c(0.4,0.4,0.2),niter=100,outfile="gibbs_sampler_out.rds") #returns RDS object of output: gibbs_sampler_out.rds
EVICORD::sampler(dists,ac=4,alphaRec=40,betaRec=0.1,alphaIBD=10,alpha0=1,beta0=10,piPriors=c(0.4,0.4,0.2),niter=100,outfile="gibbs_sampler_out.rds") #returns RDS object of output: gibbs_sampler_out.rds
```

3. Visualize results from gibbs_sampler_out.rds. Required parameters: **filename**, e.g. gibbs\_sampler\_out.rds; **ac**, allele count; **tf** integer to thin chains by for visualization. Output plot panels:
Expand All @@ -41,15 +41,15 @@ ibdibsR::sampler(dists,ac=4,alphaRec=40,betaRec=0.1,alphaIBD=10,alpha0=1,beta0=1
4. empirical cumulative density of the fraction of z samples for each variant greater than x (with x=1,2,3,...); i.e., zGT1 = the fraction of posterior samples that the variant is assigned to be non-IBD.

```R
ibdibsR::viz_res(filename="gibbs_sampler_out.rds",ac=4,tf=5,outfile="gibbs_sampler_out_viz.pdf") #saves visualizations in pdf gibbs_sampler_out_viz.pdf
EVICORD::viz_res(filename="gibbs_sampler_out.rds",ac=4,tf=5,outfile="gibbs_sampler_out_viz.pdf") #saves visualizations in pdf gibbs_sampler_out_viz.pdf
```

## Estimating beta for non-IBD variants from data

Posterior values of alpha and beta for the distribution of TMRCAs for non-IBD variants need to be supplied to run the sampler function above. These can be estimated from the pairwise recombination distances of non-IBD allele pairs from multiallelic or known non-IBD variants using the function nonIBD\_param\_values. Requires input of non-IBD variants' pairwise recombination distance (columns: varID, distL, distR). Several parameter values must be provided: the number of non-IBD allele pair distances provided for each variant (**n_nonIBDpairs**), the value of alpha for non-IBD allele pairs (**alpha**); priors for beta (**alpha0**, **beta0**), the number of iterations to run (**niter**), and the output file name (**outfile**).

```R
ibdibsR::nonIBD_param_values(dists,n_nonIBDpairs=3,alpha=40,alpha0=1,beta0=10,niter=500,outfile="sampled_nonIBD_param_values_alpha40") #saves sampled values to RDS file sampled_nonIBD_param_values.rds
EVICORD::nonIBD_param_values(dists,n_nonIBDpairs=3,alpha=40,alpha0=1,beta0=10,niter=500,outfile="sampled_nonIBD_param_values_alpha40") #saves sampled values to RDS file sampled_nonIBD_param_values.rds
```


Expand Down
6 changes: 3 additions & 3 deletions example/sample_nonIBD_params_example.ipynb
Expand Up @@ -4,7 +4,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"# Example for ibdibsR function nonIBD_param_values\n",
"# Example for EVICORD function nonIBD_param_values\n",
"\n",
"## Set up"
]
Expand All @@ -17,10 +17,10 @@
},
"outputs": [],
"source": [
"library(ibdibsR)\n",
"library(EVICORD)\n",
"library(ggplot2)\n",
"library(dplyr)\n",
"setwd(\"./ibdibsR/example\")"
"setwd(\"./EVICORD/example\")"
]
},
{
Expand Down
6 changes: 3 additions & 3 deletions example/sampler_example.ipynb
Expand Up @@ -4,7 +4,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"# Example for ibdibsR sampler"
"# Example for EVICORD sampler"
]
},
{
Expand All @@ -15,10 +15,10 @@
},
"outputs": [],
"source": [
"library(ibdibsR)\n",
"library(EVICORD)\n",
"library(ggplot2)\n",
"library(dplyr)\n",
"setwd(\"./ibdibsR/example\")"
"setwd(\"./EVICORD/example\")"
]
},
{
Expand Down

0 comments on commit 4ad3755

Please sign in to comment.