Delete pCHiC directory

# Chromatin interaction anlaysis on pCHiC data

## pCHiC data mapping:

hicup.sh

runHiCU_pCHiCmapping.sh

## Chromatin compartment analysis:

runHOMER_compartment.sh

## Interaction calling:

runChicago_interactions_calling.sh

## Interactions viewed as 4C data

plotBaits.R

## Differential anlaysis at individual promoters:

runEdgeR_pCHiC_promoter_reads.R

## Interactions rewiring analysis:

pCHiC_rewire_analysis.R

#!/bin/bash

GSE=$1

GSM=$2

fastq=$3

cd $GSE/$GSM/

file1="$GSM.r_1.fastq"

file2="$GSM.r_2.fastq"

digestFile="/serenity/data/reference-genomes/Digest_mm10_HindIII_None_15-51-59_29-04-2016.txt"

path="/usr/bin/bowtie2"

idx="/serenity/data/reference-genomes/mm10"

hicup --bowtie2 $path  --index $idx --digest $digestFile --format Sanger --longest 800 --shortest 150 --threads 12 $file1 $file2

setwd("~/Documents/Chicago_new/Differential.analysis/rewire_analysis")

All.score <- read.delim("R1-2_chicago.scores.matrix.txt", sep="\t",header=TRUE,row.names = 1)

### high means gained interations, low means lost interactions

############### WT vs NPM1

WT.pos.NPM1.neg <- read.delim("WT.pos.NPM1.neg.txt", sep="\t",header=FALSE)

WT.neg.NPM1.pos <- read.delim("WT.neg.NPM1.pos.txt", sep="\t",header=FALSE)

WT.NPM1.all <- rbind(WT.pos.NPM1.neg,WT.neg.NPM1.pos)

WT.NPM1.all.score <- All.score[row.names(All.score) %in% WT.NPM1.all[,1],1:2]

a = as.integer(nrow(WT.NPM1.all.score)*0.25) 

b = as.integer(nrow(WT.NPM1.all.score)*0.5) 

c = as.integer(nrow(WT.NPM1.all.score)*0.75)

WT.NPM1.tmp1 <- WT.NPM1.all.score[order(-WT.NPM1.all.score[,1]),] 

WT.NPM1.tmp1$WT.group <- c(rep("A",a),rep("B",b-a),rep("C",c-b),rep("D",nrow(WT.NPM1.tmp1)-c))

WT.NPM1.tmp2 <- WT.NPM1.tmp1[order(-WT.NPM1.tmp1[,2]),] 

WT.NPM1.tmp2$NPM1.group <- c(rep("A",a),rep("B",b-a),rep("C",c-b),rep("D",nrow(WT.NPM1.tmp2)-c))

WT.high.NPM1.low <- subset(WT.NPM1.tmp2,WT.group == "A" & NPM1.group == "D")

WT.low.NPM1.high <- subset(WT.NPM1.tmp2,WT.group == "D" & NPM1.group == "A")

write.table(WT.high.NPM1.low,file="WT.high.NPM1.low.txt",row.names=TRUE,quote = FALSE, sep = "\t",col.names = TRUE)

write.table(WT.low.NPM1.high,file="WT.low.NPM1.high.txt",row.names=TRUE,quote = FALSE, sep = "\t",col.names = TRUE)

############### WT vs FLT3

WT.pos.FLT3.neg <- read.delim("WT.pos.FLT3.neg.txt", sep="\t",header=FALSE)

WT.neg.FLT3.pos <- read.delim("WT.neg.FLT3.pos.txt", sep="\t",header=FALSE)

WT.FLT3.all <- rbind(WT.pos.FLT3.neg,WT.neg.FLT3.pos)

WT.FLT3.all.score <- All.score[row.names(All.score) %in% WT.FLT3.all[,1],c(1,3)] ## select specific columns

a = as.integer(nrow(WT.FLT3.all.score)*0.25) 

