Last updated: 2018-02-02

Code version: 01f846f

Deeptools

Deeptools is a way to make plots similar to genomation but they have enrichment plots and heatmaps together. Also this is in bash so I do not need to get the bam files into R.

Install

I added the deeptools package to my net-seq conda environment.

Running:

Step 1: Create bigwig coverage files with bamcoverage

  • bamCoverage -b reads.bam -o coverage.bw

Step 2: computeMatrix

  • I will need my normalized bigwig reads and the bed interval file (in my case PAS clusters)

  • ex: computeMatrix scale-regions -S -R -b 1000 -a 1000 -out

  • –skipZeros (option- not included in first try)

Step 3: Plot heatmap

  • required –matrixFile, -m (from the compute matrix), -out (file name to save image.png)

  • –sortRegions descending

  • –plotTitle, -T

I will run this first on the hES 3’ seq becasue we expect heavy enrichment.

#!/bin/bash


#SBATCH --job-name=deeptools_hES
#SBATCH --time=8:00:00
#SBATCH --partition=gilad
#SBATCH --mem=16G
#SBATCH --mail-type=END


module load Anaconda3  

source activate net-seq  

bamCoverage -b /project2/gilad/briana/Lianoglou.data/hES.hg38.sorted.bam -o hES.cov.bw  

computeMatrix reference-point -S hES.cov.bw -R /project2/gilad/briana/apa_sites/clusters.hg38.bed -b 1000 -a 1000 -out hES.hg38.gz

plotHeatmap --sortRegions descend -m hES.hg38.gz  -out hES.hg38.apa.png

This plot shows enrichement of the 3’ seq on called PAS sites. The axis are wrong.

hES 3’ coverage at APA sites

hES 3’ coverage at APA sites

I will make this plot excluding exons. I will use bedtools intersect to get APA sites not in exons.

I want A interects B -v. This will give me regions in A not in B. A is the cluster file. /project2/gilad/briana/apa_sites/clusters.hg38.bed and B is /project2/gilad/briana/apa_sites/deeptools/exons_hg38.bed

This bash file is /project2/gilad/briana/apa_sites/deeptools/heatmatrix_hES_apa_excludeExons.sh

#!/bin/bash


#SBATCH --job-name=deeptools_hES_noexon
#SBATCH --time=8:00:00
#SBATCH --partition=gilad
#SBATCH --mem=16G
#SBATCH --mail-type=END


module load Anaconda3  

source activate net-seq  

bedtools intersect -v -a /project2/gilad/briana/apa_sites/clusters.hg38.bed -b /project2/gilad/yangili/hg38_exons.bed > ../clusters.hg38.noExons.bed 


#bamCoverage -b /project2/gilad/briana/Lianoglou.data/hES.hg38.sorted.bam -o hES.cov.bw
# dont need to redo this. The file exsists

computeMatrix reference-point -S hES.cov.bw -R /project2/gilad/briana/apa_sites/clusters.hg38.noExons.bed -b 1000 -a 1000 -out hES.hg38.noExons.gz

plotHeatmap --sortRegions descend -T "hES 3` APA enrichment exclude exons" -m hES.hg38.noExons.gz  -out hES.hg38.apa.noExons.png
hES 3’ coverage at APA sites not in exons

hES 3’ coverage at APA sites not in exons

Session information

sessionInfo()
R version 3.4.2 (2017-09-28)
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.4/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.4/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     

loaded via a namespace (and not attached):
 [1] compiler_3.4.2  backports_1.1.2 magrittr_1.5    rprojroot_1.3-2
 [5] tools_3.4.2     htmltools_0.3.6 yaml_2.1.16     Rcpp_0.12.15   
 [9] stringi_1.1.6   rmarkdown_1.8.5 knitr_1.18      git2r_0.21.0   
[13] stringr_1.2.0   digest_0.6.14   evaluate_0.10.1

This R Markdown site was created with workflowr