Chromatin marks and PAS

Last updated: 2019-03-25

Checks: 6 0

Knit directory: threeprimeseq/analysis/

This reproducible R Markdown analysis was created with workflowr (version 1.2.0). The Report tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.

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.

Cache: none

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

Repository version: f3740ea

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility. The version displayed above was the version of the Git repository at the time these results were generated.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use 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:    data/perm_QTL_trans_noMP_5percov/
    Ignored:    output/.DS_Store

Untracked files:
    Untracked:  KalistoAbundance18486.txt
    Untracked:  analysis/4suDataIGV.Rmd
    Untracked:  analysis/AdaptTonyExampleplots.Rmd
    Untracked:  analysis/DirectionapaQTL.Rmd
    Untracked:  analysis/EvaleQTLs.Rmd
    Untracked:  analysis/YL_QTL_test.Rmd
    Untracked:  analysis/groSeqAnalysis.Rmd
    Untracked:  analysis/ncbiRefSeq_sm.sort.mRNA.bed
    Untracked:  analysis/peaksWithMisprimming.Rmd
    Untracked:  analysis/snake.config.notes.Rmd
    Untracked:  analysis/verifyBAM.Rmd
    Untracked:  analysis/verifybam_dubs.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/AllPeak_counts/
    Untracked:  data/ApaQTLs/
    Untracked:  data/ApaQTLs_otherPhen/
    Untracked:  data/CTCF/
    Untracked:  data/ChromHmmOverlap/
    Untracked:  data/DistTXN2Peak_genelocAnno/
    Untracked:  data/EmpiricalDists/
    Untracked:  data/FeatureoverlapPeaks/
    Untracked:  data/GM12878.chromHMM.bed
    Untracked:  data/GM12878.chromHMM.txt
    Untracked:  data/GWAS_overlap/
    Untracked:  data/LianoglouLCL/
    Untracked:  data/LocusZoom/
    Untracked:  data/LocusZoom_Unexp/
    Untracked:  data/LocusZoom_proc/
    Untracked:  data/MatchedSnps/
    Untracked:  data/NucSpecQTL/
    Untracked:  data/NuclearApaQTLs.txt
    Untracked:  data/PeakCounts/
    Untracked:  data/PeakCounts_noMP_5perc/
    Untracked:  data/PeakCounts_noMP_genelocanno/
    Untracked:  data/PeakUsage/
    Untracked:  data/PeakUsage_noMP/
    Untracked:  data/PeakUsage_noMP_GeneLocAnno/
    Untracked:  data/PeaksUsed/
    Untracked:  data/PeaksUsed_noMP_5percCov/
    Untracked:  data/PolyA_DB/
    Untracked:  data/QTL_overlap/
    Untracked:  data/RNAdecay/
    Untracked:  data/RNAkalisto/
    Untracked:  data/RefSeq_annotations/
    Untracked:  data/Replicates_usage/
    Untracked:  data/Signal_Loc/
    Untracked:  data/TotalApaQTLs.txt
    Untracked:  data/Totalpeaks_filtered_clean.bed
    Untracked:  data/UnderstandPeaksQC/
    Untracked:  data/WASP_STAT/
    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/apaExamp_proc/
    Untracked:  data/apaQTL_examp_noMP/
    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_GeneLocAnno/
    Untracked:  data/diff_iso_proc/
    Untracked:  data/diff_iso_trans/
    Untracked:  data/eQTL_inAPA/
    Untracked:  data/eQTLs_Lietal/
    Untracked:  data/ensemble_to_genename.txt
    Untracked:  data/example_gene_peakQuant/
    Untracked:  data/explainProtVar/
    Untracked:  data/filtPeakOppstrand_cov_noMP_GeneLocAnno_5perc/
    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/molPheno_noMP/
    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_inAPA/
    Untracked:  data/pQTL_otherphen/
    Untracked:  data/pacbio_cov/
    Untracked:  data/peakPerRefSeqGene/
    Untracked:  data/peaks4DT/
    Untracked:  data/perm_QTL/
    Untracked:  data/perm_QTL_GeneLocAnno_noMP_5percov/
    Untracked:  data/perm_QTL_GeneLocAnno_noMP_5percov_3UTR/
    Untracked:  data/perm_QTL_diffWindow/
    Untracked:  data/perm_QTL_opp/
    Untracked:  data/perm_QTL_trans/
    Untracked:  data/perm_QTL_trans_filt/
    Untracked:  data/protAndAPAAndExplmRes.Rda
    Untracked:  data/protAndAPAlmRes.Rda
    Untracked:  data/protAndExpressionlmRes.Rda
    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:  data/threePrimeSeqMetaData55Ind.txt
    Untracked:  data/threePrimeSeqMetaData55Ind.xlsx
    Untracked:  data/threePrimeSeqMetaData55Ind_noDup.txt
    Untracked:  data/threePrimeSeqMetaData55Ind_noDup.xlsx
    Untracked:  data/threePrimeSeqMetaData55Ind_noDup_WASPMAP.txt
    Untracked:  data/threePrimeSeqMetaData55Ind_noDup_WASPMAP.xlsx
    Untracked:  data/threePrimeSeqMetaData55Ind_redobatch4.xlsx
    Untracked:  data/~$threePrimeSeqMetaData55Ind_redobatch4.xlsx
    Untracked:  manuscript/
    Untracked:  output/LZ/
    Untracked:  output/deeptools_plots/
    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/NewPeakPostMP.Rmd
    Modified:   analysis/NuclearSpecQTL.Rmd
    Modified:   analysis/PeakToXper.Rmd
    Modified:   analysis/RNAdecayAndAPA.Rmd
    Modified:   analysis/apaQTLoverlapGWAS.Rmd
    Modified:   analysis/characterize_apaQTLs.Rmd
    Modified:   analysis/cleanupdtseq.internalpriming.Rmd
    Modified:   analysis/coloc_apaQTLs_protQTLs.Rmd
    Modified:   analysis/dif.iso.usage.leafcutter.Rmd
    Modified:   analysis/diffIsoAnalysisNewMapping.Rmd
    Modified:   analysis/diff_iso_pipeline.Rmd
    Modified:   analysis/explainpQTLs.Rmd
    Modified:   analysis/explore.filters.Rmd
    Modified:   analysis/fixBWChromNames.Rmd
    Modified:   analysis/flash2mash.Rmd
    Modified:   analysis/initialPacBioQuant.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/pipeline_55Ind.Rmd
    Modified:   analysis/swarmPlots_QTLs.Rmd
    Modified:   analysis/test.max2.Rmd
    Modified:   analysis/test.smash.Rmd
    Modified:   analysis/understandPeaks.Rmd
    Modified:   analysis/unexplainedeQTL_analysis.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.