b = as.integer(nrow(WT.FLT3.all.score)*0.5) 

c = as.integer(nrow(WT.FLT3.all.score)*0.75)

WT.FLT3.tmp1 <- WT.FLT3.all.score[order(-WT.FLT3.all.score[,1]),] 

WT.FLT3.tmp1$WT.group <- c(rep("A",a),rep("B",b-a),rep("C",c-b),rep("D",nrow(WT.FLT3.tmp1)-c))

WT.FLT3.tmp2 <- WT.FLT3.tmp1[order(-WT.FLT3.tmp1[,2]),] 

WT.FLT3.tmp2$FLT3.group <- c(rep("A",a),rep("B",b-a),rep("C",c-b),rep("D",nrow(WT.FLT3.tmp2)-c))

WT.high.FLT3.low <- subset(WT.FLT3.tmp2,WT.group == "A" & FLT3.group == "D")

WT.low.FLT3.high <- subset(WT.FLT3.tmp2,WT.group == "D" & FLT3.group == "A")

write.table(WT.high.FLT3.low,file="WT.high.FLT3.low.txt",row.names=TRUE,quote = FALSE, sep = "\t",col.names = TRUE)

write.table(WT.low.FLT3.high,file="WT.low.FLT3.high.txt",row.names=TRUE,quote = FALSE, sep = "\t",col.names = TRUE)

############### WT vs DM

WT.pos.DM.neg <- read.delim("WT.pos.DM.neg.txt", sep="\t",header=FALSE)

WT.neg.DM.pos <- read.delim("WT.neg.DM.pos.txt", sep="\t",header=FALSE)

WT.DM.all <- rbind(WT.pos.DM.neg,WT.neg.DM.pos)

WT.DM.all.score <- All.score[row.names(All.score) %in% WT.DM.all[,1],c(1,4)] 

a = as.integer(nrow(WT.DM.all.score)*0.25) 

b = as.integer(nrow(WT.DM.all.score)*0.5) 

c = as.integer(nrow(WT.DM.all.score)*0.75)

WT.DM.tmp1 <- WT.DM.all.score[order(-WT.DM.all.score[,1]),] 

WT.DM.tmp1$WT.group <- c(rep("A",a),rep("B",b-a),rep("C",c-b),rep("D",nrow(WT.DM.tmp1)-c))

WT.DM.tmp2 <- WT.DM.tmp1[order(-WT.DM.tmp1[,2]),] 

WT.DM.tmp2$DM.group <- c(rep("A",a),rep("B",b-a),rep("C",c-b),rep("D",nrow(WT.DM.tmp2)-c))

WT.high.DM.low <- subset(WT.DM.tmp2,WT.group == "A" & DM.group == "D")

WT.low.DM.high <- subset(WT.DM.tmp2,WT.group == "D" & DM.group == "A")

write.table(WT.high.DM.low,file="WT.high.DM.low.txt",row.names=TRUE,quote = FALSE, sep = "\t",col.names = TRUE)

write.table(WT.low.DM.high,file="WT.low.DM.high.txt",row.names=TRUE,quote = FALSE, sep = "\t",col.names = TRUE)

library(devtools)

library(Chicago)

library(PCHiCdata)

########## plotBaits ##########

library(Chicago)

setwd("~/analyses/Chicago/enhancer_project")

WT.run <- readRDS("WT.run.rds")

NPM1.run <- readRDS("NPM1.run.rds")

FLT3.run <- readRDS("FLT3.run.rds")

DM.run <- readRDS("DM.run.rds")

## Irf8 729806

plotBaits(WT.run, n=1, Ncol = "N", baits = 729806,ylim=c(0,150),plotBaitNames = TRUE, plotBprof = TRUE, outfile = "Irf8_4C_norm_WT.pdf", width = 7, height = 3.5, removeBait2bait = FALSE, maxD = 0.1e6, bgCol = "black", plevel1 = 5, plevel2 = 3, lev1Col = "red", lev2Col = "blue",bgPch = 20, lev1Pch = 18, lev2Pch = 18)

