Install and load packages
library(RADstackshelpR)
Optimize ‘m’
#optimize_m function will generate summary stats on your 5 iterative runs
#input can be full path to each file, or just the file name if the vcf files are in your working directory
m.out<-optimize_m(m3="/Users/devder/Desktop/phil.dicaeum/m3.vcf",
m4="/Users/devder/Desktop/phil.dicaeum/m4.vcf",
m5="/Users/devder/Desktop/phil.dicaeum/m5.vcf",
m6="/Users/devder/Desktop/phil.dicaeum/m6.vcf",
m7="/Users/devder/Desktop/phil.dicaeum/m7.vcf")
#visualize depth of coverage
vis_depth(output = m.out)
## [1] "Visualize how different values of m affect average depth in each sample"
## Picking joint bandwidth of 23.7
#visualize the effect of varying m on the number of SNPs retained
vis_snps(output = m.out, stacks_param = "m")
## Visualize how different values of m affect number of SNPs retained.
## Density plot shows the distribution of the number of SNPs retained in each sample,
## while the asterisk denotes the total number of SNPs retained at an 80% completeness cutoff.
## Picking joint bandwidth of 5850
#visualize the effect of varying m on the number of loci retained
vis_loci(output = m.out, stacks_param = "m")
## Visualize how different values of m affect number of polymorphic loci retained.
## Density plot shows the distribution of the number of loci retained in each sample,
## while the asterisk denotes the total number of loci retained at an 80% completeness cutoff. The optimal value is denoted by red color.
## Picking joint bandwidth of 1270
Optimize ‘M’
#optimize_bigM function will generate summary stats on your 5 iterative runs
#input can be full path to each file, or just the file name if the vcf files are in your working directory
bigM.out<-optimize_bigM(M1="/Users/devder/Desktop/phil.dicaeum/bigM1.vcf",
M2="/Users/devder/Desktop/phil.dicaeum/bigM2.vcf",
M3="/Users/devder/Desktop/phil.dicaeum/bigM3.vcf",
M4="/Users/devder/Desktop/phil.dicaeum/bigM4.vcf",
M5="/Users/devder/Desktop/phil.dicaeum/bigM5.vcf",
M6="/Users/devder/Desktop/phil.dicaeum/bigM6.vcf",
M7="/Users/devder/Desktop/phil.dicaeum/bigM7.vcf",
M8="/Users/devder/Desktop/phil.dicaeum/bigM8.vcf")
#use this function to visualize the effect of varying 'M' on the number of SNPs retained
vis_snps(output = bigM.out, stacks_param = "M")
## Visualize how different values of M affect number of SNPs retained.
## Density plot shows the distribution of the number of SNPs retained in each sample,
## while the asterisk denotes the total number of SNPs retained at an 80% completeness cutoff.
## Picking joint bandwidth of 6290
#visualize the effect of varying 'M' on the number of polymorphic loci retained
vis_loci(output = bigM.out, stacks_param = "M")
## Visualize how different values of M affect number of polymorphic loci retained.
## Density plot shows the distribution of the number of loci retained in each sample,
## while the asterisk denotes the total number of loci retained at an 80% completeness cutoff. The optimal value is denoted by red color.
## Picking joint bandwidth of 1230
Optimize ‘n’
#optimize_bigM function will generate summary stats on your 5 iterative runs
#input can be full path to each file, or just the file name if the vcf files are in your working directory
n.out<-optimize_n(nequalsMminus1="/Users/devder/Desktop/phil.dicaeum/n7.vcf",
nequalsM="/Users/devder/Desktop/phil.dicaeum/n8.vcf",
nequalsMplus1="/Users/devder/Desktop/phil.dicaeum/n9.vcf")
#use this function to visualize the effect of varying 'n' on the number of SNPs retained
vis_snps(output = n.out, stacks_param = "n")
## Visualize how different values of n affect number of SNPs retained.
## Density plot shows the distribution of the number of SNPs retained in each sample,
## while the asterisk denotes the total number of SNPs retained at an 80% completeness cutoff.
## Picking joint bandwidth of 8570
#visualize the effect of varying 'n' on the number of polymorphic loci retained
vis_loci(output = n.out, stacks_param = "n")
## Visualize how different values of n affect number of polymorphic loci retained.
## Density plot shows the distribution of the number of loci retained in each sample,
## while the asterisk denotes the total number of loci retained at an 80% completeness cutoff. The optimal value is denoted by red color.
## Picking joint bandwidth of 1350
make summary figure
#load gridExtra package to combine ggplot visualizations
library(gridExtra)
#combine all of these prior visulizations in a single list
gl<-list()
gl[[1]]<-vis_depth(output = m.out)
## [1] "Visualize how different values of m affect average depth in each sample"
gl[[2]]<-vis_snps(output = m.out, stacks_param = "m")
## Visualize how different values of m affect number of SNPs retained.
## Density plot shows the distribution of the number of SNPs retained in each sample,
## while the asterisk denotes the total number of SNPs retained at an 80% completeness cutoff.
gl[[3]]<-vis_loci(output = m.out, stacks_param = "m")
## Visualize how different values of m affect number of polymorphic loci retained.
## Density plot shows the distribution of the number of loci retained in each sample,
## while the asterisk denotes the total number of loci retained at an 80% completeness cutoff. The optimal value is denoted by red color.
gl[[4]]<-vis_snps(output = bigM.out, stacks_param = "M")
## Visualize how different values of M affect number of SNPs retained.
## Density plot shows the distribution of the number of SNPs retained in each sample,
## while the asterisk denotes the total number of SNPs retained at an 80% completeness cutoff.
gl[[5]]<-vis_loci(output = bigM.out, stacks_param = "M")
## Visualize how different values of M affect number of polymorphic loci retained.
## Density plot shows the distribution of the number of loci retained in each sample,
## while the asterisk denotes the total number of loci retained at an 80% completeness cutoff. The optimal value is denoted by red color.
gl[[6]]<-vis_snps(output = n.out, stacks_param = "n")
## Visualize how different values of n affect number of SNPs retained.
## Density plot shows the distribution of the number of SNPs retained in each sample,
## while the asterisk denotes the total number of SNPs retained at an 80% completeness cutoff.
gl[[7]]<-vis_loci(output = n.out, stacks_param = "n")
## Visualize how different values of n affect number of polymorphic loci retained.
## Density plot shows the distribution of the number of loci retained in each sample,
## while the asterisk denotes the total number of loci retained at an 80% completeness cutoff. The optimal value is denoted by red color.
#visualize each item of the list as part of a single grid
grid.arrange(grobs = gl, widths = c(1,1,1,1,1,1),
layout_matrix = rbind(c(1,1,2,2,3,3),
c(4,4,4,5,5,5),
c(6,6,6,7,7,7))
)
## Picking joint bandwidth of 23.7
## Picking joint bandwidth of 5850
## Picking joint bandwidth of 1270
## Picking joint bandwidth of 6290
## Picking joint bandwidth of 1230
## Picking joint bandwidth of 8570
## Picking joint bandwidth of 1350
#remotes::install_github("idmn/ggview")
#library(ggview)
#ggview(units="in", width=4, height=4)
#save
g<-arrangeGrob(grobs = gl, widths = c(1,1,1,1,1,1),
layout_matrix = rbind(c(1,1,2,2,3,3),
c(4,4,4,5,5,5),
c(6,6,6,7,7,7))
)
## Picking joint bandwidth of 23.7
## Picking joint bandwidth of 5850
## Picking joint bandwidth of 1270
## Picking joint bandwidth of 6290
## Picking joint bandwidth of 1230
## Picking joint bandwidth of 8570
## Picking joint bandwidth of 1350
#library(ggplot2)
#ggsave(file="~/Desktop/phil.dicaeum/denovo.optimize.pdf", g, units="in",width=8,height=8) #saves g