Condition on Types ApaQTLs

Make QQplots

overlapQTLplot_totalQTL=function(cut, plotfile){

    #helper functions  
    sigsnp=function(cutoff){
      permTot$bh=p.adjust(permTot$bpval, method="fdr")
      file_sig=permTot %>% filter(-log10(bh)> cutoff) %>% select(sid)
      print(paste("Sig snps=", nrow(file_sig), sep=" "))
      return(file_sig)
    }
    randomsnps=function(SigSnpList){
      nsnp=nrow(SigSnpList)
      randomSnpDF= permTot %>% sample_n(nsnp) %>% arrange(sid) %>% select(sid)
      return(randomSnpDF)
    }
    top_Nuclear=function(snp_list){
      filt_nuc=permNuc %>% semi_join(snp_list, by="sid") %>% group_by(sid) %>% add_tally() %>% ungroup() %>% mutate(corrPval=bpval)
      filt_nuc_top= filt_nuc %>% group_by(sid) %>% top_n(-1, corrPval)
      print(paste("Nuclear overlap=", nrow(filt_nuc_top), sep=" "))
      return(filt_nuc_top)
    }
    makeQQ=function(test, baseline){
      plot=qqplot(-log10(runif(nrow(baseline))), -log10(baseline$corrPval), ylab="Observed", xlab="Expected", main="Significant Total QTLs- nuclear Pval")
      points(sort(-log10(runif(nrow(test)))), sort(-log10(test$corrPval)), col= alpha("Red"))
      abline(0,1)
      return(plot)
    }
  
    TL=sigsnp(cut)
    BL=randomsnps(TL)
    #top snps test and base total
    topN_T=top_Nuclear(TL)
    topN_B=top_Nuclear(BL)
    #plot Total
    png(plotfile)
    totalPlot=makeQQ(topN_T,topN_B)
    dev.off()

}

overlapQTLplot_totalQTL(1, "../output/plots/TotalQTLinNuclear.png")

[1] "Sig snps= 118"
[1] "Nuclear overlap= 31"
[1] "Nuclear overlap= 9"

quartz_off_screen 
                2 

overlapQTLplot_totalQTL=function(cut, plotfile){

    #helper functions  
    sigsnp=function(cutoff){
      permNuc$bh=p.adjust(permNuc$bpval, method="fdr")
      file_sig=permNuc %>% filter(-log10(bh)> cutoff) %>% select(sid)
      print(paste("Sig snps=", nrow(file_sig), sep=" "))
      return(file_sig)
    }
    randomsnps=function(SigSnpList){
      nsnp=nrow(SigSnpList)
      randomSnpDF= permNuc %>% sample_n(nsnp) %>% arrange(sid) %>% select(sid)
      return(randomSnpDF)
    }
    top_Total=function(snp_list){
      filt_tot=permTot %>% semi_join(snp_list, by="sid") %>% group_by(sid) %>% add_tally() %>% ungroup() %>% mutate(corrPval=bpval)
      filt_tot_top= filt_tot %>% group_by(sid) %>% top_n(-1, corrPval)
      print(paste("Total overlap=", nrow(filt_tot_top), sep=" "))
      return(filt_tot_top)
    }
    makeQQ=function(test, baseline){
      plot=qqplot(-log10(runif(nrow(baseline))), -log10(baseline$corrPval), ylab="Observed", xlab="Expected", main="Significant Nuclear QTLs- Total Pval")
      points(sort(-log10(runif(nrow(test)))), sort(-log10(test$corrPval)), col= alpha("Red"))
      abline(0,1)
      return(plot)
    }
  
    TL=sigsnp(cut)
    BL=randomsnps(TL)
    #top snps test and base total
    topN_T=top_Total(TL)
    topN_B=top_Total(BL)
    #plot Total
    png(plotfile)
    totalPlot=makeQQ(topN_T,topN_B)
    dev.off()

}
overlapQTLplot_totalQTL(1, "../output/plots/NuclearQTLinTotal.png")

[1] "Sig snps= 880"
[1] "Total overlap= 85"
[1] "Total overlap= 48"

quartz_off_screen 
                2 

I should change this to focus on peak. I can say give me the genes with significant QTLs in total or nuclear then look at the pvalues for those peaks in the other file. As I did before, I am going to work on all of the functions seperatly then put them together.

Get the peaks with significant QTLs, and the same number of random peaks.

