Last updated: 2018-10-02
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: f66679b
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: 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/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/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/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/PeakToGeneAssignment.Rmd
Modified: analysis/callMolQTLS.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/test.max2.Rmd
Modified: code/Snakefile
Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.
File | Version | Author | Date | Message |
---|---|---|---|---|
Rmd | f66679b | Briana Mittleman | 2018-10-02 | overlap plots at peak level |
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':
smiths
library(workflowr)
This is workflowr version 1.1.1
Run ?workflowr for help getting started
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= 677"
[1] "Nuclear overlap= 250"
[1] "Nuclear overlap= 147"
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= 3049"
[1] "Total overlap= 689"
[1] "Total overlap= 561"
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")) %>% 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")) %>% 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){
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)
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= 3049"
[1] "Total overlap= 2352"
[1] "Total overlap= 2369"
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")) %>% 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")) %>% 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){
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)
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= 677"
[1] "Nuclear overlap= 574"
[1] "Nuclear overlap= 527"
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 workflowr_1.1.1 reshape2_1.4.3 forcats_0.3.0
[5] stringr_1.3.1 dplyr_0.7.6 purrr_0.2.5 readr_1.1.1
[9] tidyr_0.8.1 tibble_1.4.2 ggplot2_3.0.0 tidyverse_1.2.1
loaded via a namespace (and not attached):
[1] tidyselect_0.2.4 haven_1.1.2 lattice_0.20-35
[4] colorspace_1.3-2 htmltools_0.3.6 yaml_2.2.0
[7] rlang_0.2.2 R.oo_1.22.0 pillar_1.3.0
[10] glue_1.3.0 withr_2.1.2 R.utils_2.7.0
[13] modelr_0.1.2 readxl_1.1.0 bindr_0.1.1
[16] plyr_1.8.4 munsell_0.5.0 gtable_0.2.0
[19] cellranger_1.1.0 rvest_0.3.2 R.methodsS3_1.7.1
[22] evaluate_0.11 knitr_1.20 broom_0.5.0
[25] Rcpp_0.12.18 scales_1.0.0 backports_1.1.2
[28] jsonlite_1.5 hms_0.4.2 digest_0.6.16
[31] stringi_1.2.4 grid_3.5.1 rprojroot_1.3-2
[34] cli_1.0.0 tools_3.5.1 magrittr_1.5
[37] lazyeval_0.2.1 crayon_1.3.4 whisker_0.3-2
[40] pkgconfig_2.0.2 xml2_1.2.0 lubridate_1.7.4
[43] assertthat_0.2.0 rmarkdown_1.10 httr_1.3.1
[46] rstudioapi_0.7 R6_2.2.2 nlme_3.1-137
[49] git2r_0.23.0 compiler_3.5.1
This reproducible R Markdown analysis was created with workflowr 1.1.1