These are the previous versions of the R Markdown and HTML files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view them.

File	Version	Author	Date	Message
Rmd	f3740ea	Briana Mittleman	2019-03-25	add top and bottom dnase
html	3eac464	Briana Mittleman	2019-03-23	Build site.
Rmd	4939dfe	Briana Mittleman	2019-03-23	look at mnase at all categories
html	1b7f088	Briana Mittleman	2019-03-21	Build site.
Rmd	c02e927	Briana Mittleman	2019-03-21	add mnase merge chipseq
html	a6b0fe4	Briana Mittleman	2019-03-20	Build site.
Rmd	54168fd	Briana Mittleman	2019-03-20	add histone mod analysis

library(tidyverse)

── Attaching packages ──────────────────────────────────────────────────────── tidyverse 1.2.1 ──

✔ ggplot2 3.1.0       ✔ purrr   0.3.1  
✔ tibble  2.0.1       ✔ dplyr   0.8.0.1
✔ tidyr   0.8.3       ✔ stringr 1.4.0  
✔ readr   1.3.1       ✔ forcats 0.4.0

Warning: package 'tibble' was built under R version 3.5.2

Warning: package 'tidyr' was built under R version 3.5.2

Warning: package 'purrr' was built under R version 3.5.2

Warning: package 'dplyr' was built under R version 3.5.2

Warning: package 'stringr' was built under R version 3.5.2

Warning: package 'forcats' was built under R version 3.5.2

── Conflicts ─────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()

library(workflowr)

This is workflowr version 1.2.0
Run ?workflowr for help getting started

library(data.table)

Warning: package 'data.table' was built under R version 3.5.2


Attaching package: 'data.table'

The following objects are masked from 'package:dplyr':

    between, first, last

The following object is masked from 'package:purrr':

    transpose

library(cowplot)

Warning: package 'cowplot' was built under R version 3.5.2


Attaching package: 'cowplot'

The following object is masked from 'package:ggplot2':

    ggsave

http://science.sciencemag.org/content/352/6291/aad9926.full?ijkey=fkp/DIzVNS9RY&keytype=ref&siteid=sci

This article talks about chromatin modifications for heterochromatin downstream of PAS. I will look at enrichment for repressive histone marks downstream of my called PAS.

Repressive marks H3K27me3, H3K9me3

http://hgdownload.soe.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeUwHistone/

H3k27me3 H3k36me3 H3k4me3

Deeptools plot

h3k27me3DTmypeaks.sh

#!/bin/bash