plotBaits(NPM1.run, n=1, Ncol = "N", baits = 729806,ylim=c(0,170),plotBaitNames = TRUE, plotBprof = TRUE, outfile = "Irf8_4C_norm_NPM1.pdf", width = 7, height = 3.5, removeBait2bait = FALSE, maxD = 0.1e6, bgCol = "black", plevel1 = 5, plevel2 = 3, lev1Col = "red", lev2Col = "blue",bgPch = 20, lev1Pch = 18, lev2Pch = 18)

plotBaits(FLT3.run, n=1, Ncol = "N", baits = 729806,ylim=c(0,155),plotBaitNames = TRUE, plotBprof = TRUE, outfile = "Irf8_4C_norm_FLT3.pdf", width = 7, height = 3.5, removeBait2bait = FALSE, maxD = 0.1e6, bgCol = "black", plevel1 = 5, plevel2 = 3, lev1Col = "red", lev2Col = "blue",bgPch = 20, lev1Pch = 18, lev2Pch = 18)

plotBaits(DM.run, n=1, Ncol = "N", baits = 729806,ylim=c(0,120),plotBaitNames = TRUE, plotBprof = TRUE, outfile = "Irf8_4C_norm_DM.pdf", width = 7, height = 3.5, removeBait2bait = FALSE, maxD = 0.1e6, bgCol = "black", plevel1 = 5, plevel2 = 3, lev1Col = "red", lev2Col = "blue",bgPch = 20, lev1Pch = 18, lev2Pch = 18)

### Spi1 425092

plotBaits(WT.run, n=1, Ncol = "N", baits = 425092,ylim=c(0,160),plotBaitNames = TRUE, plotBprof = TRUE, outfile = "Spi1_4C_norm_WT.pdf", width = 7, height = 3.5, removeBait2bait = FALSE, maxD = 0.03e6, bgCol = "black", plevel1 = 5, plevel2 = 3, lev1Col = "red", lev2Col = "blue",bgPch = 20, lev1Pch = 18, lev2Pch = 18)

plotBaits(NPM1.run, n=1, Ncol = "N", baits = 425092,ylim=c(0,200),plotBaitNames = TRUE, plotBprof = TRUE, outfile = "Spi1_4C_norm_NPM1.pdf", width = 7, height = 3.5, removeBait2bait = FALSE, maxD = 0.03e6, bgCol = "black", plevel1 = 5, plevel2 = 3, lev1Col = "red", lev2Col = "blue",bgPch = 20, lev1Pch = 18, lev2Pch = 18)

plotBaits(FLT3.run, n=1, Ncol = "N", baits = 425092,ylim=c(0,200),plotBaitNames = TRUE, plotBprof = TRUE, outfile = "Spi1_4C_norm_FLT3.pdf", width = 7, height = 3.5, removeBait2bait = FALSE, maxD = 0.03e6, bgCol = "black", plevel1 = 5, plevel2 = 3, lev1Col = "red", lev2Col = "blue",bgPch = 20, lev1Pch = 18, lev2Pch = 18)

plotBaits(DM.run, n=1, Ncol = "N", baits = 425092,ylim=c(0,200),plotBaitNames = TRUE, plotBprof = TRUE, outfile = "Spi1_4C_norm_DM.pdf", width = 7, height = 3.5, removeBait2bait = FALSE, maxD = 0.03e6, bgCol = "black", plevel1 = 5, plevel2 = 3, lev1Col = "red", lev2Col = "blue",bgPch = 20, lev1Pch = 18, lev2Pch = 18)

### Hoxa9 617336

plotBaits(WT.run, n=1, Ncol = "N", baits = 617336, ylim=c(0,100), plotBaitNames = TRUE, plotBprof = TRUE, outfile = "Hoxa9_4C_norm_WT.pdf", width = 7, height = 3.5, removeBait2bait = FALSE, maxD = 1.2e6, bgCol = "black", plevel1 = 5, plevel2 = 3, lev1Col = "red", lev2Col = "blue",bgPch = 20, lev1Pch = 18, lev2Pch = 18)

