Last updated: 2019-01-18
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: ee0c781
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/DirectionapaQTL.Rmd
Untracked: analysis/EvaleQTLs.Rmd
Untracked: analysis/PreAshExplore.Rmd
Untracked: analysis/YL_QTL_test.Rmd
Untracked: analysis/ncbiRefSeq_sm.sort.mRNA.bed
Untracked: analysis/snake.config.notes.Rmd
Untracked: analysis/verifyBAM.Rmd
Untracked: code/PeaksToCoverPerReads.py
Untracked: code/strober_pc_pve_heatmap_func.R
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/LianoglouLCL/
Untracked: data/LocusZoom/
Untracked: data/NuclearApaQTLs.txt
Untracked: data/PeakCounts/
Untracked: data/PeakUsage/
Untracked: data/PeaksUsed/
Untracked: data/RNAkalisto/
Untracked: data/TotalApaQTLs.txt
Untracked: data/Totalpeaks_filtered_clean.bed
Untracked: data/UnderstandPeaksQC/
Untracked: data/YL-SP-18486-T-combined-genecov.txt
Untracked: data/YL-SP-18486-T_S9_R1_001-genecov.txt
Untracked: data/YL_QTL_test/
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/example_gene_peakQuant/
Untracked: data/explainProtVar/
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/nuc_10up/
Untracked: data/other_qtls/
Untracked: data/pQTL_otherphen/
Untracked: data/peakPerRefSeqGene/
Untracked: data/perm_QTL/
Untracked: data/perm_QTL_opp/
Untracked: data/perm_QTL_trans/
Untracked: data/perm_QTL_trans_filt/
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: data/threePrimeSeqMetaData.csv
Untracked: output/picard/
Untracked: output/plots/
Untracked: output/qual.fig2.pdf
Unstaged changes:
Modified: analysis/28ind.peak.explore.Rmd
Modified: analysis/CompareLianoglouData.Rmd
Modified: analysis/apaQTLoverlapGWAS.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/explainpQTLs.Rmd
Modified: analysis/explore.filters.Rmd
Modified: analysis/flash2mash.Rmd
Modified: analysis/mispriming_approach.Rmd
Modified: analysis/overlapMolQTL.Rmd
Modified: analysis/overlapMolQTL.opposite.Rmd
Modified: analysis/overlap_qtls.Rmd
Modified: analysis/peakOverlap_oppstrand.Rmd
Modified: analysis/peakQCPPlots.Rmd
Modified: analysis/pheno.leaf.comb.Rmd
Modified: analysis/swarmPlots_QTLs.Rmd
Modified: analysis/test.max2.Rmd
Modified: analysis/understandPeaks.Rmd
Modified: code/Snakefile
I want to show RNA seq vs 3’ seq in peaks that are internal.
To do this I need to pull in the peaks (use the filtered ones from peakQCPlots) to get peaks that are used but are not the most distal for the gene. I need to split by strand when i do this becasue most distal is different in both cases.
Load Libraries
library(data.table)
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::between() masks data.table::between()
✖ dplyr::filter() masks stats::filter()
✖ dplyr::first() masks data.table::first()
✖ dplyr::lag() masks stats::lag()
✖ dplyr::last() masks data.table::last()
✖ purrr::transpose() masks data.table::transpose()library(workflowr)This is workflowr version 1.1.1
Run ?workflowr for help getting startedlibrary(cowplot)
Attaching package: 'cowplot'The following object is masked from 'package:ggplot2':
ggsaveLoad peaks for total and nuclear post filter:
total_PeakUsMean=read.table("../data/PeakUsage/PeakUsageMeanPost5percFilter.Total.txt", head = T)
nuclear_PeakUsMean=read.table("../data/PeakUsage/PeakUsageMeanPost5percFilter.Nuclear.txt", head=T)Seperate positive and negative:
total_PeakUsMean_pos=total_PeakUsMean %>% filter(strand=="+")
total_PeakUsMean_neg=total_PeakUsMean %>% filter(strand=="-")nuclear_PeakUsMean_pos=nuclear_PeakUsMean %>% filter(strand=="+")
nuclear_PeakUsMean_neg=nuclear_PeakUsMean %>% filter(strand=="-")Group by gene and keep internal (remove genes with only 1)
total_PeakUsMean_pos_internal=total_PeakUsMean_pos %>% group_by(gene) %>% mutate(n=n()) %>% filter(n>1) %>% top_n(1,peak)
nuclear_PeakUsMean_pos_internal=nuclear_PeakUsMean_pos %>% group_by(gene) %>% mutate(n=n()) %>% filter(n>1) %>% top_n(1,peak)*For negative strand use top_n(-1)
total_PeakUsMean_neg_internal=total_PeakUsMean_neg %>% group_by(gene) %>% mutate(n=n()) %>% filter(n>1) %>% top_n(-1,peak)
nuclear_PeakUsMean_neg_internal=nuclear_PeakUsMean_neg %>% group_by(gene) %>% mutate(n=n()) %>% filter(n>1) %>% top_n(-1,peak)Bind the total and nuclear rows back together:
Total
total_PeakUsMean_internal=as.data.frame(rbind(total_PeakUsMean_pos_internal,total_PeakUsMean_neg_internal)) %>% arrange(chr, start, end)
write.table(total_PeakUsMean_internal, file="../data/PeakUsage/InternalFilteredPeak.total.txt", row.names = F, col.names = F, quote = F)Nuclear:
nuclear_PeakUsMean_internal=as.data.frame(rbind(nuclear_PeakUsMean_pos_internal,nuclear_PeakUsMean_neg_internal)) %>% arrange(chr, start, end)
write.table(nuclear_PeakUsMean_internal, file="../data/PeakUsage/InternalFilteredPeak.nuclear.txt", row.names = F, col.names = F, quote = F)Use these peak numbers to filter the bed file that I use for deep tools. I can do this in python with a dictionary of the peaks to keep.
The file I want to subset is /project2/gilad/briana/threeprimeseq/data/mergedPeaks_comb/filtered_APApeaks_merged_allchrom_refseqTrans.noties_sm.fixed.bed
filterPeaksInner.py
#python
peaksIn="/project2/gilad/briana/threeprimeseq/data/mergedPeaks_comb/filtered_APApeaks_merged_allchrom_refseqTrans.noties_sm.fixed.bed"
innerPeaksTotal="/project2/gilad/briana/threeprimeseq/data/PeakUsage/InternalFilteredPeak.total.txt"
innerPeaksNuclear="/project2/gilad/briana/threeprimeseq/data/PeakUsage/InternalFilteredPeak.nuclear.txt"
OkPeakTot={}
for ln in open(innerPeaksTotal, "r"):
peaknum=ln.split()[5]
peaknumOnly=peaknum[4:]
OkPeakTot[peaknumOnly] = ""
OkPeakNuc={}
for ln in open(innerPeaksNuclear, "r")
peaknum=ln.split()[5]
peaknumOnly=peaknum[4:]
OkPeakNuc[peaknumOnly]=""
TotOut=open("/project2/gilad/briana/threeprimeseq/data/mergedPeaks_comb/filtered_APApeaks_merged_allchrom_refseqTrans.noties_sm.fixed_TotalInternal.bed", "w")
NucOut=open("/project2/gilad/briana/threeprimeseq/data/mergedPeaks_comb/filtered_APApeaks_merged_allchrom_refseqTrans.noties_sm.fixed_NuclearInternal.bed", "w")
#read peaks and write out to each file
for ln in open(peaksIn,"r"):
peakName=ln.split()[3]
peakNameStr=str(peakName)
if peakName in OkPeakTot.keys():
TotOut.write(ln)
if peakName in OkPeakNuc.keys():
NucOut.write(ln)
TotOut.close()
NucOut.close()
Now I can make the deep tools plots:
Three prime seq with inner total
TotalDTPlotmyPeaks_Internal.sh
#!/bin/bash
#SBATCH --job-name=TotalDTPlotmyPeaks_Internal
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=TotalDTPlotmyPeaks_Internal.out
#SBATCH --error=TotalDTPlotmyPeaks_Internal.err
#SBATCH --partition=bigmem2
#SBATCH --mem=100G
#SBATCH --mail-type=END
module load Anaconda3
source activate three-prime-env
computeMatrix reference-point -S /project2/gilad/briana/threeprimeseq/data/mergedBW/Total_MergedBamCoverage.bw -R /project2/gilad/briana/threeprimeseq/data/mergedPeaks_comb/filtered_APApeaks_merged_allchrom_refseqTrans.noties_sm.fixed_TotalInternal.bed -b 1000 -a 1000 --outFileName /project2/gilad/briana/threeprimeseq/data/LianoglouDeepTools/Total_myPeaks_Internal.gz
plotHeatmap --sortRegions descend -m /project2/gilad/briana/threeprimeseq/data/LianoglouDeepTools/Total_myPeaks_Internal.gz --refPointLabel "Internal Peaks" --plotTitle "Total Combined Reads Internal Peaks" --heatmapHeight 7 --colorMap YlGnBu -out /project2/gilad/briana/threeprimeseq/data/LianoglouDeepTools/Total_myPeaks_internal.png
NuclearDTPlotmyPeaks_Internal.sh
#!/bin/bash
#SBATCH --job-name=NuclearDTPlotmyPeaks_Internal
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=NucelarDTPlotmyPeaks_Internal.out
#SBATCH --error=NucelarDTPlotmyPeaks_Internal.err
#SBATCH --partition=bigmem2
#SBATCH --mem=100G
#SBATCH --mail-type=END
module load Anaconda3
source activate three-prime-env
computeMatrix reference-point -S /project2/gilad/briana/threeprimeseq/data/mergedBW/Nucelar_MergedBamCoverage.bw -R /project2/gilad/briana/threeprimeseq/data/mergedPeaks_comb/filtered_APApeaks_merged_allchrom_refseqTrans.noties_sm.fixed_NucelarInternal.bed -b 1000 -a 1000 --outFileName /project2/gilad/briana/threeprimeseq/data/LianoglouDeepTools/Nuclear_myPeaks_Internal.gz
plotHeatmap --sortRegions descend -m /project2/gilad/briana/threeprimeseq/data/LianoglouDeepTools/Nuclear_myPeaks_Internal.gz --refPointLabel "Internal Peaks" --plotTitle "Nuclear Combined Reads Internal Peaks" --heatmapHeight 7 --colorMap YlGnBu -out /project2/gilad/briana/threeprimeseq/data/LianoglouDeepTools/Nuclear_myPeaks_internal.png
RNA seq plot with the total internal peaks
RNAseqDTPlotmyPeaks_Internal.sh
#!/bin/bash
#SBATCH --job-name=RNAseqDTPlotmyPeaks_Internal
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=RNAseqDTPlotmyPeaks_Internal.out
#SBATCH --error=RNAseqDTPlotmyPeaks_Internal.err
#SBATCH --partition=bigmem2
#SBATCH --mem=100G
#SBATCH --mail-type=END
module load Anaconda3
source activate three-prime-env
computeMatrix reference-point -S /project2/gilad/briana/threeprimeseq/data/rnaseq_bw/RNAseqGeuvadis_STAR_6samp_MergedBams.sort.bw -R /project2/gilad/briana/threeprimeseq/data/mergedPeaks_comb/filtered_APApeaks_merged_allchrom_refseqTrans.noties_sm.fixed_TotalInternal.bed -b 1000 -a 1000 -out /project2/gilad/briana/threeprimeseq/data/rnaseq_deeptools/RNAseq_myPeaks_Internal.gz
plotHeatmap --sortRegions descend -m /project2/gilad/briana/threeprimeseq/data/rnaseq_deeptools/RNAseq_myPeaks_Internal.gz --refPointLabel "Total Internal Peaks" --plotTitle "Combined RNAseq Reads at Total Internal" --heatmapHeight 7 --colorMap YlGnBu -out /project2/gilad/briana/threeprimeseq/data/rnaseq_deeptools/RNAseq_myPeaks_Internal.pngRNAseqDTPlotmyPeaks_NuclearInternal.sh
#!/bin/bash
#SBATCH --job-name=RNAseqDTPlotmyPeaks_NuclearInternal
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=RNAseqDTPlotmyPeaks_NuclearInternal.out
#SBATCH --error=RNAseqDTPlotmyPeaks_NuclearInternal.err
#SBATCH --partition=bigmem2
#SBATCH --mem=100G
#SBATCH --mail-type=END
module load Anaconda3
source activate three-prime-env
computeMatrix reference-point -S /project2/gilad/briana/threeprimeseq/data/rnaseq_bw/RNAseqGeuvadis_STAR_6samp_MergedBams.sort.bw -R /project2/gilad/briana/threeprimeseq/data/mergedPeaks_comb/filtered_APApeaks_merged_allchrom_refseqTrans.noties_sm.fixed_NuclearInternal.bed -b 1000 -a 1000 -out /project2/gilad/briana/threeprimeseq/data/rnaseq_deeptools/RNAseq_myPeaks_NuclearInternal.gz
plotHeatmap --sortRegions descend -m /project2/gilad/briana/threeprimeseq/data/rnaseq_deeptools/RNAseq_myPeaks_NuclearInternal.gz --refPointLabel "Nucelar Internal Peaks" --plotTitle "Combined RNAseq Reads at Nuclear Internal" --heatmapHeight 7 --colorMap YlGnBu -out /project2/gilad/briana/threeprimeseq/data/rnaseq_deeptools/RNAseq_myPeaks_NuclearInternal.pngsessionInfo()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 cowplot_0.9.3 workflowr_1.1.1
[4] forcats_0.3.0 stringr_1.3.1 dplyr_0.7.6
[7] purrr_0.2.5 readr_1.1.1 tidyr_0.8.1
[10] tibble_1.4.2 ggplot2_3.0.0 tidyverse_1.2.1
[13] data.table_1.11.8
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.19 backports_1.1.2 scales_1.0.0
[28] jsonlite_1.5 hms_0.4.2 digest_0.6.17
[31] stringi_1.2.4 grid_3.5.1 rprojroot_1.3-2
[34] cli_1.0.1 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.8 R6_2.3.0 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