#SBATCH --job-name=h3k27me3DTmypeaks
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=h3k27me3DTmypeaks.out
#SBATCH --error=h3k27me3DTmypeaks.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/ChipSeq/wgEncodeUwHistoneGm06990H3k27me3StdRawRep1.bigWig -R /project2/gilad/briana/threeprimeseq/data/peaks4DT/APAPeaks_5percCov_fixedStrand.bed -b 1000 -a 1000  -out /project2/gilad/briana/threeprimeseq/data/ChipSeq/H3K27me3_myPeaksNompfilt.gz


plotHeatmap --sortRegions descend -m /project2/gilad/briana/threeprimeseq/data/ChipSeq/H3K27me3_myPeaksNompfilt.gz --refPointLabel "Called PAS" --plotTitle "H3K27me3 at All Called PAS" --heatmapHeight 7 --colorMap YlGnBu  -out /project2/gilad/briana/threeprimeseq/data/ChipSeq/H3K27me3_myPeaksNompfilt.png

Put all of the marks on one plot:

I also want to just use the last base of the peak APAPAS_5percCov_fixedStrand.bed histonemarksDTmypeaks.sh

#!/bin/bash

#SBATCH --job-name=histonemarksDTmypeaks
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=histonemarksDTmypeaks.out
#SBATCH --error=histonemarksDTmypeaks.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/ChipSeq/wgEncodeUwHistoneGm06990H3k27me3StdRawRep1.bigWig /project2/gilad/briana/threeprimeseq/data/ChipSeq/wgEncodeUwHistoneGm06990H3k36me3StdRawRep1.bigWig  /project2/gilad/briana/threeprimeseq/data/ChipSeq/wgEncodeUwHistoneGm06990H3k4me3StdRawRep1.bigWig -R /project2/gilad/briana/threeprimeseq/data/peaks4DT/APAPAS_5percCov_fixedStrand.bed -b 1000 -a 1000  -out /project2/gilad/briana/threeprimeseq/data/ChipSeq/HistMarks_myPeaksNompfilt.gz


plotHeatmap --sortRegions descend -m /project2/gilad/briana/threeprimeseq/data/ChipSeq/HistMarks_myPeaksNompfilt.gz --refPointLabel "Called PAS" --plotTitle "Histone marks at All Called PAS" --heatmapHeight 7 --colorMap YlGnBu  -out /project2/gilad/briana/threeprimeseq/data/ChipSeq/HistMarks_myPeaksNompfilt.png

Scales are too different to put these on the same spot:

H3k27me3DTmyPAS.sh

#!/bin/bash

#SBATCH --job-name=H3k27me3DTmypeaks.sh
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=H3k27me3DTmypeaks.out
#SBATCH --error=H3k27me3DTmypeaks.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/ChipSeq/MergedGm06990H3k27me3.bw  -R /project2/gilad/briana/threeprimeseq/data/peaks4DT/APAPAS_5percCov_fixedStrand.bed -b 5000 -a 5000  -out /project2/gilad/briana/threeprimeseq/data/ChipSeq/H3k27me3_myPeaksNompfilt.gz --outFileNameMatrix /project2/gilad/briana/threeprimeseq/data/ChipSeq/H3k27me3_myPeaksNompfilt_matrix.txt --outFileSortedRegions /project2/gilad/briana/threeprimeseq/data/ChipSeq/H3k27me3_myPeaksNompfilt_SortedRegions.txt 


plotHeatmap --sortRegions descend -m /project2/gilad/briana/threeprimeseq/data/ChipSeq/H3k27me3_myPeaksNompfilt.gz --refPointLabel "Called PAS" --plotTitle "H3k27me3 at All Called PAS" --heatmapHeight 7 --colorMap YlGnBu  -out /project2/gilad/briana/threeprimeseq/data/ChipSeq/H3k27me3_myPeaksNompfilt.png

H3k36me3DTmyPAS.sh

#!/bin/bash

#SBATCH --job-name=H3k27me3DTmypeaks.sh
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=H3k36me3DTmypeaks.out
#SBATCH --error=H3k36me3DTmypeaks.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/ChipSeq/MergedGm06990H3k36me3.bw   -R /project2/gilad/briana/threeprimeseq/data/peaks4DT/APAPAS_5percCov_fixedStrand.bed -b 5000 -a 5000  -out /project2/gilad/briana/threeprimeseq/data/ChipSeq/H3k36me3_myPeaksNompfilt.gz --outFileNameMatrix /project2/gilad/briana/threeprimeseq/data/ChipSeq/H3k36me3_myPeaksNompfilt_matrix.txt --outFileSortedRegions /project2/gilad/briana/threeprimeseq/data/ChipSeq/H3k36me3_myPeaksNompfilt_SortedRegions.txt 