sigpeak=function(cutoff){
    permNuc$bh=p.adjust(permNuc$bpval, method="fdr")
    file_sig=permNuc %>% filter(-log10(bh)> cutoff) %>%separate(pid, sep = ":", into=c("chr", "start", "end", "id")) %>% separate(id, sep = "_", into=c("gene", "strand", "peak")) %>% select(peak)
    print(paste("Sig peaks=", nrow(file_sig), sep=" "))
    return(file_sig)
    }
    
x=sigpeak(1)

randompeak=function(SigSnpList){
  nsnp=nrow(SigSnpList)
  randomPeakDF= permNuc %>% sample_n(nsnp) %>% separate(pid, sep = ":", into=c("chr", "start", "end", "id")) %>% separate(id, sep = "_", into=c("gene", "strand", "peak")) %>% select(peak)
  return(randomPeakDF)
}

y=randompeak(x)

I can now get the top pvalue for each of these using the total permuted pval.

Peak_overlap=function(snp_list){
      filt_tot=permTot %>% separate(pid, sep = ":", into=c("chr", "start", "end", "id")) %>% separate(id, sep = "_", into=c("gene", "strand", "peak")) %>% semi_join(snp_list, by="peak")
      print(paste("Total overlap=", nrow(filt_tot), sep=" "))
      return(filt_tot)
}


#run on real sig peaks and random peaks  
Test=Peak_overlap(x)
base=Peak_overlap(y)

Plot:

  makeQQ_peak=function(test, baseline){
      plot=qqplot(-log10(runif(nrow(test))), -log10(test$bpval), ylab="Observed", xlab="Expected", main="Peaks with Significant QTLs in Nuc \n pvalues in Tot")
      points(sort(-log10(runif(nrow(baseline)))), sort(-log10(baseline$bpval)), col=alpha("Red"))
      abline(0,1)
      return(plot)
}
plot=makeQQ_peak(Test,base)

SigNucPeakOverlapTot=function(cutoff, plotfile){
    sigpeak=function(cutoff){
      permNuc$bh=p.adjust(permNuc$bpval, method="fdr")
      file_sig=permNuc %>% filter(-log10(bh)> cutoff) %>%separate(pid, sep = ":", into=c("chr", "start", "end", "id")) %>% separate(id, sep = "_", into=c("gene", "strand", "peak")) %>% dplyr::select(peak)
      print(paste("Sig peaks=", nrow(file_sig), sep=" "))
      return(file_sig)
      }
  
  randompeak=function(SigSnpList){
    nsnp=nrow(SigSnpList)
    randomPeakDF= permNuc %>% sample_n(nsnp) %>% separate(pid, sep = ":", into=c("chr", "start", "end", "id")) %>% separate(id, sep = "_", into=c("gene", "strand", "peak")) %>% dplyr::select(peak)
    return(randomPeakDF)
  }
  Peak_overlap=function(snp_list){
        filt_tot=permTot %>% separate(pid, sep = ":", into=c("chr", "start", "end", "id")) %>% separate(id, sep = "_", into=c("gene", "strand", "peak")) %>% semi_join(snp_list, by="peak")
        print(paste("Total overlap=", nrow(filt_tot), sep=" "))
        return(filt_tot)
  }

  makeQQ_peak=function(test, baseline){
      p0test=pi0est(test$bpval)
      p1test=1-p0test$pi0
      plot=qqplot(-log10(runif(nrow(baseline))), -log10(baseline$bpval), ylab="Observed", xlab="Expected", main="Peaks with Significant QTLs in Nuclear \n pvalues in Total")
      points(sort(-log10(runif(nrow(test)))), sort(-log10(test$bpval)), col= alpha("Red"))
      abline(0,1)
      text(1.5,3, paste("pi_1=", round(p1test, digit=3), sep=" "))
      return(plot)
  }
  
  
testPeaks=sigpeak(1)
basePeaks=randompeak(testPeaks)
testSet=Peak_overlap(testPeaks)
baselineSet=Peak_overlap(basePeaks)
png(plotfile)
plot=makeQQ_peak(testSet,baselineSet )
dev.off()
}

SigNucPeakOverlapTot(1, "../output/plots/SigNucPeakTotpval.png")

[1] "Sig peaks= 880"
[1] "Total overlap= 353"
[1] "Total overlap= 348"

quartz_off_screen 
                2 

