RayDebashree
Inquiry on PLACO Method from PLoS Genetics Paper （[email protected]）
I have used [email protected] to send an email to you containing all the questions mentioned earlier. Perhaps you could copy the content of the email response so that it can be visible to everyone on GitHub as well.
I am reaching out to discuss the PLACO technique introduced in your publication: Ray, D., Chatterjee, N. (2020) "A powerful method for pleiotropic analysis under composite null hypothesis identifies novel shared loci between Type 2 Diabetes and Prostate Cancer". PLoS Genetics 16(12): e1009218. I find the technique very useful and have several questions I hope you could clarify:
1.Can PLACO be used to analyze the relationship between a binary trait (such as ovarian cancer with case control) and a continuous variable (such as epigenetic age acceleration) with complete GWAS data, including beta, se, and pval?
2.Can PLACO be used to analyze between one continuous variable (such as accelerated epigenetic aging, with the unit being years) and another continuous variable (such as telomere length, with the unit being standard deviation change in log-transformed telomere length), both having complete GWAS data including beta, standard error (se), and p-value (pval)?
3.After de-correlating the Z-scores, can PLACO be used on datasets with substantial sample overlap, or almost complete overlap, such as two datasets from the UK Biobank?
4.In your GitHub, you mentioned that PLACO should not be used for highly correlated GWAS, stating "PLACO does not work well if the two traits are strongly correlated." Does this conflict with the approach taken in the paper published in JAMA Psychiatry titled "Role of the Gut-Brain Axis in the Shared Genetic Etiology Between Gastrointestinal Tract Diseases and Psychiatric Disorders: A Genome-Wide Pleiotropic Analysis," where PLACO analysis was conducted on trait pairs that exhibited high genetic correlation? Of course, they performed de-correlation of the Z-scores.
5.You mentioned in GitHub, "When samples are related, PLACO can use the summary statistics from EMMAX (or other univariate mixed model frameworks) to appropriately test for genetic associations." However, it seems that this step is not included in the PLACO R code. Is that the case, or is it that this step is not necessary?
6.You mentioned, "Harmonize the same effect allele across the two studies/traits so that Z-scores from the two datasets can be jointly analyzed appropriately using PLACO." For example, if I have traits A and B, and their effect alleles happen to be opposite, which trait's beta value should be flipped to be negative, and does the EAF need to be changed to 1-EAF? Or is it fine either way, just flipping the beta value of one trait? Is there anything else that needs to be flipped? 7.Because the computation time for PLACO is very lengthy, I believe the coding approach can be altered a bit. Since the range of the square of the Z score is from 0 to 80, we only need to first compute the VarZ for this pair of traits, then calculate each Z1Z2 (for calculating p.placo, only the absolute value of Z1Z2 is used, which ranges from 0 to 80, and I've set the interval at 0.001), and determine the corresponding p.placo value for this VarZ. This requires only 80,000 calculations, and the result will be the same as the original code (which requires computing the p.placo for each SNP in the entire GWAS, totaling millions of calculations). If you find any errors, please let me know. I have only changed the order of computation, not the calculation rules of PLACO.
Here is my code:
#First calculate VarZ, which is different for each pair of traits.
VarZ <- var.placo(Z.matrix.decor, P.matrix, p.threshold = 1e-04)
varz<-VarZ
print("Now let's start deriving the Zplus quantity.")
#Redefine the function, the meaning remains unchanged, z[1]*z[2] has been replaced with zplus
.pdfx<-function(x) besselK(x=abs(x),nu=0)/pi
.p.besselchange<-function(zplus, varz, AbsTol=1e-13){
p1<-2as.double(integrate(Vectorize(.pdfx),abs(zplus/sqrt(varz[1])),Inf,abs.tol=AbsTol)$value) p2<-2as.double(integrate(Vectorize(.pdfx),abs(zplus/sqrt(varz[2])),Inf,abs.tol=AbsTol)$value) p0<-2*as.double(integrate(Vectorize(.pdfx), abs(zplus),Inf, abs.tol=AbsTol)$value) pval.compnull<-p1+p2-p0 return(pval.compnull) }
zplus_values <- seq(0, 80, by=0.001)
pb <- txtProgressBar(min = 0, max = length(zplus_values), style = 3)
list_results <- lapply(seq_along(zplus_values), function(i) {
setTxtProgressBar(pb, i)
p_value <- .p.besselchange(zplus = zplus_values[i], varz = varz)
c(p.placo = p_value, Zplus = zplus_values[i])
})
close(pb)
results <- do.call(rbind, list_results)
results<-as.data.frame(results)
head(results)
Z.matrix.decor <- cbind(Z.matrix.decor, abs(Z.matrix.decor[,1] * Z.matrix.decor[,2]))
colnames(Z.matrix.decor)[3] <- "Zplus"
Z.matrix.decor[, "Zplus"] <- round(Z.matrix.decor[, "Zplus"], 3)
Z.matrix.decordataframe<-as.data.frame(Z.matrix.decor)
head(Z.matrix.decordataframe)
head(t)
combined_df <- cbind(t, Z.matrix.decordataframe)
head(combined_df)
combined_df$Zplus <- as.character(combined_df$Zplus)
#Convert to string type, otherwise there will be errors in matching floating-point numbers. results$Zplus <- as.character(results$Zplus)
out_2 <- merge(combined_df, results, by = "Zplus", all.x = TRUE)
out_2$T.placo<-out_2$Zplus
head(out_2)
I appreciate your time and assistance in answering these questions and am looking forward to implementing the PLACO method in my research.
6.You mentioned, "Harmonize the same effect allele across the two studies/traits so that Z-scores from the two datasets can be jointly analyzed appropriately using PLACO." For example, if I have traits A and B, and their effect alleles happen to be opposite, which trait's beta value should be flipped to be negative, and does the EAF need to be changed to 1-EAF? Or is it fine either way, just flipping the beta value of one trait? Is there anything else that needs to be flipped?
Hi @nvice111, when harmonizing two summary statistics I personally just flip one of the beta values to its negative, and EAF does not need to be flipped since EAF is not included in the matrix. Here is how I performed harmonization:
dat = merge(trait1, trait2, by=c("SNP", "chr"))
reversed_rows = (dat$A1.x == dat$A2.y) & (dat$A2.x == dat$A1.y)
dat$beta.y[reversed_rows] = -dat$beta.y[reversed_rows]
dat$A1.y[reversed_rows] = dat$A1.x[reversed_rows]
dat$A2.y[reversed_rows] = dat$A2.x[reversed_rows]
discard_rows = !(dat$A1.x %in% c(dat$A1.y, dat$A2.y) & dat$A2.x %in% c(dat$A1.y, dat$A2.y))
dat = dat[!discard_rows, ]
I look forward to your implementation.