plotHeatmap --sortRegions descend -m /project2/gilad/briana/threeprimeseq/data/ChipSeq/H3k36me3_myPeaksNompfilt.gz --refPointLabel "Called PAS" --plotTitle "H3k36me3 at All Called PAS" --heatmapHeight 7 --colorMap YlGnBu  -out /project2/gilad/briana/threeprimeseq/data/ChipSeq/H3k36me3_myPeaksNompfilt.png

H3k4me3DTmyPAS.sh

#!/bin/bash

#SBATCH --job-name=H3k4me3DTmypeaks.sh
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=H3k36me3DTmypeaks.out
#SBATCH --error=H3k36me3DTmypeaks.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/ChipSeq/MergedGm06990H3k4me3.bw  -R /project2/gilad/briana/threeprimeseq/data/peaks4DT/APAPAS_5percCov_fixedStrand.bed -b 5000 -a 5000  -out /project2/gilad/briana/threeprimeseq/data/ChipSeq/H3k4me3_myPeaksNompfilt.gz --outFileNameMatrix project2/gilad/briana/threeprimeseq/data/ChipSeq/H3k4me3_myPeaksNompfilt_matrix.txt --outFileSortedRegions /project2/gilad/briana/threeprimeseq/data/ChipSeq/H3k4me3_myPeaksNompfilt_SortedRegions.txt 


plotHeatmap --sortRegions descend -m /project2/gilad/briana/threeprimeseq/data/ChipSeq/H3k4me3_myPeaksNompfilt.gz --refPointLabel "Called PAS" --plotTitle "H3k4me3 at All Called PAS" --heatmapHeight 7 --colorMap YlGnBu  -out /project2/gilad/briana/threeprimeseq/data/ChipSeq/H3k4me3_myPeaksNompfilt.png

Download both replicates for these and merge:

mergeH3k27me3.sh


#!/bin/bash

#SBATCH --job-name=mergeH3k27me3
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=mergeH3k27me3.out
#SBATCH --error=mergeH3k27me3.err
#SBATCH --partition=broadwl
#SBATCH --mem=36G
#SBATCH --mail-type=END
module load Anaconda3
source activate three-prime-env


bigWigMerge /project2/gilad/briana/threeprimeseq/data/ChipSeq/wgEncodeUwHistoneGm06990H3k27me3StdRawRep1.bigWig /project2/gilad/briana/threeprimeseq/data/ChipSeq/wgEncodeUwHistoneGm06990H3k27me3StdRawRep2.bigWig /project2/gilad/briana/threeprimeseq/data/ChipSeq/MergedGm06990H3k27me3.bedGraph


sort -k1,1 -k2,2n /project2/gilad/briana/threeprimeseq/data/ChipSeq/MergedGm06990H3k27me3.bedGraph /project2/gilad/briana/threeprimeseq/data/ChipSeq/MergedGm06990H3k27me3.sort.bedGraph

bedGraphToBigWig /project2/gilad/briana/threeprimeseq/data/ChipSeq/MergedGm06990H3k27me3.sort.bedGraph /project2/gilad/briana/genome_anotation_data/chrom.length.chr.txt  /project2/gilad/briana/threeprimeseq/data/ChipSeq/MergedGm06990H3k27me3.bw

mergeH3k36me3.sh


#!/bin/bash

#SBATCH --job-name=mergeH3k36me3
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=mergeH3k36me3.out
#SBATCH --error=mergeH3k36me3.err
#SBATCH --partition=broadwl
#SBATCH --mem=36G
#SBATCH --mail-type=END
module load Anaconda3
source activate three-prime-env


bigWigMerge /project2/gilad/briana/threeprimeseq/data/ChipSeq/wgEncodeUwHistoneGm06990H3k36me3StdRawRep1.bigWig /project2/gilad/briana/threeprimeseq/data/ChipSeq/wgEncodeUwHistoneGm06990H3k36me3StdRawRep2.bigWig /project2/gilad/briana/threeprimeseq/data/ChipSeq/MergedGm06990H3k36me3.bedGraph


sort -k1,1 -k2,2n /project2/gilad/briana/threeprimeseq/data/ChipSeq/MergedGm06990H3k36me3.bedGraph > /project2/gilad/briana/threeprimeseq/data/ChipSeq/MergedGm06990H3k36me3.sort.bedGraph

bedGraphToBigWig /project2/gilad/briana/threeprimeseq/data/ChipSeq/MergedGm06990H3k36me3.sort.bedGraph /project2/gilad/briana/genome_anotation_data/chrom.length.chr.txt  /project2/gilad/briana/threeprimeseq/data/ChipSeq/MergedGm06990H3k36me3.bw

