Pipeline for peak coverage

Create bedgraph and bigwig:

I can do this step in my snakefile. First, I added the following to my environemnt.

• ucsc-bedgraphtobigwig
  
• ucsc-bigwigmerge
  
• ucsc-wigtobigwig
  
• ucsc-bigwigtobedgraph

I want to create bedgraph for each file. I will add a rule to my snakefile that does this and puts them in the bedgraph directory.

#add to directory
dir_bedgraph= dir_data + "bedgraph/"

#add to rule_all  

expand(dir_bedgraph + "{samples}.bg", samples=samples)

#rule
rule bedgraph: 
  input:
    bam = dir_sort + "{samples}-sort.bam"
  output: dir_bedgraph + "{samples}.bg"
  shell: "bedtools genomecov -ibam {input.bam} -bg -5 > {output}"

I want to add more memory for this rule in the cluster.json

"bedgraph" :
    {
            "mem": 16000
    }

I will use the bedgraphtobigwig tool.

#add to directory
dir_bigwig= dir_data + "bigwig/"
dir_sortbg= dir_data + "bedgraph_sort/"

#add to rule_all  
expand(dir_sortbg + "{samples}.sort.bg", samples=samples)
expand(dir_bigwig + "{samples}.bw", samples=samples)

rule sort_bg:
    input: dir_bedgraph + "{samples}.bg"
    output: dir_sortbg + "{samples}.sort.bg"
    shell: "sort -k1,1 -k2,2n {input} > {output}"

rule bg_to_bw:
    input: 
        bg=dir_sortbg + "{samples}.sort.bg"
        len= chrom_length 
    output: dir_bigwig + "{samples}.bw"
    shell: "bedGraphToBigWig {input.bg} {input.len} {output}""

Merge BW

This next step will take all of the files in the bigwig directory and merge them. To do this I will create a script that creates a list of all of the files then uses this list in the merge script.

mergeBW.sh

#!/bin/bash

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

module load Anaconda3
source activate three-prime-env