SigTotPeakOverlapNuc=function(cutoff, plotfile){
    sigpeak=function(cutoff){
      permTot$bh=p.adjust(permTot$bpval, method="fdr")
      file_sig=permTot %>% filter(-log10(bh)> cutoff) %>% separate(pid, sep = ":", into=c("chr", "start", "end", "id")) %>% separate(id, sep = "_", into=c("gene", "strand", "peak")) %>% dplyr::select(peak)
      print(paste("Sig peaks=", nrow(file_sig), sep=" "))
      return(file_sig)
      }
  
  randompeak=function(SigSnpList){
    nsnp=nrow(SigSnpList)
    randomPeakDF= permTot %>% sample_n(nsnp) %>% separate(pid, sep = ":", into=c("chr", "start", "end", "id")) %>% separate(id, sep = "_", into=c("gene", "strand", "peak")) %>% dplyr::select(peak)
    return(randomPeakDF)
  }
  Peak_overlap=function(snp_list){
        filt_nuc=permNuc %>% separate(pid, sep = ":", into=c("chr", "start", "end", "id")) %>% separate(id, sep = "_", into=c("gene", "strand", "peak")) %>% semi_join(snp_list, by="peak")
        print(paste("Nuclear overlap=", nrow(filt_nuc), sep=" "))
        return(filt_nuc)
  }

  makeQQ_peak=function(test, baseline){
      p0test=pi0est(test$bpval)
      p1test=1-p0test$pi0
      plot=qqplot(-log10(runif(nrow(baseline))), -log10(baseline$bpval), ylab="Observed", xlab="Expected", main="Peaks with Significant QTLs in Total \n pvalues in Nuclear")
      points(sort(-log10(runif(nrow(test)))), sort(-log10(test$bpval)), col= alpha("Red"))
      abline(0,1)
      text(2,3, paste("pi_1=", round(p1test, digit=3), sep=" "))
      return(plot)
  }
  
  
testPeaks=sigpeak(1)
basePeaks=randompeak(testPeaks)
testSet=Peak_overlap(testPeaks)
baselineSet=Peak_overlap(basePeaks)
png(plotfile)
plot=makeQQ_peak(testSet,baselineSet )
dev.off()
}

SigTotPeakOverlapNuc(1, "../output/plots/SigTotPeakNucpval.png")

[1] "Sig peaks= 118"
[1] "Nuclear overlap= 69"
[1] "Nuclear overlap= 70"

quartz_off_screen 
                2 

Try historgram:

SigNucPeakOverlapTot_hist=function(cutoff, plotfile){
    sigpeak=function(cutoff){
      permNuc$bh=p.adjust(permNuc$bpval, method="fdr")
      file_sig=permNuc %>% filter(-log10(bh)> cutoff) %>%separate(pid, sep = ":", into=c("chr", "start", "end", "id")) %>% separate(id, sep = "_", into=c("gene", "strand", "peak")) %>% dplyr::select(peak)
      print(paste("Sig peaks=", nrow(file_sig), sep=" "))
      return(file_sig)
      }
  
  randompeak=function(SigSnpList){
    nsnp=nrow(SigSnpList)
    randomPeakDF= permNuc %>% sample_n(nsnp) %>% separate(pid, sep = ":", into=c("chr", "start", "end", "id")) %>% separate(id, sep = "_", into=c("gene", "strand", "peak")) %>% dplyr::select(peak)
    return(randomPeakDF)
  }
  Peak_overlap=function(snp_list){
        filt_tot=permTot %>% separate(pid, sep = ":", into=c("chr", "start", "end", "id")) %>% separate(id, sep = "_", into=c("gene", "strand", "peak")) %>% semi_join(snp_list, by="peak")
        print(paste("Total overlap=", nrow(filt_tot), sep=" "))
        return(filt_tot)
  }

  makeQQ_peak=function(test, baseline){
      p0test=pi0est(test$bpval)
      p1test=1-p0test$pi0
      plot=hist(test$bpval, breaks=20, main="Peaks with Significant QTLs in Nuclear \n Pvalues in Total", xlab="Total APAqtl Pvalue")
      text(.8,140, paste("pi_1=", round(p1test, digit=3), sep=" "))
      return(plot)
  }
  
  
testPeaks=sigpeak(1)
basePeaks=randompeak(testPeaks)
testSet=Peak_overlap(testPeaks)
baselineSet=Peak_overlap(basePeaks)
png(plotfile)
plot=makeQQ_peak(testSet,baselineSet )
dev.off()
}