mergeH3k4me3.sh


#!/bin/bash

#SBATCH --job-name=mergeH3k4me3
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=mergeH3k4me3.out
#SBATCH --error=mergeH3k4me3.err
#SBATCH --partition=broadwl
#SBATCH --mem=36G
#SBATCH --mail-type=END
module load Anaconda3
source activate three-prime-env


bigWigMerge /project2/gilad/briana/threeprimeseq/data/ChipSeq/wgEncodeUwHistoneGm06990H3k4me3StdRawRep1.bigWig /project2/gilad/briana/threeprimeseq/data/ChipSeq/wgEncodeUwHistoneGm06990H3k4me3StdRawRep2.bigWig /project2/gilad/briana/threeprimeseq/data/ChipSeq/MergedGm06990H3k4me3.bedGraph


sort -k1,1 -k2,2n /project2/gilad/briana/threeprimeseq/data/ChipSeq/MergedGm06990H3k4me3.bedGraph >/project2/gilad/briana/threeprimeseq/data/ChipSeq/MergedGm06990H3k4me3.sort.bedGraph

bedGraphToBigWig /project2/gilad/briana/threeprimeseq/data/ChipSeq/MergedGm06990H3k4me3.sort.bedGraph /project2/gilad/briana/genome_anotation_data/chrom.length.chr.txt  /project2/gilad/briana/threeprimeseq/data/ChipSeq/MergedGm06990H3k4me3.bw

MNASE:

MNASEmyPAS.sh

#!/bin/bash

#SBATCH --job-name=MNASEmyPAS
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=MNASEmyPAS.out
#SBATCH --error=MNASEmyPAS.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/ChipSeq/ENCFF000VME.bigWig  -R /project2/gilad/briana/threeprimeseq/data/peaks4DT/APAPAS_5percCov_fixedStrand.bed -b 1000 -a 1000  -out /project2/gilad/briana/threeprimeseq/data/ChipSeq/mnase_myPeaksNompfilt.gz 


plotHeatmap --sortRegions descend -m /project2/gilad/briana/threeprimeseq/data/ChipSeq/mnase_myPeaksNompfilt.gz --refPointLabel "Called PAS" --plotTitle "MNASE at All Called PAS" --heatmapHeight 7 --colorMap YlGnBu  -out /project2/gilad/briana/threeprimeseq/data/ChipSeq/mnase_myPeaksNompfilt.png

Try second mnase track.

MNASEmyPAS_secondfile.sh


#!/bin/bash

#SBATCH --job-name=MNASEmyPAS_secondfile
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=MNASEmyPAS_secondfile.out
#SBATCH --error=MNASEmyPAS_secondfile.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/ChipSeq/wgEncodeSydhNsomeGm12878Sig.bigWig  -R /project2/gilad/briana/threeprimeseq/data/peaks4DT/APAPAS_5percCov_fixedStrand.bed -b 1000 -a 1000  -out /project2/gilad/briana/threeprimeseq/data/ChipSeq/mnase2_myPeaksNompfilt.gz 


plotHeatmap --sortRegions descend -m /project2/gilad/briana/threeprimeseq/data/ChipSeq/mnase2_myPeaksNompfilt.gz --refPointLabel "Called PAS" --plotTitle "MNASE at All Called PAS" --heatmapHeight 7 --colorMap YlGnBu  -out /project2/gilad/briana/threeprimeseq/data/ChipSeq/mnase2_myPeaksNompfilt.png

Convert to PAS ratehr than peak: APAPeaks_5percCov_fixedStrand_INTRON.bed

python Peak2PAS.py /project2/gilad/briana/threeprimeseq/data/peaks4DT/APAPeaks_5percCov_fixedStrand_INTRON.bed /project2/gilad/briana/threeprimeseq/data/peaks4DT/APAPAS_5percCov_fixedStrand_INTRON.bed

Run this with intronic vs utr
MNASEmyPASIntron.sh

#!/bin/bash

#SBATCH --job-name=MNASEmyPASIntron
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=MNASEmyPASIntron.out
#SBATCH --error=MNASEmyPASIntron.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/ChipSeq/ENCFF000VME.bigWig  -R  /project2/gilad/briana/threeprimeseq/data/peaks4DT/APAPAS_5percCov_fixedStrand_INTRON.bed -b 1000 -a 1000  -out /project2/gilad/briana/threeprimeseq/data/ChipSeq/mnase_myPeaksIntronNompfilt.gz 


plotHeatmap --sortRegions descend -m /project2/gilad/briana/threeprimeseq/data/ChipSeq/mnase_myPeaksIntronNompfilt.gz --refPointLabel "Called PAS" --plotTitle "MNASE at Intronic Called PAS" --heatmapHeight 7 --colorMap YlGnBu  -out /project2/gilad/briana/threeprimeseq/data/ChipSeq/mnase_myPeaksIntronNompfilt.png

