Delete BagFoot_analysis directory

# Bivariate-Genomic-Footprinting

Code for footprinting and motif-flanking accessibility using BaGFoot (Baek et al)

aux.R

bagfoot.R

Baek S, et al (2017) ***Genomic Footprinting Detects Changes in Transcription Factor Activity*** ([Cell Rep 19, 1710-1722](https://pubmed.ncbi.nlm.nih.gov/28538187/)).

#install.packages("hash")

#install.packages("digest")

#install.packages("data.table")

#install.packages("Cairo")

#install.packages("aplpack")

#source("https://bioconductor.org/biocLite.R")

#biocLite("BSgenome.Hsapiens.UCSC.hg19")

#biocLite("BSgenome.Mmusculus.UCSC.mm9")

setwd("./bagfoot1")

library(Rsamtools)

library(hash)

library(data.table)

library(digest)

library(Cairo)

library(aplpack)

library(BSgenome.Hsapiens.UCSC.hg19)

library(BSgenome.Mmusculus.UCSC.mm9)

library(BSgenome.Mmusculus.UCSC.mm10)

library(bagfoot)

source("aux.R")

bamfile1= 'LN.bam';  

cc1 = countReadsBAM(bamfile1, refgenome = "mm10");  # counts the number of cuts in the sequence file.

cc1

cutcountfile1 = makeCutCountBAM(bamfile1, refgenome = "mm10");   # generate a BedGraph file with ATAC Cleavages counts

bamfile2= 'DM.bam';

cc2 = countReadsBAM(bamfile2, refgenome = "mm10");  

cc2

cutcountfile2 = makeCutCountBAM(bamfile2, refgenome = "mm10"); 

hotspotfile1 = 'A1_peaks.csv';   # A peak sites file is comma-separated and must has headers with chromosome, start, end (1-based)

hotspotfile2 = 'A1_peaks.csv';   # this can be different, or same, set of peaks

combinedhotspotfile = combineTwoHotspots( hotspotfile1, hotspotfile2, 'LN', 'DM');

#for bedtools multicov

pooledBed <- read.csv('pooled_LN_DM_hotspot.csv', header=TRUE)

write.table(x=pooledBed[,2:4], file = "pooledBed.bed", append = FALSE, quote = FALSE, sep = "\t", row.names =FALSE, col.names = FALSE)

system(' bedtools multicov -bams LN.bam -bed pooledBed.bed > NcountsFile1.bed')

system(' bedtools multicov -bams DM.bam -bed pooledBed.bed > NcountsFile2.bed')

NcountsFile1 <- sum( read.csv('NcountsFile1.bed', head=FALSE , sep="\t")[,4] ,na.rm=TRUE )

NcountsFile2 <- sum( read.csv('NcountsFile2.bed', head=FALSE , sep="\t")[,4] ,na.rm=TRUE )

NcountsFile1

NcountsFile2

rm(pooledBed)

system(' rm pooledBed.bed ' )

tabNoMappability = MakeBiasCorrectionTableBAM(  # This function call generates a hexamer bias frequency table from the given bamfile without mappability assumption

	bamfile='all.bam',

	outfile="Hexamer_all_mm10_withoutMap.txt",

	refgenome="mm10", 

	np=6, 

	mapdir='',    # if mappability file directory is set to empty, mappability is not used. 

	atac=TRUE     # Set TRUE for ATAC data. The default value is FALSE

	);

# run "BaGFoot"

LN <- CutCount(file = "./bagfoot1/LN_AC_cutcount.bgr.gz" , count = NcountsFile1, name = "LN" ); 

DM <- CutCount(file = "./bagfoot1/DM_AC_cutcount.bgr.gz" , count = NcountsFile2, name = "DM" );

motifdb_mm10 <- MotifDB(motiflistfile = 'motiflist_mm10_cisbp_filtered.txt',   # This file contains the list of FIMO outputs file

			directory='/mm10_cisbp_fimo/');	                       # Directory of FIMO outputs

gfootoption_ctrl_treat <- GFootOption(biasfile = "Hexamer_all_mm10_withoutMap.txt",                   # Hexamer bias frequency table (see "bigfoot_prep_example.R")

					hexamer_pattern= "./bagfoot1/nuccode_mm10_6mer_{chr}.dat");   # nuc. code files generated by MakeBiasCorrectionTableBAM 

# Definition of BiGFoot analysis 

GFoot_ctrl_treat= GFoot(control = LN,      ##  Control Data defined above

			treatment = DM,    ##  Treatment Data defined above

			sitefile = "./bagfoot1/pooled_LN_DM_hotspot.csv",  ##  Sites file: Peak File.  This file must have "chr", "st", "ed" in the header.

			motifDB = motifdb_mm10,                ##  Information about FIMO motifs

			gfootoption= gfootoption_ctrl_treat,   ##  Options as defined above

			outputdir = "OUTPUT_LN_vs_DM",  ## OUTPUT directory

			cachedir = "CACHE_LN_vs_DM");   ## CACHE directory 

GFoot_ctrl_treat <- run(GFoot_ctrl_treat , graphout=T,yrange=c(-2.5,1.5), mc.cores = 24 ) 

dat <- read.table('LN_DM_On_pooled_LN_DM_hotspot_footprint_depth_table.csv');

gen_bagplot_chisq(dat, dataname1='LN', dataname2='DM', factor=2.5);   # generates a bagplot from the calculated footprinting depths

#!/bin/bash

Rscript --vanilla bagfoot.R