Nothing
#' Correlation replication test based on multivariate analysis results
#'
#' This function is developed to implement correlation replication test based on MVA or cMVA results
#'
#' @param marker The SNP to be analyzed
#' @param gwa.1 GWAS summary statistics for sample 1, includes A1, A2 and two columns for each trait: beta and se
#' @param gwa.2 GWAS summary statistics for sample 2, includes A1, A2 and two columns for each trait: beta and se
#' @param R.1 Phenotypic correlation matrix for sample 1
#' @param R.2 Phenotypic correlation matrix for sample 2
#' @param traits Traits to be analyzed
#' @param nrep The number of Monte Carlo repetitions
#' @param probs Percentiles of the endpoints of confidence interval
#' @param method The method used for computing correlation coefficient
#' @param plot If the results for making correlation test figure are needed
#'
#' @return The function returns two lists of \code{res}, which includes
#' 1) \code{correlation} Estimated correlation computed from original sample;
#' 2) \code{ci.left} The value at left endpoint of confidence interval;
#' 3) \code{ci.right} The value at right endpoint of confidence interval (Note: If there are only two traits,
#' then the ratio of correlation equals to one is provided instead of ci.left and ci.right),
#' and \code{df.plot}, will be provided if \code{plot = TRUE}, includes
#' 1) \code{traits} The name of traits in analysis;
#' 2) \code{rank.1} The rank of estimated effect sizes in sample 1;
#' 3) \code{rank.2} The rank of estimated effect sizes in sample 2;
#' 4) \code{mean.conc} The mean of concordant pairs in MC generated by the trait;
#' 5) \code{sd.conc} The standard deviation of concordant pairs in MC generated by the trait;
#' 6) \code{se.beta} The standard error of estimated effect sizes computed using inverse-variance weighting.
#'
#' @author Zheng Ning, Xia Shen
#'
#' @references
#' Zheng Ning, Yakov Tsepilov, Sodbo Zh. Sharapov, Alexander K. Grishenko, Masoud Shirali, Peter K. Joshi,
#' James F. Wilson, Yudi Pawitan, Chris S. Haley, Yurii S. Aulchenko, Xia Shen (2018).
#' Multivariate discovery, replication, and interpretation of pleiotropic loci using summary association statistics. \emph{Submitted}.
#'
#' @seealso
#' \code{MultiSummary}
#'
#' @examples
#' \dontrun{
#'
#' data(example.MV.cor.test)
#'
#' ## Six-trait correlation test ##
#' traits <- c("HEIGHT", "BMI", "HIP", "WC", "WHR", "WEIGHT")
#' set.seed(510)
#' MV.cor.test(marker = "rs905938", gwa.1 = example.gwa.1, gwa.2 = example.gwa.2, R.1 = example.R.1,
#' R.2 = example.R.2, traits = traits, nrep = 10000)
#'
#' ## Make correlation correlation test figure ##
#' require(ggplot2)
#' require(cowplot)
#'
#' set.seed(510)
#' res.mv.cor <- MV.cor.test(marker = "rs905938", gwa.1 = example.gwa.1, gwa.2 = example.gwa.2, R.1 = example.R.1,
#' R.2 = example.R.2, traits = traits, nrep = 10000, plot = TRUE)
#' df.plot <- res.mv.cor$df.plot
#'
#' p1 <- ggplot()+
#' geom_point(data=df.plot, mapping=aes(x=rank.1, y=rank.2, color=traits), size=2) +
#' geom_point(data=df.plot, mapping=aes(x=rank.1, y=rank.2, color=traits, size = se.beta), alpha = 0.2) +
#' stat_smooth(data=df.plot, mapping=aes(x=rank.1, y=rank.2), method = "lm", se=FALSE, color="black", size=0.3, fullrange = TRUE) +
#' coord_cartesian(xlim = c(0.5, 6.5), ylim = c(0.5, 6.5)) + xlim(0,200) +
#' scale_size_continuous(range = c(3, 10)) +
#' theme(axis.text=element_text(size=10),
#' axis.title=element_text(size=14,face="bold"),
#' strip.text.x = element_text(size = 16))+
#' theme(axis.title.x=element_blank(),axis.text.x=element_blank(),
#' axis.ticks.x=element_blank(),axis.title.y=element_blank(), legend.position = c(0.8,0.3),
#' legend.background=element_rect(colour='NA', fill='transparent'), legend.key=element_blank(),
#' legend.title=element_text(size=14),
#' legend.text=element_text(size=12), legend.key.size = unit(1.4, 'lines')) +
#' guides(colour = guide_legend(override.aes = list(alpha = 1)), size = FALSE) +
#' scale_colour_discrete(name = "Traits")
#'
#' p2 <- ggplot(data=df.plot, aes(x=rank.1,y=mean.conc)) +
#' coord_cartesian(xlim = c(0.5, 6.5), ylim = c(0, 5.5)) +
#' geom_bar(stat = "identity", aes(fill=traits), width = 0.4) + theme(legend.position="none") + theme(
#' strip.background = element_blank(),
#' strip.text.x = element_blank()
#' ) + geom_errorbar(aes(ymin = mean.conc - sd.conc,ymax = mean.conc + sd.conc), width = 0.1) +
#' theme(axis.title.x=element_blank(),axis.text.x=element_blank(),
#' axis.ticks.x=element_blank(),axis.title.y=element_blank()) +
#' theme(plot.margin = unit(c(0, 0, 0, 0), "cm"))
#'
#' require(cowplot)
#' plot_grid(p1,p2,ncol=1,align = "v", rel_heights = c(2,1))
#'
#' }
#' @aliases MV.cor.test, mv.cor.test
#' @keywords multivariate, replication, correlation test
#'
`MV.cor.test` <- function(marker, gwa.1, gwa.2, R.1, R.2, traits, nrep = 10000, probs = c(0.025,0.975), method = "kendall", plot = FALSE){
if((method %in% c("kendall", "pearson", "spearman")) == FALSE){
stop("The method has to be one of 'kendall', 'pearson' or 'spearman'.")
}
R.1 <- R.1[traits, traits]
R.2 <- R.2[traits, traits]
beta.1 <- as.numeric(gwa.1[marker, paste0(traits, ".beta")])
se.1 <- as.numeric(gwa.1[marker, paste0(traits, ".se")])
cov.beta.1 <- diag(se.1) %*% R.1 %*% diag(se.1)
beta.2 <- as.numeric(gwa.2[marker, paste0(traits, ".beta")])
se.2 <- as.numeric(gwa.2[marker, paste0(traits, ".se")])
cov.beta.2 <- diag(se.2) %*% R.2 %*% diag(se.2)
A1.1 <- gwa.1[marker, "A1"]
A1.2 <- gwa.2[marker, "A1"]
if(A1.1 != A1.2){
beta.2 <- -beta.2
}
cor.sample <- cor(beta.1, beta.2, method = method)
n <- nrep
beta.1.mc <- rmvnorm(n, mean = beta.1, sigma = cov.beta.1)
beta.2.mc <- rmvnorm(n, mean = beta.2, sigma = cov.beta.2)
cor.mc <- numeric(n)
for(i in 1:n){
cor.mc[i] <- cor(beta.1.mc[i,], beta.2.mc[i,], method = method)
}
ci <- quantile(cor.mc, probs = probs)
df.plot <- NULL
if(plot == TRUE && length(traits) > 2){
agree <- function(df,i){
sum(apply(sign(df[,i] - df[,-i]),2, prod) == 1)
}
rank.1 <- rank(beta.1, ties.method = "random")
rank.2 <- rank(beta.2, ties.method = "random")
concordance <- matrix(0,n,length(traits))
for(i in 1:n){
concordance[i,] <- sapply(1:length(traits), FUN = agree, df = rbind(beta.1.mc[i,],beta.2.mc[i,]))
}
mean.conc <- colMeans(concordance)
sd.conc <- apply(concordance, 2, sd)
se.beta <- sqrt(1 / (1/(se.1^2) + 1/(se.2^2)))
df.plot <- cbind(rank.1,rank.2,mean.conc, sd.conc, se.beta)
df.plot <- data.frame(traits, df.plot)
colnames(df.plot) <- c("traits", "rank.1","rank.2", "mean.conc","sd.conc", "se.beta")
}
if(plot == TRUE && length(traits) == 2){
stop("The number of traits should be greater than two to make the figure meaningful.")
}
if(length(traits) == 2){
ratio.1 <- sum(cor.mc) / n
res <- c(cor.sample,ratio.1)
names(res) <- c("correlation", "ratio")
} else{
res <- c(cor.sample,ci)
names(res) <- c("correlation", "ci.left", "ci.right")
}
if(plot == TRUE){
return(list(res = res, df.plot = df.plot))
} else{
return(res)
}
}
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.