Nuclear specific:

python Peak2PAS.py /project2/gilad/briana/threeprimeseq/data/peaks4DT/APAPeaks_5percCov_fixedStrand_SigUsageNuc.bed /project2/gilad/briana/threeprimeseq/data/peaks4DT/APAPAS_5percCov_fixedStrand_SigUsageNuc.bed

MNASEmyPASNuclear.sh

#!/bin/bash

#SBATCH --job-name=MNASEmyPASNuclear
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=MNASEmyPASNuclear.out
#SBATCH --error=MNASEmyPASNuclear.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/ChipSeq/ENCFF000VME.bigWig  -R  /project2/gilad/briana/threeprimeseq/data/peaks4DT/APAPAS_5percCov_fixedStrand_SigUsageNuc.bed -b 1000 -a 1000  -out /project2/gilad/briana/threeprimeseq/data/ChipSeq/mnase_myPeaksNuclearNompfilt.gz 


plotHeatmap --sortRegions descend -m /project2/gilad/briana/threeprimeseq/data/ChipSeq/mnase_myPeaksNuclearNompfilt.gz --refPointLabel "Called PAS" --plotTitle "MNASE at Nuclear specific PAS" --heatmapHeight 7 --colorMap YlGnBu  -out /project2/gilad/briana/threeprimeseq/data/ChipSeq/mnase_myPeaksNuclearNompfilt.png

Nuclear Intronic:

APAPeaks_5percCov_fixedStrand_SigUsageNuc_Intron.bed

python Peak2PAS.py /project2/gilad/briana/threeprimeseq/data/peaks4DT/APAPeaks_5percCov_fixedStrand_SigUsageNuc_Intron.bed /project2/gilad/briana/threeprimeseq/data/peaks4DT/APAPAS_5percCov_fixedStrand_SigUsageNuc_Intron.bed

MNASEmyPASNuclearIntronic.sh

#!/bin/bash

#SBATCH --job-name=MNASEmyPASNuclearIntronic
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=MNASEmyPASNuclearIntronic.out
#SBATCH --error=MNASEmyPASNuclearIntronic.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/ChipSeq/ENCFF000VME.bigWig  -R  /project2/gilad/briana/threeprimeseq/data/peaks4DT/APAPAS_5percCov_fixedStrand_SigUsageNuc_Intron.bed -b 1000 -a 1000  -out /project2/gilad/briana/threeprimeseq/data/ChipSeq/mnase_myPeaksNuclearIntronNompfilt.gz 


plotHeatmap --sortRegions descend -m /project2/gilad/briana/threeprimeseq/data/ChipSeq/mnase_myPeaksNuclearIntronNompfilt.gz --refPointLabel "Called PAS" --plotTitle "MNASE at Nuclear specific PAS in Intron" --heatmapHeight 7 --colorMap YlGnBu  -out /project2/gilad/briana/threeprimeseq/data/ChipSeq/mnase_myPeaksNuclearIntronNompfilt.png

Sepatate by mean usage (strong/weak)

Total

Mean usage:

top 20% and bottom 20%. by mean usage

meanUsageTot=read.table("../data/PeaksUsed_noMP_5percCov/filtered_APApeaks_merged_allchrom_refseqGenes.GeneLocAnno.NoMP_sm_quant.Total_fixed.pheno.5percPeaks.txt", col.names = c("chr", "start", "end","gene", "strand", "name", "meanUsage")) %>% mutate(perc=ntile(meanUsage,n=100)) 


meanUsageTot_bot20=meanUsageTot %>% filter(perc <20) %>% dplyr::select(name)
meanUsageTot_top20=meanUsageTot %>% filter(perc >80)%>% dplyr::select(name)

Write out the peaks:

write.table(meanUsageTot_bot20, file="../data/PeaksUsed_noMP_5percCov/TotalPeaksBottom20Usage.txt", quote=F, row.names = F, col.names = F, sep="\t")

write.table(meanUsageTot_top20, file="../data/PeaksUsed_noMP_5percCov/TotalPeaksTop20Usage.txt", quote=F, row.names = F, col.names = F, sep="\t")

Copy to /project2/gilad/briana/threeprimeseq/data/PeakUsage_noMP_GeneLocAnno/

Subset PAS file:

subsetPAStottop20perc.py

top20=open("/project2/gilad/briana/threeprimeseq/data/PeakUsage_noMP_GeneLocAnno/TotalPeaksTop20Usage.txt", "r")


