deNovo peak callling

Create Bedgraph

I will call peaks de novo in the combined total and nuclear fraction 3’ Seq. The data is reletevely clean so I will start with regions that have continuous coverage. I will first create a bedgraph.

#!/bin/bash

#SBATCH --job-name=Tbedgraph
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=Tbedgraph.out
#SBATCH --error=Tbedgraph.err
#SBATCH --partition=broadwl
#SBATCH --mem=40G
#SBATCH --mail-type=END


module load Anaconda3
source activate three-prime-env 

samtools sort -o /project2/gilad/briana/threeprimeseq/data/macs2/TotalBamFiles.sort.bam /project2/gilad/briana/threeprimeseq/data/macs2/TotalBamFiles.bam

bedtools genomecov -ibam /project2/gilad/briana/threeprimeseq/data/macs2/TotalBamFiles.sort.bam -bga > /project2/gilad/briana/threeprimeseq/data/bedgraph/TotalBamFiles.bedgraph

Next I will create the file without the 0 places in the genome. I will be able to use this for the bedtools merge function.

awk '{if ($4 != 0) print}' TotalBamFiles.bedgraph >TotalBamFiles_no0.bedgraph 

I can merge the regions with consequtive reads using the bedtools merge function.

• -i input bed

• -c colomn to act on

• -o collapse, print deliminated list of the counts from -c call

• -delim “,”

This is the mergeBedgraph.sh script. It takes in the no 0 begraph filename without the path.

#!/bin/bash

#SBATCH --job-name=merge
#SBATCH --account=pi-yangili1
#SBATCH --time=8:00:00
#SBATCH --output=merge.out
#SBATCH --error=merge.err
#SBATCH --partition=broadwl
#SBATCH --mem=16G
#SBATCH --mail-type=END

module load Anaconda3
source activate three-prime-env 

bedgraph=$1
describer=$(echo ${bedgraph} | sed -e "s/.bedgraph$//")

bedtools merge -c 4,4,4 -o count,mean,collapse -delim "," -i /project2/gilad/briana/threeprimeseq/data/bedgraph/$1 > /project2/gilad/briana/threeprimeseq/data/bedgraph/${describer}.peaks.bed

Run this first on the total bedgraph, TotalBamFiles_no0.bedgraph. The file has chromosome, start, end, number of regions, mean, and a string of the values.

This is not exaclty what I want. I need to go back and do genome cov not collapsing with bedgraph.

To evaluate this I will bring the file into R and plot some statistics about it.

#!/bin/bash

#SBATCH --job-name=Tgencov
#SBATCH --account=pi-yangili1
#SBATCH --time=24:00:00
#SBATCH --output=Tgencov.out
#SBATCH --error=Tgencov.err
#SBATCH --partition=broadwl
#SBATCH --mem=40G
#SBATCH --mail-type=END


module load Anaconda3
source activate three-prime-env 


bedtools genomecov -ibam /project2/gilad/briana/threeprimeseq/data/macs2/TotalBamFiles.sort.bam -d > /project2/gilad/briana/threeprimeseq/data/bedgraph/TotalBamFiles.genomecov.bed

Evaluate regions

First I will look at the bedgraph file. This is not as imformative becuase it combined regions with the same counts.

total_bedgraph=read.table("../data/bedgraph_peaks/TotalBamFiles_no0.peaks.bed",col.names = c("chr", "start", "end", "regions", "mean", "counts"))

Plot the mean:

plot(sort(log10(total_bedgraph$mean), decreasing=T), xlab="Region", ylab="log10 of bedgraph region bin", main="Distribution of log10 region means from bedgraph")

I want to look at the distribution of how many bases are included in the regions.

library(dplyr)

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

Attaching package: 'dplyr'

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

    filter, lag

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

    intersect, setdiff, setequal, union

Tregion_bases=total_bedgraph %>% mutate(bases=end-start) %>% select(bases)

Warning: package 'bindrcpp' was built under R version 3.4.4

plot(sort(log10(Tregion_bases$bases), decreasing = T), xlab="Region", ylab="log10 of region size", main="Distribution of bases in regions- log10")

Given the reads are abotu 60bp this is probably pretty good.

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     

other attached packages:
[1] bindrcpp_0.2.2 dplyr_0.7.5   

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.17      knitr_1.18        bindr_0.1.1      
 [4] whisker_0.3-2     magrittr_1.5      workflowr_1.0.1  
 [7] tidyselect_0.2.4  R6_2.2.2          rlang_0.2.1      
[10] stringr_1.3.1     tools_3.4.2       R.oo_1.22.0      
[13] git2r_0.21.0      htmltools_0.3.6   yaml_2.1.19      
[16] rprojroot_1.3-2   digest_0.6.15     assertthat_0.2.0 
[19] tibble_1.4.2      purrr_0.2.5       R.utils_2.6.0    
[22] glue_1.2.0        evaluate_0.10.1   rmarkdown_1.8.5  
[25] stringi_1.2.2     pillar_1.1.0      compiler_3.4.2   
[28] backports_1.1.2   R.methodsS3_1.7.1 pkgconfig_2.0.1