These are the exact implementation of the processing steps used in the manuscript
These FASTX and CTK wrappers were developed into CLIPflexR
You can see vignettes of CLIPflexR processing or CTK processing executed using CLIPflexR
Additionally, see CLIPflexR documentation for:
unzip
FASTX’s fastq_quality_filter
FASTX’s fastx_quality_stats
FASTX’s fastx_collapser
# #unzip fastqs
#
# gunzipKate <- function(fileToGunZip,gunzipTom="gunzip",
# stderr=file.path(getwd(),"gunzip_stderr"),
# stdout=file.path(getwd(),"gunzip_stdout")){
# cmd <- gunzipTom
# if(!file.exists(fileToGunZip))stop("File does not exist")
#
# fileWithoutExtension <- tools::file_path_sans_ext(fileToGunZip)
# if(file.exists(fileToGunZip) & !file.exists(fileWithoutExtension)){
# args <- c(fileToGunZip,"-k")
# message(cmd)
# message(args)
#
# system2(cmd,
# args,
# stdout=stdout,
# stderr=stderr
# )
# }
# return(fileWithoutExtension)
# }
#
# #R command to unzip:
# fileToGunZip <- "/Users/kathryn/Reprocess_all_paper_datasets/fastq_gz"
# fileToGunZip <- dir(fileToGunZip, pattern= "*.fastq.gz", full.names = TRUE)
#
# for (file in fileToGunZip) {
# gunzipKate(file)
# }
#
#
# #FastX quality filter:
#
# fastq_quality_filter <- function(fileTofqf,fqf="fastq_quality_filter",qEncoding=33,
# minimumQuality=20,
# minimumPercentOfRead=80,
# stderr=file.path(getwd(),"fastq_quality_filter_stderr"),
# stdout=file.path(getwd(),"fastq_quality_filter_stdout")){
#
#
# cmd <- fqf
#
# if(!file.exists(fileTofqf))stop("File does not exist")
#
# file_fqf <- file.path(dirname(fileTofqf),paste0("QF_",basename(fileTofqf)))
# if(file.exists(fileTofqf) & !file.exists(file_fqf)){
#
#
# args <- c(
# paste0("-Q ",qEncoding),
# paste0("-q ",minimumQuality),
# paste0("-p ",minimumPercentOfRead),
# paste0("-i ",fileTofqf),
# paste0("-o ",file_fqf)
# )
# print(cmd)
# print(args)
#
# system2(cmd,
# args,
# stdout=stdout,
# stderr=stderr
# )
# }
# return(file_fqf)
# }
#
#
# QF_dir <- "/Users/kathryn/Reprocess_all_paper_datasets"
# toQF <- dir(QF_dir, pattern= "fastq", full.names = TRUE)
#
# for (file in toQF) { fastq_quality_filter(file,
# fqf = "/Users/kathryn/FastX/fastq_quality_filter")
# }
#
#
# #get quality stats on fastq
#
# fastx_quality_stats <- function(fileTofqs,fqs="fastx_quality_stats",qEncoding=33,
# stderr=paste0(getwd(),"fastq_quality_stats_stderr"),
# stdout=paste0(getwd(),"fastq_quality_stats_stdout")){
# cmd <- fqs
#
# if(!file.exists(fileTofqs))stop("File does not exist")
#
# file_fqs <- file.path(dirname(fileTofqs),gsub("\\.fastq",".txt",basename(fileTofqs)))
# if(file.exists(fileTofqs) & !file.exists(file_fqs)){
#
#
# args <- c(
# paste0("-Q ",qEncoding),
# paste0("-i ",fileTofqs),
# paste0("-o ",file_fqs)
# )
# print(cmd)
# print(args)
#
# system2(cmd,
# args,
# stdout=stdout,
# stderr=stderr
# )
# }
# return(file_fqs)
# }
#
#
# stat_dir <- "/Users/kathryn/Reprocess_all_paper_datasets"
# to_stat <- dir(stat_dir, pattern= "fastq", full.names = TRUE)
#
# for (file in to_stat) { fastq_quality_filter(file,
# fqf = "/Users/kathryn/FastX/fastq_quality_filter")
# }
#
# #FastX collapse exact duplicates
#
# fastx_collapser <- function(fileTofxc,fxc="fastx_collapser",qEncoding=33,
# stderr=file.path(getwd(),"fastq_collapse_stderr"),
# stdout=file.path(getwd(),"fastq_collapse_stdout")){
# cmd <- fxc
#
# if(!file.exists(fileTofxc))stop("File does not exist")
#
# file_fxc <- file.path(dirname(fileTofxc),gsub("\\.fa","_collapse.fasta",basename(fileTofxc)))
# if(file.exists(fileTofxc) & !file.exists(file_fxc)){
#
# args <- c(
# paste0("-v "),
# paste0("-Q ",qEncoding),
# paste0("-i ",fileTofxc),
# paste0("-o ",file_fxc)
# )
# print(cmd)
# print(args)
#
# system2(cmd,
# args,
# stdout=stdout,
# stderr=stderr
# )
# }
# return(file_fxc)
# }
#
# pathtocollapse <- "/rugpfs/fs0/rice_lab/scratch/krozen/fastq_gz"
# tocollapse <- dir(pathtocollapse, pattern= "QF_*", full.names = TRUE)
#
# for (file in tocollapse) { fastq_quality_filter(file,
# fqf = "/Users/kathryn/FastX/fastx_collapser")
# }Additionally, see CLIPflexR documentation for:
FASTX’s fastx_barcode_splitter.pl
CTK’s stripBarcode.pl
FASTX’s fastx_clipper
# #FastX barcode_splitter to split samples by index:
#
# fastx_barcode_splitter <- function(fileTofxc,bcFile,mismatches=0,
# fxc="fastx_barcode_splitter.pl",
# stderr=file.path(getwd(),"fastx_barcode_splitter_stderr"),
# stdout=file.path(getwd(),"fastx_barcode_splitter_stdout")){
#
# cmd <- fxc
#
# if(!file.exists(fileTofxc))stop("File does not exist")
#
# prefix <- gsub("QF_|\\.fasta|\\.fastq","",basename(fileTofxc))
#
#
# cmd2 <- paste0("cat ",
# fileTofxc," ",
# "| ",
# cmd," ",
# " --bcfile ",bcFile," ",
# "--bol --mismatches ",mismatches," ",
# "--prefix '",prefix,"_' ")
# temp <- system(cmd2,wait = TRUE,intern = TRUE)
#
# return(temp)
# }
#
# #did this for each fasta with its paired BC.txt file
# pathtoBC <- "/Users/kathryn/Reprocess_all_paper_datasets/QF_in_collapse.fasta"
# bcFile <- "/Users/kathryn/Reprocess_all_paper_datasets/BC_files_final_only/BC.txt"
# fastx_barcode_splitter(pathtoBC, bcFile, fxc = "/Users/kathryn/FastX/fastx_barcode_splitter.pl")
#
# #CTK, strip 5' index, linker sequence/barcode:
#
# stripBarcode <- function(filesToRun,
# sb="stripBarcode.pl",
# PATHTOPERLLIB=NULL,
# stderr=file.path(getwd(),"stripBarcode_stderr"),
# stdout=file.path(getwd(),"stripBarcode_stdout"),
# linkerlength=NULL){
# fileToRun <- filesToRun[1]
# cmd <- sb
#
# if(!file.exists(fileToRun)) stop("File does not exist")
#
# exportPATH <- ifelse(!is.null(PATHTOPERLLIB),paste0("export PERL5LIB=",PATHTOPERLLIB,";"),"")
#
# cmd2 <- paste0(exportPATH," ",
# cmd," ",
# " -len ",linkerlength," -v ",
# " ",fileToRun," ",
# paste(fileToRun,"_rm5",sep=""))
# temp <- system(cmd2,wait = TRUE,intern = TRUE)
#
# return(temp)
# }
#
# Dir <- "/Users/kathryn/Reprocess_all_paper_datasets"
# strip <- dir(Dir, pattern= "KRG*", full.names = TRUE)
#
# for (file in strip) { stripBarcode(file,
# sb="/Users/kathryn/CTK/stripBarcode.pl", PATHTOPERLLIB = "/Users/kathryn/czplib", linkerlength = 27)
# }
#
# #FastX clip 3' linker and keep reads at least 18nt:
#
# fastx_clipper <- function(fileTofqs,fqc="fastx_clipper",length=18,
# adaptor="GTGTCAGTCACTTCCAGCGG",
# stderr=file.path(getwd(),"clipper_stats_stderr"),
# stdout=file.path(getwd(),"clipper_stats_stdout")){
# cmd <- fqc
#
# if(!file.exists(fileTofqs))stop("File does not exist")
#
# file_fqs <- file.path(dirname(fileTofqs),gsub("rm5","rm5_rm3.fa",basename(fileTofqs)))
# if(file.exists(fileTofqs) & !file.exists(file_fqs)){
#
# args <- c(
# paste0("-l ",length),
# paste0("-a ",adaptor),
# paste0("-o ",file_fqs),
# paste0("-i ",fileTofqs)
# )
# print(cmd)
# print(args)
#
# system2(cmd,
# args,
# stdout=stdout,
# stderr=stderr
# )
# }
# return(file_fqs)
# }
#
# Dir <- "/Users/kathryn/Reprocess_all_paper_datasets"
# clip <- dir(Dir, pattern= "*_rm5", full.names = TRUE)
# for (file in clip) { fastx_clipper(file,fqc = "/Users/kathryn/FastX/fastx_clipper")
# }Additionally, see CLIPflexR documentation for:
CTK’s stripBarcode.pl
FASTX’s fastx_barcode_splitter.pl
FASTX’s fastx_clipper
FASTX’s fastx_trimmer
# #CTK, strip 5' index, linker sequence/barcode:
# stripBarcode <- function(filesToRun,
# sb="stripBarcode.pl",
# PATHTOPERLLIB=NULL,
# stderr=file.path(getwd(),"stripBarcode_stderr"),
# stdout=file.path(getwd(),"stripBarcode_stdout"),
# linkerlength=NULL){
# fileToRun <- filesToRun[1]
# cmd <- sb
#
# if(!file.exists(fileToRun)) stop("File does not exist")
#
# exportPATH <- ifelse(!is.null(PATHTOPERLLIB),paste0("export PERL5LIB=",PATHTOPERLLIB,";"),"")
#
# cmd2 <- paste0(exportPATH," ",
# cmd," ",
# " -len ",linkerlength," -v ",
# " ",fileToRun," ",
# paste(fileToRun,"_rm5",sep=""))
# #"temp.rm")
# temp <- system(cmd2,wait = TRUE,intern = TRUE)
#
# return(temp)
# }
#
#
# Dir <- "/Users/kathryn/Reprocess_all_paper_datasets"
# strip <- dir(Dir, pattern= "BrdU*", full.names = TRUE)
#
# for (file in strip) { stripBarcode(file,
# sb="/Users/kathryn/CTK/stripBarcode.pl", PATHTOPERLLIB = "/Users/kathryn/czplib", linkerlength = 7)
# }
#
# #FastX barcode_splitter to split samples by index:
# #did this for each fasta with its paired BC.txt file
# pathtoBC <- "/Users/kathryn/Reprocess_all_paper_datasets/QF_in_collapse.fasta"
# bcFile <- "/Users/kathryn/Reprocess_all_paper_datasets/BC_files_final_only/BC.txt"
# fastx_barcode_splitter(pathtoBC, bcFile, fxc = "/Users/kathryn/FastX/fastx_barcode_splitter.pl")
#
# #original trim and clip terminal command
# #R command:
#
# fastx_trimmer <- function(fileTofqt,fqt="fastx_trimmer",length=10,
# stderr=file.path(getwd(),"trimmer_stats_stderr"),
# stdout=file.path(getwd(),"trimmer_stats_stdout")){
# cmd <- fqt
#
# if(!file.exists(fileTofqt))stop("File does not exist")
#
# file_fqt <- file.path(dirname(fileTofqt),paste0(basename(fileTofqt), "_rm5.fa"))
# if(file.exists(fileTofqt) & !file.exists(file_fqt)){
#
# args <- c(
# paste0("-f ",length),
# paste0("-o ",file_fqt),
# paste0("-i ",fileTofqt)
# )
# print(cmd)
# print(args)
#
# system2(cmd,
# args,
# stdout=stdout,
# stderr=stderr
# )
# }
# return(file_fqt)
# }
#
# Dir <- "/Users/kathryn/Reprocess_all_paper_datasets"
# trim <- dir(Dir, pattern= "BrdU*", full.names = TRUE)
# for (file in trim) { fastx_trimmer(file,fqt = "/Users/kathryn/FastX/fastx_trimmer")
# }
#
# Dir <- "/Users/kathryn/Reprocess_all_paper_datasets"
# clip <- dir(Dir, pattern= "*_rm5", full.names = TRUE)
# for (file in clip) { fastx_clipper(file,fqc = "/Users/kathryn/FastX/fastx_clipper")
# }Rbowtie2 indices were previously made: Ae. aegypti genome fasta (AaegL5) used to make index available at Vectorbase; they have since updated the genome and I’m unsure if they still have the older version, or what the differences are between these genomes; I can make this genome available upon request
NC_001564.2 Cell fusing agent virus (CFAV) strain Galveston, complete genome fasta is available from NCBI
Additionally, see CLIPflexR documentation for:
bowtie2_index
bowtie_align
bamtobed
# bowtie_align <- function(fq,index,sam=gsub("\\.fq|\\.fastq|\\.rm|\\.fa",".sam",fq)
# ) {
# require(Rbowtie2)
# if(!dir.exists(sam)){
# bowtie2(bt2Index = index,
# samOutput = sam,
# seq1 = fq,"--threads 4 -f -N 1 -L 18")
# Rsamtools::asBam(sam,gsub("\\.sam","",sam))
# }
#
# }
#
#
# Dir <- "/Users/kathryn/Reprocess_all_paper_datasets/BCsplit"
# fqfiles <- dir(Dir, pattern= "*_rm5_rm3.fa", full.names = TRUE)
# bplapply(fqfiles, bowtie_align, index = "/Users/kathryn/Bowtie_indices/CFAV")
# bplapply(fqfiles, bowtie_align, index = "/Users/kathryn/Bowtie_indices/AaegL5/Aedes-aegypti-LVP_AGWG_CHROMOSOMES_AaegL5")
#
# #convert to BED:
#
# bamtobed <- function(file,filtDup=FALSE){
# require(GenomicAlignments)
# require(rtracklayer)
# temp <- readGAlignments(file,
# param=ScanBamParam(what = "qname"))
# names(temp) <- mcols(temp)$qname
# temp <- granges(temp)
#
# if(filtDup) temp <- temp[!duplicated(temp),]
# export.bed(temp,
# con = gsub("\\.bam",".bed",file))
# }
#
# Dir <- "/Users/kathryn/Reprocess_all_paper_datasets/CFAV_mapped"
# mapped <- dir(Dir,pattern="bam$",full.names = TRUE)
# bplapply(mapped, bamtobed)
#
# Dir <- "/Users/kathryn/Reprocess_all_paper_datasets/AaegL5_mapped"
# mapped <- dir(Dir,pattern="bam$",full.names = TRUE)
# bplapply(mapped, bamtobed)Additionally, see HOMER and CLIPflexR documentation
# homer_peaks <- function(fileTofqs,maketagdir="makeTagDirectory",
# findpeaks="findpeaks",
# stderr=file.path(getwd(),"homer_stats_stderr"),
# stdout=file.path(getwd(),"homer_stats_stdout")){
# cmd <- maketagdir
#
# if(!file.exists(fileTofqs))stop("File does not exist")
# baseNAME <- make.names(basename(fileTofqs))
# tagDir <- file.path(dirname(fileTofqs),
# gsub("\\.bed","",baseNAME))
#
# if(file.exists(fileTofqs) & !dir.exists(tagDir)){
#
# dir.create(tagDir,showWarnings = TRUE,recursive = TRUE)
#
# args <- c(
# paste0(tagDir),
# paste0(fileTofqs),
# paste0("-single "),
# paste0("-format bed ")
# )
# print(cmd)
# print(args)
#
# system2(cmd,
# args,
# stdout=stdout,
# stderr=stderr
# )
# }
#
# cmd <- findpeaks
#
#
# if(dir.exists(tagDir)){
#
#
# args <- c(
# paste0(tagDir),
# paste0("-o auto -style factor -L 2 -localSize 10000 -strand separate -minDist 50 -size 10 -fragLength 25 -gsize 1278731969")
# )
# print(cmd)
# print(args)
#
# system2(cmd,
# args,
# stdout=stdout,
# stderr=stderr
# )
#
# }
# return("")
# }
#
# peakDir <- "/Users/kathryn/Reprocess_all_paper_datasets/AaegL5_mapped"
# bedfile <- dir(peakDir,pattern="*.bed$",full.names = TRUE)
# bedfile <- grep("unmatched", bedfile, invert= TRUE, value = TRUE)
#
# bplapply(bedfile, homer_peaks, maketagdir = "/Users/kathryn/homer/bin/makeTagDirectory", findpeaks = "/Users/kathryn/homer/bin/findPeaks" )
#
# #ran this command in terminal to rename peak txt files according to directory name
# ##for subdir in *; do mv $subdir/peaks.txt $subdir.txt; done;
#
# ##Note that this is updated in CLIPflexR (https://kathrynrozengagnon.github.io/CLIPflexR/), so there is no longer a need to reformat peaks.txt files# #make matrix with all datasets
# #read in peak files and make true/false matrix of peaks
# peakDir <- "/Users/kathryn/Reprocess_all_paper_datasets/AaegL5_mapped"
# peakToCall <- dir(peakDir,pattern="*.txt",full.names = TRUE)
# peakToCall <- grep("unmatched", peakToCall, invert= TRUE, value = TRUE)
# peakToCall <- grep("33243", peakToCall, invert= TRUE, value = TRUE)
# peakToCall <- grep("KRG040518", peakToCall, invert= TRUE, value = TRUE)
#
#
# sfsf <- lapply(peakToCall,function(x)read.delim(x,sep="\t",comment.char = "#",h=FALSE))
# sfsf <- lapply(sfsf,
# function(x)GRanges(x$V2,IRanges(x$V3,x$V4),strand = x$V5)
# )
# names(sfsf) <- basename(peakToCall)
# sfsf <- GRangesList(sfsf)
# rSet <- unlist(sfsf)
# nrSet <- GenomicRanges::reduce(rSet)
# mat <- do.call(cbind,lapply(sfsf,
# function(x)
# nrSet %over% x))
# mcols(nrSet) <- matAdditionally, see CLIPflexR documentation for:
countFromBed
# #load bedfiles to get read counts in peaks and bind with true/false peak matrix
# bedfile <- dir(peakDir,pattern="*.bed",full.names = TRUE)
# bedfile <- grep("unmatched", bedfile, invert= TRUE, value = TRUE)
# bedfile <- grep("33243", bedfile, invert= TRUE, value = TRUE)
# bedfile <- grep("KRG040518", bedfile, invert= TRUE, value = TRUE)
#
# countFromBed <- function(Bed,GR,notStranded=TRUE,interFeature=FALSE){
# require(rtracklayer)
# require(GenomicRanges)
# require(GenomicAlignments)
# reads <- import.bed(Bed)
# fk <- summarizeOverlaps(GR,reads,ignore.strand = notStranded,inter.feature=interFeature)
# assay(fk)
# }
#
# kks <- lapply(bedfile,countFromBed,GR=nrSet,notStranded=FALSE)
# # will get warnings here that bedfile containes sequence levels not in peak file, which is to be expected - as long as there are still sequence levels in common it's ok
#
# kksMat <- do.call(cbind,kks)
# colnames(kksMat) <- basename(bedfile)
# mcols(nrSet) <- cbind(as.data.frame(mcols(nrSet)),kksMat)
#
# #convert to dataframe for PCA/visualization
# peaksDF <- as.data.frame(nrSet)
# peaksDF$peakID <- paste0(peaksDF$seqnames, ":", peaksDF$start, "_", peaksDF$end, ":", peaksDF$strand)
#
# #write.table(peaksDF, file = "/Users/kathryn/Reprocess_all_paper_datasets/KRG_all_final_peaksDF_finalGRs_rename.txt", quote = FALSE, sep = "\t", col.names = TRUE, row.names = FALSE)output file “KRG_all_final_peaksDF_finalGRs_rename.txt” is available in my Github
input file “final_samples_KRG.txt” is available in my Github
Additionally, see CLIPflexR documentation for:
CLIP_bw2
# #cat beds to make wigs
# samples <- read.delim("/rugpfs/fs0/rice_lab/scratch/krozen/final_samples_KRG.txt", header = TRUE, sep = "\t")
# samples$group <- paste0(samples$Antibody, samples$Lysate)
#
# Dir <- "/rugpfs/fs0/rice_lab/scratch/krozen/AaegL5_mapped"
# beds <- dir(Dir, pattern= "*.bed", full.names = FALSE)
#
# samples$X <- as.character(samples$X)
# samples$X <- paste0(samples$X, "_rm5_rm3.bed")
# aegyptiAgo1 <- subset(samples$X, samples$group == "Ago1aegypti")
# aegyptiAgo2 <- subset(samples$X, samples$group == "Ago2aegypti")
# aegyptirIgG <- subset(samples$X, samples$group == "rIgGaegypti")
# aegyptimIgG <- subset(samples$X, samples$group == "mIgGaegypti")
# Aag2Ago1 <- subset(samples$X, samples$group == "Ago1Aag2")
# Aag2Ago2 <- subset(samples$X, samples$group == "Ago2Aag2")
# Aag2rIgG <- subset(samples$X, samples$group == "rIgGAag2")
# Aag2mIgG <- subset(samples$X, samples$group == "mIgGAag2")
#
# aegyptiAgo1 <- paste0(Dir, "/", aegyptiAgo1)
# aegyptiAgo2 <- paste0(Dir, "/", aegyptiAgo2)
# aegyptirIgG <- paste0(Dir, "/", aegyptirIgG)
# aegyptimIgG <- paste0(Dir, "/", aegyptimIgG)
# Aag2Ago1 <- paste0(Dir, "/", Aag2Ago1)
# Aag2Ago2 <- paste0(Dir, "/", Aag2Ago2)
# Aag2rIgG <- paste0(Dir, "/", Aag2rIgG)
# Aag2mIgG <- paste0(Dir, "/", Aag2mIgG)
#
# aegyptiAgo1bed <- lapply(aegyptiAgo1, read.delim, header = FALSE, sep = "\t")
# aegyptiAgo2bed <- lapply(aegyptiAgo2, read.delim, header = FALSE, sep = "\t")
# aegyptirIgGbed <- lapply(aegyptirIgG , read.delim, header = FALSE, sep = "\t")
# aegyptimIgGbed <- lapply(aegyptimIgG, read.delim, header = FALSE, sep = "\t")
#
# Aag2Ago1bed <- lapply(Aag2Ago1, read.delim, header = FALSE, sep = "\t")
# Aag2Ago2bed <- lapply(Aag2Ago2, read.delim, header = FALSE, sep = "\t")
# Aag2rIgGbed <- lapply(Aag2rIgG, read.delim, header = FALSE, sep = "\t")
# Aag2mIgGbed <- lapply(Aag2mIgG, read.delim, header = FALSE, sep = "\t")
#
# cat_aegyptiAgo1 <- do.call("rbind", aegyptiAgo1bed)
# cat_aegyptirIgG <- do.call("rbind", aegyptirIgGbed)
# cat_aegyptiAgo2 <- do.call("rbind", aegyptiAgo2bed)
# cat_aegyptimIgG <- do.call("rbind", aegyptimIgGbed)
# cat_Aag2Ago1 <- do.call("rbind", Aag2Ago1bed)
# cat_Aag2Ago2 <- do.call("rbind", Aag2Ago2bed)
# cat_Aag2rIgG <- do.call("rbind", Aag2rIgGbed)
# cat_Aag2mIgG <- do.call("rbind", Aag2mIgGbed)
#
#
# h <-grep("cat",names(.GlobalEnv),value=TRUE)
# l <- do.call("list",mget(h))
#
# #for(i in seq_along(l)) {
# # write.table(l[i], paste(names(l)[i], ".bed", sep = ""),
# # col.names = FALSE, row.names = FALSE, sep = "\t", quote = FALSE)
# #}
#
#
# pathToBed <- "/rugpfs/fs0/rice_lab/scratch/krozen/AaegL5_mapped"
# pathToBed <- dir(pathToBed, pattern= "cat_*", full.names = TRUE)
# print(pathToBed)
#
# #make GRanges from peaks to get reads over peak ranges
# peaksGR <- makeGRangesFromDataFrame(peaksDF,
# keep.extra.columns=FALSE,
# ignore.strand=FALSE,
# seqinfo=NULL,
# seqnames.field="seqnames",
# start.field="start",
# end.field="end",
# strand.field="strand",
# starts.in.df.are.0based=FALSE)
#
#
# makeCov <- function(pathToBed,peaks=NULL){
# test <- pathToBed %>% import.bed
# if(is.null(peaks)){
# coverage(test,weight = (1/length(test))*(10^6)) %>% export.bw(con=gsub("\\.bed","\\.bw",pathToBed))
# coverage(test[strand(test) == "+"],weight = (1/length(test))*(10^6)) %>% export.bw(con=gsub("\\.bed","Pos\\.bw",pathToBed))
# coverage(test[strand(test) == "-"],weight = (1/length(test))*(10^6)) %>% export.bw(con=gsub("\\.bed","Neg\\.bw",pathToBed))
# }else{
# coverage(test[test %over% peaks],weight = (1/length(test))*(10^6)) %>% export.bw(con=gsub("\\.bed","\\.bw",pathToBed))
# testPos <- test[strand(test) == "+"]
# testNeg <-test[strand(test) == "-"]
# coverage(testPos[testPos %over% peaks],weight = (1/length(test))*(10^6)) %>% export.bw(con=gsub("\\.bed","Pos\\.bw",pathToBed))
# coverage(testNeg[testNeg %over% peaks],weight = (1/length(test))*(10^6)) %>% export.bw(con=gsub("\\.bed","Neg\\.bw",pathToBed))
#
#
# }
# }
#
# setwd("rugpfs/fs0/rice_lab/scratch/krozen/AaegL5_mapped/bigwigs")
# bplapply(pathToBed,makeCov,peaks=peaksGR)