AllPas=open("/project2/gilad/briana/threeprimeseq/data/peaks4DT/APAPAS_5percCov_fixedStrand.bed", "r")

Top20PAS=open("/project2/gilad/briana/threeprimeseq/data/peaks4DT/APAPAS_5percCov_fixedStrand_top20inTotal.bed", "w")





def subsetPAS(use, outpas, PAS):
    okPAS={}
    for ln in use:
        peak=ln.strip()
        okPAS[peak]=""
    for ln in PAS:
        peaknum=ln.split()[3].split(":")[-1]
        print
        if peaknum in okPAS.keys():
            outpas.write(ln)
    outpas.close()


subsetPAS(top20, Top20PAS, AllPas)

subsetPAStotbottom20perc.py


bottom20=open("/project2/gilad/briana/threeprimeseq/data/PeakUsage_noMP_GeneLocAnno/TotalPeaksBottom20Usage.txt", "r")

AllPas=open("/project2/gilad/briana/threeprimeseq/data/peaks4DT/APAPAS_5percCov_fixedStrand.bed", "r")


Bottom20PAS=open("/project2/gilad/briana/threeprimeseq/data/peaks4DT/APAPAS_5percCov_fixedStrand_bottom20inTotal.bed", "w")




def subsetPAS(use, outpas, PAS):
    okPAS={}
    for ln in use:
        peak=ln.strip()
        okPAS[peak]=""
    for ln in PAS:
        peaknum=ln.split()[3].split(":")[-1]
        print
        if peaknum in okPAS.keys():
            outpas.write(ln)
    outpas.close()


subsetPAS(bottom20, Bottom20PAS, AllPas)

Deeptools plots for these:

MNASEmyPAStop20tot.sh

#!/bin/bash

#SBATCH --job-name=MNASEmyPAStop20tot
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=MNASEmyPAStop20tot.out
#SBATCH --error=MNASEmyPAStop20tot.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/ChipSeq/ENCFF000VME.bigWig  -R  /project2/gilad/briana/threeprimeseq/data/peaks4DT/APAPAS_5percCov_fixedStrand_top20inTotal.bed -b 1000 -a 1000  -out /project2/gilad/briana/threeprimeseq/data/ChipSeq/mnase_myPeakstop20tot.gz 


plotHeatmap --sortRegions descend -m /project2/gilad/briana/threeprimeseq/data/ChipSeq/mnase_myPeakstop20tot.gz --refPointLabel "Called PAS" --plotTitle "MNASE at Top 20% Total Usage" --heatmapHeight 7 --colorMap YlGnBu  -out /project2/gilad/briana/threeprimeseq/data/ChipSeq/mnase_myPeakstop20tot.png

MNASEmyPASbottom20tot.sh

#!/bin/bash

#SBATCH --job-name=MNASEmyPASbottom20tot
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=MNASEmyPASbottom20tot.out
#SBATCH --error=MNASEmyPASbottom20tot.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/ChipSeq/ENCFF000VME.bigWig  -R  /project2/gilad/briana/threeprimeseq/data/peaks4DT/APAPAS_5percCov_fixedStrand_bottom20inTotal.bed -b 1000 -a 1000  -out /project2/gilad/briana/threeprimeseq/data/ChipSeq/mnase_myPeaksbottom20tot.gz 


plotHeatmap --sortRegions descend -m /project2/gilad/briana/threeprimeseq/data/ChipSeq/mnase_myPeaksbottom20tot.gz --refPointLabel "Called PAS" --plotTitle "MNASE at Bottom 20% Total Usage" --heatmapHeight 7 --colorMap YlGnBu  -out /project2/gilad/briana/threeprimeseq/data/ChipSeq/mnase_myPeaksbottom20tot.png

Nuclear

Mean usage:

top 20% and bottom 20%. by mean usage

meanUsageNuc=read.table("../data/PeaksUsed_noMP_5percCov/filtered_APApeaks_merged_allchrom_refseqGenes.GeneLocAnno.NoMP_sm_quant.Nuclear_fixed.pheno.5percPeaks.txt", col.names = c("chr", "start", "end","gene", "strand", "name", "meanUsage")) %>% mutate(perc=ntile(meanUsage,n=100)) 


meanUsageNuc_bot20=meanUsageNuc %>% filter(perc <20) %>% dplyr::select(name)
meanUsageNuc_top20=meanUsageNuc %>% filter(perc >80)%>% dplyr::select(name)

Write out the peaks:

write.table(meanUsageNuc_bot20, file="../data/PeaksUsed_noMP_5percCov/NuclearPeaksBottom20Usage.txt", quote=F, row.names = F, col.names = F, sep="\t")