SigNucPeakOverlapTot_hist(1, "../output/plots/SigNucPeakTotpval_hist.png")

[1] "Sig peaks= 880"
[1] "Total overlap= 353"
[1] "Total overlap= 352"

quartz_off_screen 
                2 

SigTotPeakOverlapNuc_hist=function(cutoff, plotfile){
    sigpeak=function(cutoff){
      permTot$bh=p.adjust(permTot$bpval, method="fdr")
      file_sig=permTot %>% filter(-log10(bh)> cutoff) %>% separate(pid, sep = ":", into=c("chr", "start", "end", "id")) %>% separate(id, sep = "_", into=c("gene", "strand", "peak")) %>% dplyr::select(peak)
      print(paste("Sig peaks=", nrow(file_sig), sep=" "))
      return(file_sig)
      }
  
  randompeak=function(SigSnpList){
    nsnp=nrow(SigSnpList)
    randomPeakDF= permTot %>% sample_n(nsnp) %>% separate(pid, sep = ":", into=c("chr", "start", "end", "id")) %>% separate(id, sep = "_", into=c("gene", "strand", "peak")) %>% dplyr::select(peak)
    return(randomPeakDF)
  }
  Peak_overlap=function(snp_list){
        filt_nuc=permNuc %>% separate(pid, sep = ":", into=c("chr", "start", "end", "id")) %>% separate(id, sep = "_", into=c("gene", "strand", "peak")) %>% semi_join(snp_list, by="peak")
        print(paste("Nuclear overlap=", nrow(filt_nuc), sep=" "))
        return(filt_nuc)
  }

  makeQQ_peak=function(test, baseline){
      p0test=pi0est(test$bpval)
      p1test=1-p0test$pi0
      plot=hist(test$bpval, breaks=20, main="Peaks with Significant QTLs in Total \n Pvalues in Nuclear",xlab="Nuclear APAqtl Pvalue")
      text(.8,40, paste("pi_1=", round(p1test, digit=3), sep=" "))
      return(plot)
  }
  
  
testPeaks=sigpeak(1)
basePeaks=randompeak(testPeaks)
testSet=Peak_overlap(testPeaks)
baselineSet=Peak_overlap(basePeaks)
png(plotfile)
plot=makeQQ_peak(testSet,baselineSet )
dev.off()
}

SigTotPeakOverlapNuc_hist(1, "../output/plots/SigTotPeakNucpval_hist.png")

[1] "Sig peaks= 118"
[1] "Nuclear overlap= 69"
[1] "Nuclear overlap= 60"

quartz_off_screen 
                2 

Session information

sessionInfo()

R version 3.5.1 (2018-07-02)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS  10.14.1

Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] bindrcpp_0.2.2  qvalue_2.12.0   workflowr_1.2.0 reshape2_1.4.3 
 [5] forcats_0.3.0   stringr_1.4.0   dplyr_0.7.6     purrr_0.2.5    
 [9] readr_1.1.1     tidyr_0.8.1     tibble_1.4.2    ggplot2_3.0.0  
[13] tidyverse_1.2.1

loaded via a namespace (and not attached):
 [1] tidyselect_0.2.4 splines_3.5.1    haven_1.1.2      lattice_0.20-35 
 [5] colorspace_1.3-2 htmltools_0.3.6  yaml_2.2.0       rlang_0.2.2     
 [9] pillar_1.3.0     glue_1.3.0       withr_2.1.2      modelr_0.1.2    
[13] readxl_1.1.0     bindr_0.1.1      plyr_1.8.4       munsell_0.5.0   
[17] gtable_0.2.0     cellranger_1.1.0 rvest_0.3.2      evaluate_0.13   
[21] knitr_1.20       broom_0.5.0      Rcpp_0.12.19     scales_1.0.0    
[25] backports_1.1.2  jsonlite_1.6     fs_1.2.6         hms_0.4.2       
[29] digest_0.6.17    stringi_1.2.4    grid_3.5.1       rprojroot_1.3-2 
[33] cli_1.0.1        tools_3.5.1      magrittr_1.5     lazyeval_0.2.1  
[37] crayon_1.3.4     whisker_0.3-2    pkgconfig_2.0.2  xml2_1.2.0      
[41] lubridate_1.7.4  assertthat_0.2.0 rmarkdown_1.11   httr_1.3.1      
[45] rstudioapi_0.9.0 R6_2.3.0         nlme_3.1-137     git2r_0.24.0    
[49] compiler_3.5.1