ls -d -1 /project2/gilad/briana/threeprimeseq/data/bigwig/* | tail -n +2 > /project2/gilad/briana/threeprimeseq/data/list_bw/list_of_bigwig.txt

bigWigMerge -inList /project2/gilad/briana/threeprimeseq/data/list_bw/list_of_bigwig.txt /project2/gilad/briana/threeprimeseq/data/mergedBW/merged_combined_YL-SP-threeprimeseq.bg

The result of this script will be a merged bedgraph of all of the files.

Convert to coverage

library(workflowr)

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

library(ggplot2)
library(dplyr)

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

update the following: callPeaksYL_combdata.py

Merge until the coverage is less than a number (ex. 5). Then start again when it is greater than that number. Previously we used the cuttoff 5. I can look at the distribution of the coverage to pick an informative cuttoff.

merged_bg=read.table("../data/merged_combined_YL-SP-threeprimeseq.bg", col.names=c("chr", "start", "end", "cov"))
summary(merged_bg$cov)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      1       1       3      52       9 5304660 

cov_plot=plot(sort(log10(merged_bg$cov)))
abline(h=1)

From this it looks like 10 is a good cuttoff for the peaks. I want to make a script that will convert the bedgraph to a coverage file.

#!/usr/bin/env python


main(inFile, outFile):
    fout = open(outFile,'w')
    for ind,ln in enumerate(open(inFile)):
      print(ind)
      chrom, start, end, count = ln.split()
      i2=int(start)
      while i2 < int(end):
        fout.write("%s\t%d\t%s\n"%(chrom, i2 + 1, count))
        fout.flush()
        i2 += 1
    fout.close()    
    

if __name__ == "__main__":
    import numpy as np
    from misc_helper import *
    import sys
    inFile = sys.argv[1]
    outFile = sys.argv[2]
    main(inFile, outFile)

Create a bash script to run this. I want the input and output files to be arguments in the python script.

#!/bin/bash

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

module load Anaconda3
source activate three-prime-env 

python bg_to_cov.py "/project2/gilad/briana/threeprimeseq/data/mergedBW/merged_combined_YL-SP-threeprimeseq.bg" "/project2/gilad/briana/threeprimeseq/data/mergedBW/merged_combined_YL-SP-threeprimeseq.coverage.txt"

Get peaks

Next I need to run Yangs peak script. This script is callPeaksYL_GEN.py

infile: “/project2/gilad/briana/threeprimeseq/data/mergedBW/merged_combined_YL-SP-threeprimeseq.coverage.txt”

Outfile: “/project2/gilad/briana/threeprimeseq/data/mergedPeaks/merged_combined_YL-SP-threeprimeseq.coverage.peaks.chr%s.bed”%chrom

The bash script runs this on each chromosome.

#!/bin/bash

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

module load Anaconda3
source activate three-prime-env


for i in $(seq 1 22); do 
  python callPeaksYL_GEN.py $i
done

Problem: all peak files are empty because my genome coverage file does not have 0s. In need to merge a genome coverage file with only 0s with the file I just made.

awk '{print $1 "\t" $2 "\t" "0"}' /project2/gilad/briana/threeprimeseq/data/bedgraph_comb/NuclearBamFiles.split.genomecov.bed > /project2/gilad/briana/threeprimeseq/data/mergedBW/genomecov_zero.txt 

The next script will sort the coverage file and merge it:

add_zero.R

#!/bin/rscripts

# usage: ./addZeros.R in_cov, zero_cov, outfile

#this script takes in the coverage file without 0s and adds the rows with 0s

library(optparse)
library(dplyr)
library(tidyr)
library(ggplot2)



option_list = list(
  make_option(c("-f", "--file"), action="store", default=NA, type='character',
              help="input coverage file"),
  make_option(c("-z", "--zero"), action="store", default=NA, type='character',
              help="gencov file 0"),
  make_option(c("-o", "--output"), action="store", default=NA, type='character',
              help="output file")
)


opt_parser <- OptionParser(option_list=option_list)
opt <- parse_args(opt_parser)


#interrupt execution if no file is  supplied
if (is.null(opt$file)){
  print_help(opt_parser)
  stop("Need input file", call.=FALSE)
}

names=c("chr", "base", "cov")
cov_file=read.table(file = opt$file, sep="\t",  col.names=names)
zero_file=read.table(file=opt$zero, sep="\t",  col.names=names)

anti=anti_join(zero_file, cov_file, by=c("chr", "base"))  

final=full_join(cov_file, anti, by=c("chr", "base", "cov"))

write.table(final, file=opt$output, quote = F, row.names = F, col.names = F)  

Test this method:

names=c("chr", "base", "cov")
cov=read.table("../data/gencov.test.csv", col.names=names, header=F, sep=",")

Warning in read.table("../data/gencov.test.csv", col.names = names,
header = F, : incomplete final line found by readTableHeader on '../data/
gencov.test.csv'

zero=read.table("../data/gencov_zero.test.csv",col.names=names, header=F, sep=",")


anti=anti_join(zero, cov, by=c("chr", "base"))  

final=full_join(cov, anti, by=c("chr", "base", "cov")) 

run_add_zero.sh

#!/bin/bash

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

module load Anaconda3
source activate three-prime-env

sort -k1,1 -k2,2n /project2/gilad/briana/threeprimeseq/data/mergedBW/merged_combined_YL-SP-threeprimeseq.coverage.txt > /project2/gilad/briana/threeprimeseq/data/mergedBW/merged_combined_YL-SP-threeprimeseq.coverage.sort.txt

Rscript add_zero.R -f /project2/gilad/briana/threeprimeseq/data/mergedBW/merged_combined_YL-SP-threeprimeseq.coverage.sort.txt -z /project2/gilad/briana/threeprimeseq/data/mergedBW/genomecov_zero.txt -o /project2/gilad/briana/threeprimeseq/data/mergedBW/merged_combined_YL-SP-threeprimeseq.coverage.sort.with0.txt

This next script will sort the final gencov file sort_gencov.sh :

#!/bin/bash

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

module load Anaconda3
source activate three-prime-env

sort -k1,1 -k2,2n  /project2/gilad/briana/threeprimeseq/data/mergedBW/merged_combined_YL-SP-threeprimeseq.coverage.sort.with0.txt > /project2/gilad/briana/threeprimeseq/data/mergedBW/merged_combined_YL-SP-threeprimeseq.coverage.sort.with0.sort.txt 

Try this with bash:

#!/bin/bash

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

module load Anaconda3
source activate three-prime-env

sort -k1,1 -k2,2n /project2/gilad/briana/threeprimeseq/data/mergedBW/merged_combined_YL-SP-threeprimeseq.coverage.txt > /project2/gilad/briana/threeprimeseq/data/mergedBW/merged_combined_YL-SP-threeprimeseq.coverage.sort.txt

less /project2/gilad/briana/threeprimeseq/data/mergedBW/merged_combined_YL-SP-threeprimeseq.coverage.sort.txt | awk '{print($1"."$2"\t"$3)}' >  /project2/gilad/briana/threeprimeseq/data/mergedBW/temp1.txt

less /project2/gilad/briana/threeprimeseq/data/mergedBW/genomecov_zero.txt | awk '{print($1"."$2"\t"$3)}' >  /project2/gilad/briana/threeprimeseq/data/mergedBW/temp2.txt

join  -a1  -a2 -o '0,1.2' -e 0   /project2/gilad/briana/threeprimeseq/data/mergedBW/temp1.txt /project2/gilad/briana/threeprimeseq/data/mergedBW/temp.2.txt | tr '.' '\t' | tr ' ' '\t' | cut -f1-4 > /project2/gilad/briana/threeprimeseq/data/mergedBW/merged_combined_YL-SP-threeprimeseq.coverage.sort.with0.bash.txt

Run yangs scrip on /project2/gilad/briana/threeprimeseq/data/mergedBW/merged_combined_YL-SP-threeprimeseq.coverage.sort.with0.sort.txt by making this the input file in the callPeaksYL_GEN.py

Session information

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] dplyr_0.7.6     ggplot2_3.0.0   workflowr_1.1.1

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.18      compiler_3.5.1    pillar_1.3.0     
 [4] git2r_0.23.0      plyr_1.8.4        bindr_0.1.1      
 [7] R.methodsS3_1.7.1 R.utils_2.6.0     tools_3.5.1      
[10] digest_0.6.15     evaluate_0.11     tibble_1.4.2     
[13] gtable_0.2.0      pkgconfig_2.0.1   rlang_0.2.1      
[16] rstudioapi_0.7    yaml_2.1.19       bindrcpp_0.2.2   
[19] withr_2.1.2       stringr_1.3.1     knitr_1.20       
[22] rprojroot_1.3-2   grid_3.5.1        tidyselect_0.2.4 
[25] glue_1.3.0        R6_2.2.2          rmarkdown_1.10   
[28] purrr_0.2.5       magrittr_1.5      whisker_0.3-2    
[31] backports_1.1.2   scales_0.5.0      htmltools_0.3.6  
[34] assertthat_0.2.0  colorspace_1.3-2  stringi_1.2.4    
[37] lazyeval_0.2.1    munsell_0.5.0     crayon_1.3.4     
[40] R.oo_1.22.0