diff --git a/R/qc_report_loci.R b/R/qc_report_loci.R index 9af6213e..d87d5590 100644 --- 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") }