plotBaits(NPM1.run, n=1, Ncol = "N", baits = 617336, ylim=c(0,100), plotBaitNames = TRUE, plotBprof = TRUE, outfile = "Hoxa9_4C_norm_NPM1.pdf", width = 7, height = 3.5, removeBait2bait = FALSE, maxD = 1.2e6, bgCol = "black", plevel1 = 5, plevel2 = 3, lev1Col = "red", lev2Col = "blue",bgPch = 20, lev1Pch = 18, lev2Pch = 18)

plotBaits(FLT3.run, n=1, Ncol = "N", baits = 617336, ylim=c(0,100),plotBaitNames = TRUE, plotBprof = TRUE, outfile = "Hoxa9_4C_norm_FLT3.pdf", width = 7, height = 3.5, removeBait2bait = FALSE, maxD = 1.2e6, bgCol = "black", plevel1 = 5, plevel2 = 3, lev1Col = "red", lev2Col = "blue",bgPch = 20, lev1Pch = 18, lev2Pch = 18)

plotBaits(DM.run, n=1, Ncol = "N", baits = 617336, ylim=c(0,100), plotBaitNames = TRUE, plotBprof = TRUE, outfile = "Hoxa9_4C_norm_DM.pdf", width = 7, height = 3.5, removeBait2bait = FALSE, maxD = 1.2e6, bgCol = "black", plevel1 = 5, plevel2 = 3, lev1Col = "red", lev2Col = "blue",bgPch = 20, lev1Pch = 18, lev2Pch = 18)

### Hoxa10 617342

plotBaits(WT.run, n=1, Ncol = "N", baits = 617342, ylim=c(0,100),plotBaitNames = TRUE, plotBprof = TRUE, outfile = "Hoxa10_4C_norm_WT.pdf", width = 7, height = 3.5, removeBait2bait = FALSE, maxD = 1.2e6, bgCol = "black", plevel1 = 5, plevel2 = 3, lev1Col = "red", lev2Col = "blue",bgPch = 20, lev1Pch = 18, lev2Pch = 18)

plotBaits(NPM1.run, n=1, Ncol = "N", baits = 617342, ylim=c(0,100),plotBaitNames = TRUE, plotBprof = TRUE, outfile = "Hoxa10_4C_norm_NPM1.pdf", width = 7, height = 3.5, removeBait2bait = FALSE, maxD = 1.2e6, bgCol = "black", plevel1 = 5, plevel2 = 3, lev1Col = "red", lev2Col = "blue",bgPch = 20, lev1Pch = 18, lev2Pch = 18)

plotBaits(FLT3.run, n=1, Ncol = "N", baits = 617342, ylim=c(0,100),plotBaitNames = TRUE, plotBprof = TRUE, outfile = "Hoxa10_4C_norm_FLT3.pdf", width = 7, height = 3.5, removeBait2bait = FALSE, maxD = 1.2e6, bgCol = "black", plevel1 = 5, plevel2 = 3, lev1Col = "red", lev2Col = "blue",bgPch = 20, lev1Pch = 18, lev2Pch = 18)

plotBaits(DM.run, n=1, Ncol = "N", baits = 617342, ylim=c(0,100),plotBaitNames = TRUE, plotBprof = TRUE, outfile = "Hoxa10_4C_norm_DM.pdf", width = 7, height = 3.5, removeBait2bait = FALSE, maxD = 1.2e6, bgCol = "black", plevel1 = 5, plevel2 = 3, lev1Col = "red", lev2Col = "blue",bgPch = 20, lev1Pch = 18, lev2Pch = 18)

### Gfi1b 404376

