From 55484440066724eaa85a04460b33f8a23f0aa48c Mon Sep 17 00:00:00 2001 From: piyalkarum Date: Tue, 19 Mar 2024 12:50:01 +0100 Subject: [PATCH] vst permutation added --- NAMESPACE | 1 + NEWS.md | 3 +++ R/post_detect.R | 30 ++++++++++++++++++++---------- man/vstPermutation.Rd | 13 +++++++++++-- 4 files changed, 35 insertions(+), 12 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index 87bfd6e..1b11067 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -46,6 +46,7 @@ importFrom(graphics,polygon) importFrom(graphics,rasterImage) importFrom(graphics,text) importFrom(qgraph,qgraph) +importFrom(stats,as.dist) importFrom(stats,chisq.test) importFrom(stats,density) importFrom(stats,fisher.test) diff --git a/NEWS.md b/NEWS.md index 28db83c..cbb25f4 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,6 @@ +# rCNV 1.3.0 (third update) +* vstPermutation function added + # rCNV 1.2.0 (second update) * relatedness function optimized * bugs fixed in cpm.normal function diff --git a/R/post_detect.R b/R/post_detect.R index 9e62e69..e48c2d0 100644 --- a/R/post_detect.R +++ b/R/post_detect.R @@ -194,14 +194,16 @@ vst<-function(AD,pops,id.list=NULL,qGraph=TRUE,verbose=TRUE,...){ #' @param qGraph logical. Plot the network plot based on observed Vst values #' (see \code{vst()} help page for more details) #' @param histogram logical. plots the distribution histogram of permuted vst values vs. observed values +#' @param stat numeric. The stat to be plotted in histogram. 1 for Mean Absolute Distance or 2 (\code{default}) for Root Mean Square Distance #' #' @importFrom qgraph qgraph #' @importFrom graphics hist +#' @importFrom stats as.dist #' #' @return Returns a list with observed vst values, an array of permuted vst values and the p-values for the permutation test #' #' -#' @author Jorge Cortés-Miranda (email:), Piyal Karunarathne +#' @author Jorge Cortés-Miranda (email:), Piyal Karunarathne #' #' @examples #' \dontrun{data(alleleINF) @@ -212,7 +214,7 @@ vst<-function(AD,pops,id.list=NULL,qGraph=TRUE,verbose=TRUE,...){ #' vstPermutation(ad,pops=substr(colnames(ad)[-c(1:4)],1,11))} #' #' @export -vstPermutation<-function(AD,pops,nperm=100,histogram=TRUE,qGraph=TRUE){ +vstPermutation<-function(AD,pops,nperm=100,histogram=TRUE,stat=2,qGraph=TRUE){ npop<-length(unique(pops)) message("Permuting Vst") perm_vst<-sapply_pb(1:nperm,function(x){ @@ -226,24 +228,32 @@ vstPermutation<-function(AD,pops,nperm=100,histogram=TRUE,qGraph=TRUE){ obs_vst<-vst(AD, pops, qGraph = qGraph,verbose = TRUE) out_mat<-matrix(0, nrow = npop, ncol = npop) + h_dat<-data.frame(matrix(NA,ncol=2,nrow=nperm)) for(i in 1:nperm){ matrix1=perm_vst[,,i] comparison_matrix <- matrix(0, nrow = nrow(matrix1), ncol = ncol(matrix1)) comparison_matrix[!is.na(matrix1) & !is.na(obs_vst)] <- as.numeric(matrix1[!is.na(matrix1)] > obs_vst[!is.na(obs_vst)]) out_mat<-out_mat+comparison_matrix + h_dat[i,1]<-mean(abs(matrix1), na.rm = TRUE) + h_dat[i,2]<-sqrt(mean(matrix1^2, na.rm = TRUE)) } p_mat<-out_mat/nperm - p_mat[upper.tri(p_mat)] <- NA + colnames(p_mat)<-colnames(obs_vst) + rownames(p_mat)<-rownames(obs_vst) + p_mat<-as.dist(p_mat) if(histogram){ - h_dat<-perm_vst[2,1,1:nperm] - hist(h_dat, main = "Vst Permutation Test Histogram", - xlab = "VST", ylab = "Frequency", - col = "dodgerblue", border = F) - abline(v = obs_vst[2,1], col = 2, lwd = 2) - legend("topright", legend = c("Observed Vst", "Permutation Vst"), - col = c(2, "dodgerblue"), lwd = 2) + #stat<-1 + ob_stat<-c(mean(abs(obs_vst), na.rm = TRUE),sqrt(mean(obs_vst^2, na.rm = TRUE))) + xlims<-range(h_dat[,stat],ob_stat[stat],ifelse(ob_stat[stat]>0,(ob_stat[stat]+ob_stat[stat]/10),(ob_stat[stat]-(ob_stat[stat]/10)))) + hist(h_dat[,stat], main = "Vst Permutation Test Histogram", + xlab = "VST", ylab = "Frequency",xlim = xlims, + col = "dodgerblue", border ="dodgerblue" ) + abline(v = ob_stat[stat], col = 2, lwd = 2) + legend("bottomright", legend = c("Observed Vst", "Permutation Vst"), + col = c(2, "dodgerblue"), lwd = 2, bty="n",xpd=TRUE, horiz=TRUE,inset=c(0,1),cex=0.8) } return(list(observed_vst=obs_vst,permuted_vst=perm_vst,pvalue=p_mat)) } + diff --git a/man/vstPermutation.Rd b/man/vstPermutation.Rd index 04f502b..4331fa4 100644 --- a/man/vstPermutation.Rd +++ b/man/vstPermutation.Rd @@ -4,7 +4,14 @@ \alias{vstPermutation} \title{Run permutation on Vst} \usage{ -vstPermutation(AD, pops, nperm = 100, histogram = TRUE, qGraph = TRUE) +vstPermutation( + AD, + pops, + nperm = 100, + histogram = TRUE, + stat = 2, + qGraph = TRUE +) } \arguments{ \item{AD}{data frame of total allele depth values of SNPs} @@ -16,6 +23,8 @@ Must be the same length as the number of samples in AD} \item{histogram}{logical. plots the distribution histogram of permuted vst values vs. observed values} +\item{stat}{numeric. The stat to be plotted in histogram. 1 for Mean Absolute Distance or 2 (\code{default}) for Root Mean Square Distance} + \item{qGraph}{logical. Plot the network plot based on observed Vst values (see \code{vst()} help page for more details)} } @@ -35,5 +44,5 @@ vstPermutation(ad,pops=substr(colnames(ad)[-c(1:4)],1,11))} } \author{ -Jorge Cortés-Miranda (email:), Piyal Karunarathne +Jorge Cortés-Miranda (email:\href{mailto:jorge.cortes.m@ug.uchile.cl}{jorge.cortes.m@ug.uchile.cl}), Piyal Karunarathne }