Merge pull request #42 from EvolEcolGroup/qc_fix

Edit qc_report loci plots
EvolEcolGroup · Jun 11, 2024 · e9c3ce2 · e9c3ce2
2 parents f6747f3 + bf3cfae
commit e9c3ce2
Showing 1 changed file with 10 additions and 7 deletions.
diff --git a/R/qc_report_loci.R b/R/qc_report_loci.R
@@ -78,7 +78,7 @@ autoplot.qc_report_loci <- function(object,
   logp <- -log10(hwe_p)
 
   if (type == "overview") {
-    final_plot <- autoplot_l_qc_overview(object, maf_threshold, miss_threshold)
+    final_plot <- autoplot_l_qc_overview(object, maf_threshold, miss_threshold, hwe_p)
   } else if (type == "all") {
     final_plot <- autoplot_l_qc_all(object, maf_threshold, miss_threshold, hwe_p,logp)
   } else if (type == "missing") {
@@ -125,9 +125,12 @@ autoplot_l_qc_all <- function(object, maf_threshold = maf_threshold, miss_thresh
 
   #Hardy weinberg exact test p-val distribution
   qc_report$hwe_p_log <- -log10(qc_report$hwe_p)
-  qc_lowhwe <- subset(qc_report,qc_report$hwe_p < hwe_p)
+  #qc_lowhwe <- subset(qc_report,qc_report$hwe_p < hwe_p)
+  qc_lowhwe <- qc_report[qc_report$hwe_p < hwe_p, ]
+  #qc_lowhwe$hwe_p_log <- -log10(qc_lowhwe$hwe_p)
 
   hwe_all <- ggplot2::ggplot(qc_report,ggplot2::aes(x=.data$hwe_p_log))+ggplot2::geom_histogram(binwidth = 0.5,fill="#66C2A5")+ ggplot2::labs(x=expression("-log"[10]*" of HWE exact p-value"),y="Number of SNPs", title = "Hardy-Weinberg exact test")+ ggplot2::geom_vline(xintercept= logp, lty=2, col="red")
+
   hwe_low <- ggplot2::ggplot(qc_lowhwe,ggplot2::aes(x=.data$hwe_p_log))+ggplot2::geom_histogram(binwidth = 0.5,fill="#66C2A5")+ ggplot2::labs(x=expression("-log"[10]* " of HWE exact p-value"),y="Number of SNPs", title = paste("HWE exact p-value <", hwe_p))+ ggplot2::geom_vline(xintercept= logp, lty=2, col="red")
 
   hwes <- patchwork::wrap_plots(hwe_all,hwe_low)
@@ -137,18 +140,18 @@ autoplot_l_qc_all <- function(object, maf_threshold = maf_threshold, miss_thresh
 }
 
 
-autoplot_l_qc_overview <- function(object,maf_threshold=maf_threshold, miss_threshold = miss_threshold,...){
+autoplot_l_qc_overview <- function(object,hwe_p = hwe_p, maf_threshold=maf_threshold, miss_threshold = miss_threshold,...){
 
   qc_report <- object
 
-  qc_lowmaf <- subset(qc_report, qc_report$maf <= maf_threshold)
-  qc_lowhwe <- subset(qc_report,qc_report$hwe_p < 0.01)
+  qc_lowmaf <- subset(qc_report, qc_report$maf >= maf_threshold)
+  qc_lowhwe <- subset(qc_report,qc_report$hwe_p >= hwe_p)
 
 
   maf_fail <- c(qc_lowmaf$snp_id)
   hwe_fail <- c(qc_lowhwe$snp_id)
 
-  qc_highmissing <- subset(qc_report,qc_report$missingness>=0.05)
+  qc_highmissing <- subset(qc_report,qc_report$missingness <= miss_threshold)
   missing_fail <- c(qc_highmissing$snp_id)
 
   list <- list(MAF = maf_fail, HWE = hwe_fail, Missing = missing_fail)
@@ -183,7 +186,7 @@ autoplot_l_qc_sig_hwe <- function(object,hwe_p=hwe_p,logp,...){
 
   #Hardy weinberg exact test p-val distribution
   qc_report$hwe_p_log <- -log10(qc_report$hwe_p)
-  qc_lowhwe <- subset(qc_report,qc_report$hwe_p < hwe_p)
+  qc_lowhwe <- subset(qc_report,qc_report$hwe_p < logp)
 
   hwe_low <- ggplot2::ggplot(qc_lowhwe,ggplot2::aes(x=.data$hwe_p_log))+ggplot2::geom_histogram(binwidth = 0.5,fill="#66C2A5")+ ggplot2::labs(x=expression("-log"[10]* " of HWE exact p-value"),y="Number of SNPs", title = paste("HWE exact p-value <", hwe_p))+ ggplot2::geom_vline(xintercept= logp, lty=2, col="red")
 }