plotBaits(WT.run, n=1, Ncol = "N", baits = 404376, ylim=c(0,115), plotBaitNames = TRUE, plotBprof = TRUE, outfile = "Gfi1b_4C_norm_WT.pdf", width = 7, height = 3.5, removeBait2bait = FALSE, maxD = 0.04e6, bgCol = "black", plevel1 = 5, plevel2 = 3, lev1Col = "red", lev2Col = "blue",bgPch = 20, lev1Pch = 18, lev2Pch = 18)

plotBaits(NPM1.run, n=1, Ncol = "N", baits = 404376, ylim=c(0,160), plotBaitNames = TRUE, plotBprof = TRUE, outfile = "Gfi1b_4C_norm_NPM1.pdf", width = 7, height = 3.5, removeBait2bait = FALSE, maxD = 0.04e6, bgCol = "black", plevel1 = 5, plevel2 = 3, lev1Col = "red", lev2Col = "blue",bgPch = 20, lev1Pch = 18, lev2Pch = 18)

plotBaits(FLT3.run, n=1, Ncol = "N", baits = 404376, ylim=c(0,160), plotBaitNames = TRUE, plotBprof = TRUE, outfile = "Gfi1b_4C_norm_FLT3.pdf", width = 7, height = 3.5, removeBait2bait = FALSE, maxD = 0.04e6, bgCol = "black", plevel1 = 5, plevel2 = 3, lev1Col = "red", lev2Col = "blue",bgPch = 20, lev1Pch = 18, lev2Pch = 18)

plotBaits(DM.run, n=1, Ncol = "N", baits = 404376, ylim=c(0,115), plotBaitNames = TRUE, plotBprof = TRUE, outfile = "Gfi1b_4C_norm_DM.pdf", width = 7, height = 3.5, removeBait2bait = FALSE, maxD = 0.04e6, bgCol = "black", plevel1 = 5, plevel2 = 3, lev1Col = "red", lev2Col = "blue",bgPch = 20, lev1Pch = 18, lev2Pch = 18)

#!/bin/bash

# convert bam to chicago chinput files

bam2chicago.sh WT.CHiC.R1.bam CHiC.mm10.baitmap Digest.mm10.rmap WT.CHiC.R1 nodelete

bam2chicago.sh WT.CHiC.R2.bam CHiC.mm10.baitmap Digest.mm10.rmap WT.CHiC.R2 nodelete

bam2chicago.sh DM.CHiC.R1.bam CHiC.mm10.baitmap Digest.mm10.rmap DM.CHiC.R1 nodelete

bam2chicago.sh DM.CHiC.R2.bam CHiC.mm10.baitmap Digest.mm10.rmap DM.CHiC.R2 nodelete

bam2chicago.sh FLT3.CHiC.R1.bam CHiC.mm10.baitmap Digest.mm10.rmap FLT3.CHiC.R1 nodelete

bam2chicago.sh FLT3.CHiC.R2.bam CHiC.mm10.baitmap Digest.mm10.rmap FLT3.CHiC.R2 nodelete

bam2chicago.sh NPM1.CHiC.R1.bam CHiC.mm10.baitmap Digest.mm10.rmap NPM1.CHiC.R1 nodelete

bam2chicago.sh NPM1.CHiC.R2.bam CHiC.mm10.baitmap Digest.mm10.rmap NPM1.CHiC.R2 nodelete

# Call interactions on chinput files

Rscript runChicago.R --design-dir .-design-dir WT.CHiC.R1.chinput,WT.CHiC.R2.chinput WT.CHiC.R1-2

Rscript runChicago.R --design-dir .-design-dir DM.CHiC.R1.chinput,DM.CHiC.R2.chinput DM.CHiC.R1-2

Rscript runChicago.R --design-dir .-design-dir NPM1.CHiC.R1.chinput,NPM1.CHiC.R2.chinput NPM1.CHiC.R1-2

Rscript runChicago.R --design-dir .-design-dir FLT3.CHiC.R1.chinput,FLT3.CHiC.R2.chinput FLT3.CHiC.R1-2