write.table(meanUsageNuc_top20, file="../data/PeaksUsed_noMP_5percCov/NuclearPeaksTop20Usage.txt", quote=F, row.names = F, col.names = F, sep="\t")

Copy to /project2/gilad/briana/threeprimeseq/data/PeakUsage_noMP_GeneLocAnno/

Subset PAS file:

subsetPASnuctop20perc.py

top20=open("/project2/gilad/briana/threeprimeseq/data/PeakUsage_noMP_GeneLocAnno/NuclearPeaksTop20Usage.txt", "r")


AllPas=open("/project2/gilad/briana/threeprimeseq/data/peaks4DT/APAPAS_5percCov_fixedStrand.bed", "r")

Top20PAS=open("/project2/gilad/briana/threeprimeseq/data/peaks4DT/APAPAS_5percCov_fixedStrand_top20inNuclear.bed", "w")





def subsetPAS(use, outpas, PAS):
    okPAS={}
    for ln in use:
        peak=ln.strip()
        okPAS[peak]=""
    for ln in PAS:
        peaknum=ln.split()[3].split(":")[-1]
        print
        if peaknum in okPAS.keys():
            outpas.write(ln)
    outpas.close()


subsetPAS(top20, Top20PAS, AllPas)

subsetPASnucbottom20perc.py


bottom20=open("/project2/gilad/briana/threeprimeseq/data/PeakUsage_noMP_GeneLocAnno/NuclearPeaksBottom20Usage.txt", "r")

AllPas=open("/project2/gilad/briana/threeprimeseq/data/peaks4DT/APAPAS_5percCov_fixedStrand.bed", "r")


Bottom20PAS=open("/project2/gilad/briana/threeprimeseq/data/peaks4DT/APAPAS_5percCov_fixedStrand_bottom20inNuclear.bed", "w")




def subsetPAS(use, outpas, PAS):
    okPAS={}
    for ln in use:
        peak=ln.strip()
        okPAS[peak]=""
    for ln in PAS:
        peaknum=ln.split()[3].split(":")[-1]
        print
        if peaknum in okPAS.keys():
            outpas.write(ln)
    outpas.close()


subsetPAS(bottom20, Bottom20PAS, AllPas)

Deeptools plots for these:

MNASEmyPAStop20nuc.sh

#!/bin/bash

#SBATCH --job-name=MNASEmyPAStop20nuc
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=MNASEmyPAStop20nuc.out
#SBATCH --error=MNASEmyPAStop20nuc.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/ChipSeq/ENCFF000VME.bigWig  -R  /project2/gilad/briana/threeprimeseq/data/peaks4DT/APAPAS_5percCov_fixedStrand_top20inNuclear.bed -b 1000 -a 1000  -out /project2/gilad/briana/threeprimeseq/data/ChipSeq/mnase_myPeakstop20nuc.gz 


plotHeatmap --sortRegions descend -m /project2/gilad/briana/threeprimeseq/data/ChipSeq/mnase_myPeakstop20nuc.gz --refPointLabel "Called PAS" --plotTitle "MNASE at Top 20% Nuclear Usage" --heatmapHeight 7 --colorMap YlGnBu  -out /project2/gilad/briana/threeprimeseq/data/ChipSeq/mnase_myPeakstop20nuc.png

MNASEmyPASbottom20nuc.sh

#!/bin/bash

#SBATCH --job-name=MNASEmyPASbottom20nuc
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=MNASEmyPASbottom20nuc.out
#SBATCH --error=MNASEmyPASbottom20nuc.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/ChipSeq/ENCFF000VME.bigWig  -R  /project2/gilad/briana/threeprimeseq/data/peaks4DT/APAPAS_5percCov_fixedStrand_bottom20inNuclear.bed -b 1000 -a 1000  -out /project2/gilad/briana/threeprimeseq/data/ChipSeq/mnase_myPeaksbottom20nuc.gz 


plotHeatmap --sortRegions descend -m /project2/gilad/briana/threeprimeseq/data/ChipSeq/mnase_myPeaksbottom20nuc.gz --refPointLabel "Called PAS" --plotTitle "MNASE at Bottom 20% Nuclear Usage" --heatmapHeight 7 --colorMap YlGnBu  -out /project2/gilad/briana/threeprimeseq/data/ChipSeq/mnase_myPeaksbottom20nuc.png

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] cowplot_0.9.4     data.table_1.12.0 workflowr_1.2.0  
 [4] forcats_0.4.0     stringr_1.4.0     dplyr_0.8.0.1    
 [7] purrr_0.3.1       readr_1.3.1       tidyr_0.8.3      
[10] tibble_2.0.1      ggplot2_3.1.0     tidyverse_1.2.1  

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