Last updated: 2018-11-06
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: b5f744f
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: data/.DS_Store
Ignored: output/.DS_Store
Untracked files:
Untracked: KalistoAbundance18486.txt
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/ChromHmmOverlap/
Untracked: data/GM12878.chromHMM.bed
Untracked: data/GM12878.chromHMM.txt
Untracked: data/NuclearApaQTLs.txt
Untracked: data/PeaksUsed/
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/diff_iso_trans/
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/mol_pheno/
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/apaQTLoverlapGWAS.Rmd
Modified: analysis/characterizeNuclearApaQtls.Rmd
Modified: analysis/cleanupdtseq.internalpriming.Rmd
Modified: analysis/coloc_apaQTLs_protQTLs.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
| File | Version | Author | Date | Message |
|---|---|---|---|---|
| Rmd | b5f744f | Briana Mittleman | 2018-11-06 | initiate flash2mash |
I will use this analysis to implement the flash and mash packages developed by the stephens lab to better understand molecular QTL sharing and to see if adding APA to a model can help with power in protein QTLs.
Steps: 1. FLASH to see tissue patterns (https://willwerscheid.github.io/MASHvFLASH/MASHvFLASHnn.html and https://willwerscheid.github.io/MASHvFLASH/MASHvFLASHnn2.html)
2. Conditional analysis with residuals to see if I can call APA qtls on the residuals from an RNA~protein analysis 3. run MASH
Data stucture: I need to have a matrix with all of my QTL results. I want to get a snp-gene by phenotype matrix with the effect sizes and standard errors. First I will do this with the genes we have all data for (unless it is too small). To deal with the APA isoform problem I will use the peak with the most significant peak-snp pair. This should be ok because given the peaks are ratios they are all correlated with eachother.
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(workflowr)This is workflowr version 1.1.1
Run ?workflowr for help getting startedFirst I can use the permuted results to look at the genes that are tested in all of the phenotypes.
read_permfile=function(file, mol){
perm_names=c("pid" ,"nvar","shape1" ,"shape2", "dummy","sid" ,"dist","npval", "slope" , "ppval" ,"bpval")
geneNames=read.table("../data/ensemble_to_genename.txt", sep="\t", header=T,stringsAsFactors = F)
res=read.table(file, col.names = perm_names, stringsAsFactors = F)
if (mol == "protein"){
res_f= res %>% rename("Gene.stable.ID"=pid)
res_final= res_f %>% inner_join(geneNames, by="Gene.stable.ID") %>% select(c("Gene.name"))
}
else{
res_final =res %>% separate(pid, into=c("Gene.stable.ID", "ver"), sep ="[.]") %>% inner_join(geneNames, by="Gene.stable.ID") %>% select(c("Gene.name"))
}
return(res_final)
}prot_res=read_permfile("../data/other_qtls/fastqtl_qqnorm_prot.fixed.perm.out", "protein")
rna_res=read_permfile("../data/other_qtls/fastqtl_qqnorm_RNAseq_phase2.fixed.perm.out", "RNA")
rnaG_res=read_permfile("../data/other_qtls/fastqtl_qqnorm_RNAseqGeuvadis.fixed.perm.out", "RNAG")
su30_res=read_permfile("../data/other_qtls/fastqtl_qqnorm_4su30.fixed.perm.out", "su30")
su60_res=read_permfile("../data/other_qtls/fastqtl_qqnorm_4su60.fixed.perm.out", "su60")
ribo_res=read_permfile("../data/other_qtls/fastqtl_qqnorm_ribo_phase2.fixed.perm.out", "ribo")Now I need to look at the apa file genes.
NuclearAPA=read.table("../data/perm_QTL_trans/filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Nuclear_transcript_permResBH.txt", stringsAsFactors = F, header=T) %>% separate(pid, sep = ":", into=c("chr", "start", "end", "id")) %>% separate(id, sep = "_", into=c("gene", "strand", "peak")) %>% rename("Gene.name"=gene) %>% select(Gene.name)%>% distinct()
totalAPA=read.table("../data/perm_QTL_trans/filtered_APApeaks_merged_allchrom_refseqGenes_pheno_Total_transcript_permResBH.txt", stringsAsFactors = F, header=T) %>% separate(pid, sep = ":", into=c("chr", "start", "end", "id")) %>% separate(id, sep = "_", into=c("gene", "strand", "peak")) %>% rename("Gene.name"=gene) %>% select(Gene.name) %>% distinct()Look hoqw many genes are in all sets:
allgenes= NuclearAPA %>% inner_join(totalAPA,by="Gene.name") %>% inner_join(totalAPA,by="Gene.name") %>% inner_join(su30_res,by="Gene.name") %>% inner_join(su60_res,by="Gene.name") %>% inner_join(rna_res,by="Gene.name") %>% inner_join(rnaG_res,by="Gene.name")%>% inner_join(ribo_res,by="Gene.name")%>% inner_join(prot_res,by="Gene.name")
print(nrow(allgenes))[1] 904allgenes_minusprot= NuclearAPA %>% inner_join(totalAPA,by="Gene.name") %>% inner_join(totalAPA,by="Gene.name") %>% inner_join(su30_res,by="Gene.name") %>% inner_join(su60_res,by="Gene.name") %>% inner_join(rna_res,by="Gene.name") %>% inner_join(rnaG_res,by="Gene.name")%>% inner_join(ribo_res,by="Gene.name")
print(nrow(allgenes_minusprot))[1] 2195allgenes_minusribo= NuclearAPA %>% inner_join(totalAPA,by="Gene.name") %>% inner_join(totalAPA,by="Gene.name") %>% inner_join(su30_res,by="Gene.name") %>% inner_join(su60_res,by="Gene.name") %>% inner_join(rna_res,by="Gene.name") %>% inner_join(rnaG_res,by="Gene.name")%>% inner_join(prot_res,by="Gene.name")
print(nrow(allgenes_minusribo))[1] 904genes_ApaRnaProt= NuclearAPA %>% inner_join(totalAPA,by="Gene.name") %>%inner_join(rna_res,by="Gene.name") %>%inner_join(prot_res,by="Gene.name")
print(nrow(genes_ApaRnaProt))[1] 904genes_RNAProt= rna_res%>%inner_join(prot_res,by="Gene.name")
print(nrow(genes_RNAProt))[1] 4131Only have 904 genes that are tested in both APA and protein data.
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] workflowr_1.1.1 forcats_0.3.0 stringr_1.3.1 dplyr_0.7.6
[5] purrr_0.2.5 readr_1.1.1 tidyr_0.8.1 tibble_1.4.2
[9] ggplot2_3.0.0 tidyverse_1.2.1
loaded via a namespace (and not attached):
[1] Rcpp_0.12.19 cellranger_1.1.0 plyr_1.8.4
[4] compiler_3.5.1 pillar_1.3.0 git2r_0.23.0
[7] bindr_0.1.1 R.methodsS3_1.7.1 R.utils_2.7.0
[10] tools_3.5.1 digest_0.6.17 lubridate_1.7.4
[13] jsonlite_1.5 evaluate_0.11 nlme_3.1-137
[16] gtable_0.2.0 lattice_0.20-35 pkgconfig_2.0.2
[19] rlang_0.2.2 cli_1.0.1 rstudioapi_0.8
[22] yaml_2.2.0 haven_1.1.2 bindrcpp_0.2.2
[25] withr_2.1.2 xml2_1.2.0 httr_1.3.1
[28] knitr_1.20 hms_0.4.2 rprojroot_1.3-2
[31] grid_3.5.1 tidyselect_0.2.4 glue_1.3.0
[34] R6_2.3.0 readxl_1.1.0 rmarkdown_1.10
[37] modelr_0.1.2 magrittr_1.5 whisker_0.3-2
[40] backports_1.1.2 scales_1.0.0 htmltools_0.3.6
[43] rvest_0.3.2 assertthat_0.2.0 colorspace_1.3-2
[46] stringi_1.2.4 lazyeval_0.2.1 munsell_0.5.0
[49] broom_0.5.0 crayon_1.3.4 R.oo_1.22.0
This reproducible R Markdown analysis was created with workflowr 1.1.1