Last updated: 2018-10-08
workflowr checks: (Click a bullet for more information) ✔ R Markdown file: up-to-date
Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.
✔ Environment: empty
Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.
✔ Seed:
set.seed(12345)
The command set.seed(12345) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.
✔ Session information: recorded
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
✔ Repository version: 0c2fffb
wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:
Ignored files:
Ignored: .DS_Store
Ignored: .Rhistory
Ignored: .Rproj.user/
Ignored: output/.DS_Store
Untracked files:
Untracked: KalistoAbundance18486.txt
Untracked: analysis/genometrack_figs.Rmd
Untracked: analysis/ncbiRefSeq_sm.sort.mRNA.bed
Untracked: analysis/snake.config.notes.Rmd
Untracked: analysis/verifyBAM.Rmd
Untracked: data/18486.genecov.txt
Untracked: data/APApeaksYL.total.inbrain.bed
Untracked: data/NuclearApaQTLs.txt
Untracked: data/RNAkalisto/
Untracked: data/TotalApaQTLs.txt
Untracked: data/Totalpeaks_filtered_clean.bed
Untracked: data/YL-SP-18486-T-combined-genecov.txt
Untracked: data/YL-SP-18486-T_S9_R1_001-genecov.txt
Untracked: data/apaExamp/
Untracked: data/bedgraph_peaks/
Untracked: data/bin200.5.T.nuccov.bed
Untracked: data/bin200.Anuccov.bed
Untracked: data/bin200.nuccov.bed
Untracked: data/clean_peaks/
Untracked: data/comb_map_stats.csv
Untracked: data/comb_map_stats.xlsx
Untracked: data/comb_map_stats_39ind.csv
Untracked: data/combined_reads_mapped_three_prime_seq.csv
Untracked: data/ensemble_to_genename.txt
Untracked: data/filtered_APApeaks_merged_allchrom_refseqTrans.closest2End.bed
Untracked: data/filtered_APApeaks_merged_allchrom_refseqTrans.closest2End.noties.bed
Untracked: data/first50lines_closest.txt
Untracked: data/gencov.test.csv
Untracked: data/gencov.test.txt
Untracked: data/gencov_zero.test.csv
Untracked: data/gencov_zero.test.txt
Untracked: data/gene_cov/
Untracked: data/joined
Untracked: data/leafcutter/
Untracked: data/merged_combined_YL-SP-threeprimeseq.bg
Untracked: data/mol_overlap/
Untracked: data/nom_QTL/
Untracked: data/nom_QTL_opp/
Untracked: data/nom_QTL_trans/
Untracked: data/nuc6up/
Untracked: data/other_qtls/
Untracked: data/peakPerRefSeqGene/
Untracked: data/perm_QTL/
Untracked: data/perm_QTL_opp/
Untracked: data/perm_QTL_trans/
Untracked: data/reads_mapped_three_prime_seq.csv
Untracked: data/smash.cov.results.bed
Untracked: data/smash.cov.results.csv
Untracked: data/smash.cov.results.txt
Untracked: data/smash_testregion/
Untracked: data/ssFC200.cov.bed
Untracked: data/temp.file1
Untracked: data/temp.file2
Untracked: data/temp.gencov.test.txt
Untracked: data/temp.gencov_zero.test.txt
Untracked: output/picard/
Untracked: output/plots/
Untracked: output/qual.fig2.pdf
Unstaged changes:
Modified: analysis/28ind.peak.explore.Rmd
Modified: analysis/39indQC.Rmd
Modified: analysis/PeakToGeneAssignment.Rmd
Modified: analysis/cleanupdtseq.internalpriming.Rmd
Modified: analysis/dif.iso.usage.leafcutter.Rmd
Modified: analysis/diff_iso_pipeline.Rmd
Modified: analysis/explore.filters.Rmd
Modified: analysis/overlapMolQTL.Rmd
Modified: analysis/overlap_qtls.Rmd
Modified: analysis/peakOverlap_oppstrand.Rmd
Modified: analysis/pheno.leaf.comb.Rmd
Modified: analysis/swarmPlots_QTLs.Rmd
Modified: analysis/test.max2.Rmd
Modified: code/Snakefile
In this analysis I want to use the resuls from the total and nuclear APA qtl calling. I will ask if conditioning on a nuclear QTL increases the signal in the total QTL and vice versa. I will start with the significant snp-peak pairs from the permuted files. I will then overlap with the nominal pvalues from the other fraction. I will do this similar to how I did in the overlaMolQTL analysis. However in this analysis I do not have the multiple peaks per gene problem that I have when I overlap the pvalues. I can map the same peak to snp pair.
Due to file size I will do this only with the permuted files.
library(tidyverse)── Attaching packages ──────────────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──✔ ggplot2 3.0.0 ✔ purrr 0.2.5
✔ tibble 1.4.2 ✔ dplyr 0.7.6
✔ tidyr 0.8.1 ✔ stringr 1.3.1
✔ readr 1.1.1 ✔ forcats 0.3.0── Conflicts ─────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()library(reshape2)
Attaching package: 'reshape2'The following object is masked from 'package:tidyr':
smithslibrary(workflowr)This is workflowr version 1.1.1
Run ?workflowr for help getting started library(qvalue)Permuted
permTot=read.table("../data/perm_QTL_trans/filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Total_transcript_permResBH.txt",header=T ,stringsAsFactors = F)
permNuc=read.table("../data/perm_QTL_trans/filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Nuclear_transcript_permResBH.txt",header=T, stringsAsFactors = F)Nominal
nomnames=c("pid", "sid", "dist", "npval", "slope")
#nomTot=read_table("../data/nom_QTL_trans/filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Total_NomRes.txt", col_names=nomnames, col_types = c(col_character(), col_character(), col_double(), col_double(), col_double()))
#nomNuc=read_table("../data/nom_QTL_trans/filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Nuclear_NomRes.txt", col_names=nomnames, col_types = c(col_character(), col_character(), col_double(), col_double(), col_double()))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 sessionInfo()R version 3.5.1 (2018-07-02)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Sierra 10.12.6
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.1.1 reshape2_1.4.3
[5] forcats_0.3.0 stringr_1.3.1 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
[4] lattice_0.20-35 colorspace_1.3-2 htmltools_0.3.6
[7] yaml_2.2.0 rlang_0.2.2 R.oo_1.22.0
[10] pillar_1.3.0 glue_1.3.0 withr_2.1.2
[13] R.utils_2.7.0 modelr_0.1.2 readxl_1.1.0
[16] bindr_0.1.1 plyr_1.8.4 munsell_0.5.0
[19] gtable_0.2.0 cellranger_1.1.0 rvest_0.3.2
[22] R.methodsS3_1.7.1 evaluate_0.11 knitr_1.20
[25] broom_0.5.0 Rcpp_0.12.19 scales_1.0.0
[28] backports_1.1.2 jsonlite_1.5 hms_0.4.2
[31] digest_0.6.17 stringi_1.2.4 grid_3.5.1
[34] rprojroot_1.3-2 cli_1.0.1 tools_3.5.1
[37] magrittr_1.5 lazyeval_0.2.1 crayon_1.3.4
[40] whisker_0.3-2 pkgconfig_2.0.2 xml2_1.2.0
[43] lubridate_1.7.4 assertthat_0.2.0 rmarkdown_1.10
[46] httr_1.3.1 rstudioapi_0.8 R6_2.3.0
[49] nlme_3.1-137 git2r_0.23.0 compiler_3.5.1
This reproducible R Markdown analysis was created with workflowr 1.1.1