GENOVA相关代码

GENOVA相关代码
library(GENOVA)

centromeres = read.delim(‘C:/Users/admin/Desktop/陈*兵/GENOVA-master/data/hg19_cytobandAcen.bed’, sep = ‘\t’, h = F, stringsAsFactors = F)

head(centromeres)

Hap1_WT_20kb <- load_contacts(signal_path = ‘C:/Users/admin/Desktop/WT_DM_BL-HiC_Rep1_20000_iced.matrix’,
indices_path = ‘C:/Users/admin/Desktop/WT_DM_BL-HiC_Rep1_20000_abs.bed’,
sample_name = “DM”,
colour = “black”)

Hap1_WAPL_20kb <- load_contacts(signal_path = ‘C:/Users/admin/Desktop/WT_GM_BL-HiC_Rep1_20000_iced.matrix’,
indices_path = ‘C:/Users/admin/Desktop/WT_GM_BL-HiC_Rep1_20000_abs.bed’,
sample_name = “GM”,
colour = “red”)

p_arms <- data.frame(‘chromosome’ = centromeres[,1], ‘start’ = 0, ‘end’ = centromeres[,2])
cisChrom_out <- cis_trans( list(Hap1_WT_20kb, Hap1_WAPL_20kb),bed = p_arms)

barplot(cisChrom_out$cis,names.arg=cisChro{m}_{o}ut$sample )
abline(h = 90, col = ‘red’, lty = 3)
abline(h = 93, col = ‘red’, lty = 3)

trans.compartment.plot(Hap1_WT_20kb, remove = c(“chrM”,“chrY”))

Lets remove mitochondrial and Y-chromosomal contacts

chromosomeMatrix(Hap1_WT_1MB, remove = c(“chrM”,“chrY”))

H3K27acPeaks = read.delim(‘C:/Users/admin/Desktop/sbed/WT_DM_MyoD_ChIPseq_Rep2_trimmed_rmdup_peaks.narrowPeak’, h = F)

CS_out = compartment_score(list(Hap1_WT_20kb, Hap1_WAPL_20kb), bed = H3K27acPeaks)

visualise(CS_out, chr = “chr17”)

saddle_out = saddle(list(Hap1_WT_20kb, Hap1_WAPL_20kb), CS_discovery = CS_out, bins = 50)
visualise(saddle_out)

CSS <- quantify(saddle_out)
compared <- tidyr::spread(unique(CSS[,-c(3,4)]), key = ‘exp’, value = ‘strength’)
plot(compared$WT,compared$WAPL, xlim = c(0,4), ylim = c(0,4), pch = 20) abline(a = 0, b = 1)

cis.compartment.plot(exp1 = Hap1_WT_20kb, exp2 = Hap1_WAPL_20kb, chrom = ‘chr14’, arm = ‘q’, cs.lim = 1.75,
cut.off = 15, chip = H3K27ac_peaks)

H3K27acPeaks = read.delim(‘C:/Users/admin/Desktop/陈*兵/GSE125129data/GSE125129_H3K27ac_ChIP_wt.bw’, h = F)

RCP_out1 <- RCP(explist = list(Hap1_WT_20kb, Hap1_WAPL_20kb), chromsToUse = ‘chr1’)

Plot RCP: per-chromosome visualise(RCP_out, smooth = T)

visualise(RCP_out1, smooth = T)

CTCF = read.delim(‘C:/Users/admin/Desktop/bed/Bio760_FKDL202573273_peaks.bed’, h = F)

SMC1 = read.delim(‘C:/Users/admin/Desktop/bed/Cav58A_FKDL202573272_peaks.narrowPeak’, h = F)

RCP_out = RCP(list(Hap1_WT_20kb, Hap1_WAPL_20kb), bedlist = list(“CTCF” = CTCF, ‘Cohesin’ =SMC1), chromsToUse = ‘chr1’)
visualise(RCP_out)

H3K27acPeaks = read.delim(‘C:/Users/admin/Desktop/bed/Input_FKDL202573275_peaks.narrowPeak’, h = F)
CS_out = compartment_score(list(Hap1_WT_20kb, Hap1_WAPL_20kb), bed = H3K27acPeaks)

H3K27ac_peaks = read.delim(‘C:/Users/admin/Desktop/bed/Input_FKDL202573275_peaks.narrowPeak’, h = F)

cis.compartment.plot(exp1 = Hap1_WT_20kb, exp2 = Hap1_WAPL_20kb, chrom = ‘chr14’, arm = ‘q’, cs.lim = 1.75,
cut.off = 15, chip = H3K27ac_peaks)

WT_TADs = read.delim(‘C:/Users/admin/Desktop/杂/csd11.1_TAD.sort.bed’, h = F)
WT_Loops = read.delim(‘C:/Users/admin/Desktop/杂/postprocessed_pixels_10000.bedpe’, h = F, skip = 1)
sanborn2015_Loops = read.delim(‘C:/Users/admin/Desktop/2016hic/GSE74072_Hap1_HiCCUPS_looplist.txt.gz’)

ID <- insulation_domainogram(Hap1_WT_20kb, ‘chr7’, 25.5e6, 30e6, window_range = c(1, 101), step = 2)
hic.matrixplot(exp1 = Hap1_WT_20kb, chrom = ‘chr7’, start = 25e6, end=29e6, tads = WT_TADs,tads.type = ‘upper’,tads.colour = ‘#91cf60’, cut.off = 25, skipAnn = T)

plot(ID)

Hap1_10kb_insulation = insulation_score(list(Hap1_WT_20kb, Hap1_WAPL_20kb), window = 25)

visualise(Hap1_10kb_insulation, chr = ‘chr7’, start = 25e6, end=29e6, contrast = 1)

CTCF = read.delim(‘C:/Users/admin/Desktop/陈*兵/GSE125129data/GSE125129_CtcfPeaks_wt_merged.bed’, h = F)

SMC1 = read.delim(‘C:/Users/admin/Desktop/sbed/WT_GM_MyoD_ChIPseq_Rep2_trimmed_rmdup_peaks.narrowPeak’, h = F)

heatmap_insulation(Hap1_10kb_insulation, bed = CTCF, bed_pos = ‘center’)
TADcalls <- call_TAD_insulation(Hap1_10kb_insulation)
hic.matrixplot(exp1 = Hap1_WT_20kb, exp2 = Hap1_WAPL_20kb,
chrom = ‘chr7’, start = 24e6, end=27e6, tads = list(TADcalls$WAPL,TADcalls$WT), tads.type = list(‘lower’, ‘upper’), tads.colour = c(‘green’, ‘grey’), cut.off = 25)

ATA_Hap1_WTcalls <- ATA(list(“WT” = Hap1_WT_20kb, ‘WAPL’ = Hap1_WAPL_20kb), bed = TADcalls$WT)

visualise(ATA_Hap1_WTcalls, colour_lim = c(0,50), colour_lim_contrast = c(-5,5), focus = 1)

TAD_N_WT <- intra_inter_TAD(list(“WT” = Hap1_WT_20kb, ‘WAPL’ = Hap1_WAPL_20kb), tad_bed = TADcalls$WT, max_neighbour = 10)

visualise(TAD_N_WT, geom = ‘jitter’)
visualise(TAD_N_WT, geom = ‘violin’)

hic.matrixplot(exp1 = Hap1_WT_20kb, chrom = ‘chr7’, start = 26.75e6, end=28.5e6, loops = WT_Loops,
loops.colour = ‘#998ec3’, loops.type = ‘upper’, loops.radius = 20e3, type = ‘triangle’, chip = list(‘C:/Users/admin/Desktop/陈*兵/GSE125129data/GSE125129_H3K27ac_ChIP_wt.bw’, CTCF),symmAnn = F,cut.off = 65)

wt_Loops=WT_Loops[-1,]
bedpe = wt_Loops
bedpe[, 1]=bedpe[, 4]

WT_Loops_extended = anchors_extendedloops(Hap1_WT_20kb$IDX, res = 10e3, bedpe)

chromosomeMatrix(Hap1_WT_20kb, remove = c(“chrM”,“chrY”))

fileinput=system.file(“Test.dat.gz”,package=“CaTCH”)
library(CaTCH)
domain.call(fileinput)

C:/Users/admin/Desktop/陈兵/GSE125129_WT_HiC1/raw/40000/GSE125129_WT_HiC1_40000_abs.bed
C:/Users/admin/Desktop/陈兵/GSE125129_WT_HiC2/raw/40000/GSE125129_WT_HiC2_40000_abs.bed

C:/Users/admin/Desktop/陈兵/GSE125129_WT_HiC1/iced/40000/GSE125129_WT_HiC1_40000_iced.matrix
C:/Users/admin/Desktop/陈兵/GSE125129_WT_HiC2/iced/40000/GSE125129_WT_HiC2_40000_iced.matrix

C:/Users/admin/Desktop/陈*兵/GSM2493878_3885_9_CTCF_WT_TTAGGC_S13_L004_R1_001_peaks.narrowPeak

H3K27acPeaks = read.delim(‘C:/Users/admin/Desktop/陈*兵/GSM2493878_3885_9_CTCF_WT_TTAGGC_S13_L004_R1_001_peaks.narrowPeak.gz’, h = F)
CS_out = compartment_score(list(Hap1_WT_40kb, Hap1_WAPL_40kb), bed = H3K27acPeaks)

#!/bin/bash

#$ -cwd
#$ -S /bin/bash
#$ -l p=10
cd /data/zhangyong/yyp/yypold/hic/chenhebing/cs
source activate rna

Rscript cs.R

dir="/data/zhangyong/yyp/yypold/hic/chenhebing/cs" ##定义工作的文件夹

sedwd=(dir)
options(scipen = 20)
Args <- commandArgs(T)

library(GENOVA)

Ha40 <- load_contacts(signal_path = ‘/data/zhangyong/yyp/yypold/hic/chenhebing/cs/hic_results/matrix/GSE125129_WT_HiC1/iced/40000/GSE125129_WT_HiC1_40000_iced.matrix’,
indices_path = ‘/data/zhangyong/yyp/yypold/hic/chenhebing/cs/hic_results/matrix/GSE125129_WT_HiC1/raw/40000/GSE125129_WT_HiC1_40000_abs.bed’,
sample_name = “WT”,
colour = “black”)

Hb40 <- load_contacts(signal_path = ‘/data/zhangyong/yyp/yypold/hic/chenhebing/cs/hic_results/matrix/GSE125129_AdnpKO1/iced/40000/GSE125129_AdnpKO1_40000_iced.matrix’,
indices_path = ‘/data/zhangyong/yyp/yypold/hic/chenhebing/cs/hic_results/matrix/GSE125129_AdnpKO1/raw/40000/GSE125129_AdnpKO1_40000_abs.bed’,
sample_name = “WAPL”,
colour = “red”)

#jpeg(file = “style.jpg”) ##plot函数输出的文件名字

#cisChrom_out <- cis_trans( list(Ha40, Hb40) )
#barplot(cisChrom_out$cis,names.arg=cisChro{m}_{o}ut$sample, ylim = c(0, 100) )
#abline(h = 90, col = ‘red’, lty = 3)
#abline(h = 93, col = ‘red’, lty = 3)
#dev.off()

H3K27acPeaks = read.delim(’/data/zhangyong/yyp/yypold/hic/chenhebing/cs/GSM2493879_3885_10_CTCF_WAPL_GATCAG_S14_L004_R1_001_peaks.narrowPeak.gz’, h = F)
CS_out = compartment_score(list(Ha40, Hb40), bed = H3K27acPeaks)

pdf(file = “CSout.pdf”)

visualise(CS_out, chr = “chr17”)

dev.off()

pdf(file = “saddle.pdf”)
saddle_out = saddle(list(Ha40, Hb40), CS_discovery = CS_out, bins = 50)

visualise(saddle_out)
dev.off()

pdf(file = “CSS.pdf”)
CSS <- quantify(saddle_out)
compared <- tidyr::spread(unique(CSS[,-c(3,4)]), key = ‘exp’, value = ‘strength’)
plot(compared$WT,compared$WAPL, xlim = c(0,4), ylim = c(0,4), pch = 20)
abline(a = 0, b = 1)
dev.off()

pdf(file = “cis.compartment.pdf”)
H3K27ac_peaks = read.delim(’/data/zhangyong/yyp/yypold/hic/chenhebing/cs/GSM2493879_3885_10_CTCF_WAPL_GATCAG_S14_L004_R1_001_peaks.narrowPeak.gz’, h = F)
cis.compartment.plot(exp1 = Ha40, exp2 = Hb40, chrom = ‘chr14’, arm = ‘q’, cs.lim = 1.75, cut.off = 15, chip = H3K27ac_peaks)
dev.off()

pdf(file = “cis.compent.pdf”)
cis.compartment.plot(exp1 = Ha40, exp2 = Hb40, chrom = ‘chr20’, arm = ‘q’, cut.off = 1, obs.exp = T, chip = H3K27ac_peaks)
dev.off()

pdf(file = “hic.matrixplot.pdf”)
hic.matrixplot(exp1 = Ha40, chrom = ‘chr7’, start = 25e6, end=30e6, cut.off = 50)
dev.off()

pdf(file = “hic.twomatrixplot.pdf”)
hic.matrixplot(exp1 = Ha40, exp2 = Hb40, coplot = ‘diff’, chrom = ‘chr7’, start = 25e6, end=30e6, cut.off = 25)
dev.off()

CTCF = read.delim(’/data/zhangyong/yyp/yypold/hic/chenhebing/cs/Bio760_FKDL202573273_peaks.bed’, h = F)

SMC1 = read.delim(’/data/zhangyong/yyp/yypold/hic/chenhebing/cs/GSM2493879_3885_10_CTCF_WAPL_GATCAG_S14_L004_R1_001_peaks.narrowPeak.gz’, h = F)

WT_TADs = read.delim(’/data/zhangyong/yyp/yypold/hic/chenhebing/cs/csd11.1_TAD.sort.bed’, h = F)
WT_Loops = read.delim(’/data/zhangyong/yyp/yypold/hic/chenhebing/cs/postprocessed_pixels_10000.bedpe’, h = F, skip = 1)

sanborn2015_Loops = read.delim(’/data/zhangyong/yyp/yypold/hic/chenhebing/cs/GSE63525_GM12878_replicate_HiCCUPS_looplist.txt.gz’)

pdf(file = “hic.matrixplot.pdf”)
hic.matrixplot(exp1 = Hap1_WT_10kb, chrom = ‘chr7’, start = 25e6, end=30e6, loops = WT_Loops, loops.colour = ‘blue’, loops.type = ‘upper’, loops.radius = 20e3, tads = WT_TADs, tads.type = ‘lower’,tads.colour = ‘#ffd92f’, cut.off = 25)

dev.off()

pdf(file = “hic.2matrixplot.pdf”)
hic.matrixplot(exp1 = Hap1_WT_10kb, chrom = ‘chr7’, start = 26.75e6, end=28.5e6, loops = WT_Loops, loops.colour = ‘#998ec3’, loops.type = ‘upper’, loops.radius = 20e3, type = ‘triangle’, chip = list(’/data/zhangyong/yyp/yypold/hic/chenhebing/cs/GSE125129_H3K27ac_ChIP_wt.bw’,CTCF),symmAnn = F, cut.off = 65)
dev.off()

pdf(file = “hic.3matrixplot.pdf”)
ID <- insulation_domainogram(Ha40, ‘chr7’, 25.5e6, 30e6, window_range = c(1, 101), step = 2)

hic.matrixplot(exp1 = Hap1_WT_10kb, chrom = ‘chr7’, start = 25e6, end=29e6, tads = WT_TADs, tads.type = ‘upper’, tads.colour = ‘#91cf60’, cut.off = 25,skipAnn = T)
plot(ID)

dev.off()

pdf(file = “Hap1_10kb_insulation.pdf”)
Hap1_10kb_insulation = insulation_score(list(Ha40, Hb40), window = 25)
visualise(Hap1_10kb_insulation, chr = ‘chr7’, start = 25e6, end=29e6, contrast = 1)
dev.off()

pdf(file = “heatmapHap1_10kb_insulation.pdf”)

heatmap_insulation(Hap1_10kb_insulation, bed = CTCF, bed_pos = ‘center’)
dev.off()

pdf(file = “TADcalls.pdf”)
TADcalls <- call_TAD_insulation(Hap1_10kb_insulation)
hic.matrixplot(exp1 = Ha40, exp2 = Hb40,chrom = ‘chr7’, start = 24e6, end=27e6, tads = list(TADcalls$WAPL,TADcalls$WT), tads.type = list(‘lower’, ‘upper’), tads.colour = c(‘green’, ‘grey’), cut.off = 25)
dev.off()

WT_Loops_extended = anchors_extendedloops(Ha40$IDX, res = 10e3, bedpe = WT_Loops)

devtools::install_github(“robinweide/GENOVA”, ref = ‘dev’)

library(GENOVA)

Hap1_WT_20kb <- load_contacts(signal_path = ‘/data/zhangyong/yyp/yypold/hic/hicdata/outputrna/hic_results/matrix/WT_DM_BL-HiC_Rep1/iced/20000/WT_DM_BL-HiC_Rep1_20000_iced.matrix’,
indices_path = ‘/data/zhangyong/yyp/yypold/hic/hicdata/outputrna/hic_results/matrix/WT_DM_BL-HiC_Rep1/raw/20000/WT_DM_BL-HiC_Rep1_20000_abs.bed’,
sample_name = “WT”,
colour = “black”)

Hap1_WAPL_20kb <- load_contacts(signal_path = ‘/data/zhangyong/yyp/yypold/hic/hicdata/outputrna/hic_results/matrix/WT_DM_BL-HiC_Rep2/iced/20000/WT_DM_BL-HiC_Rep2_20000_iced.matrix’,
indices_path = ‘/data/zhangyong/yyp/yypold/hic/hicdata/outputrna/hic_results/matrix/WT_DM_BL-HiC_Rep2/raw/20000/WT_DM_BL-HiC_Rep2_20000_abs.bed’,
sample_name = “WAPL”,
colour = “red”)

Hap1_SCC4_20kb <- load_contacts(signal_path = ‘/data/zhangyong/yyp/yypold/hic/hicdata/outputrna/hic_results/matrix/WT_DM_BL-HiC_Rep3/iced/20000/WT_DM_BL-HiC_Rep3_20000_iced.matrix’,
indices_path = ‘/data/zhangyong/yyp/yypold/hic/hicdata/outputrna/hic_results/matrix/WT_DM_BL-HiC_Rep3/raw/20000/WT_DM_BL-HiC_Rep3_20000_abs.bed’,
sample_name = “SCC4”,
colour = “green”)

Hap1_WT_40kb <- load_contacts(signal_path = ‘/data/zhangyong/yyp/yypold/hic/hicdata/outputrna/hic_results/matrix/WT_DM_BL-HiC_Rep4/iced/40000/WT_DM_BL-HiC_Rep4_40000_iced.matrix’,
indices_path = ‘/data/zhangyong/yyp/yypold/hic/hicdata/outputrna/hic_results/matrix/WT_DM_BL-HiC_Rep4/raw/40000/WT_DM_BL-HiC_Rep4_40000_abs.bed’,
sample_name = “WT”,
colour = “black”)

Hap1_WAPL_40kb <- load_contacts(signal_path = ‘/data/zhangyong/yyp/yypold/hic/hicdata/outputrna/hic_results/matrix/WT_DM_BL-HiC_Rep5/iced/40000/WT_DM_BL-HiC_Rep5_40000_iced.matrix’,
indices_path = ‘/data/zhangyong/yyp/yypold/hic/hicdata/outputrna/hic_results/matrix/WT_DM_BL-HiC_Rep5/raw/40000/WT_DM_BL-HiC_Rep5_40000_abs.bed’,
sample_name = “WAPL”,
colour = “red”)

Hap1_WT_1MB <- load_contacts(signal_path = ‘/data/zhangyong/yyp/yypold/hic/hicdata/outputrna/hic_results/matrix/WT_DM_BL-HiC_Rep6/iced/1000000/WT_DM_BL-HiC_Rep6_1000000_iced.matrix’,
indices_path = ‘/data/zhangyong/yyp/yypold/hic/hicdata/outputrna/hic_results/matrix/WT_DM_BL-HiC_Rep6/raw/20000/WT_DM_BL-HiC_Rep6_20000_abs.bed’,
sample_name = “WT”, centromeres = centromeres,
colour = “black”)
Hap1_WAPL_1MB <- load_contacts(signal_path = ‘/data/zhangyong/yyp/yypold/hic/hicdata/outputrna/hic_results/matrix/WT_DM_BL-HiC_Rep7/iced/1000000/WT_DM_BL-HiC_Rep7_1000000_iced.matrix’,
indices_path = ‘/data/zhangyong/yyp/yypold/hic/hicdata/outputrna/hic_results/matrix/WT_DM_BL-HiC_Rep7/raw/1000000/WT_DM_BL-HiC_Rep7_1000000_abs.bed’,
sample_name = “WAPL”,
centromeres = centromeres,
colour = “red”)

Hap1_WT_10kb_juicer <- load_contacts(signal_path = ‘/data/zhangyong/yyp/yypold/hic/hicdata/outputrna/hicfile_WT_DM_BL-HiC_Rep1/WT_DM_BL-HiC_Rep1.hic’,
sample_name = “WT”,
resolution = 10e3,
balancing = ‘KR’, # this is the default
colour = “black”)

Hap1_WAPL_10kb_juicer <- load_contacts(signal_path = ‘/data/zhangyong/yyp/yypold/hic/hicdata/outputrna/hicfile_WT_DM_BL-HiC_Rep2/WT_DM_BL-HiC_Rep2.hic’,
sample_name = “WT”,
resolution = 10e3,
balancing = ‘KR’, # this is the default
colour = “black”)

RCP_out <- RCP(explist = list(Hap1_WT_40kb, Hap1_WAPL_40kb),
chromsToUse = ‘chr1’)

Plot RCP: per-chromosome

visualise(RCP_out,
smooth = T)

devtools::install_github(“robinweide/GENOVA”, ref = ‘dev’)

library(GENOVA)
library(strawr)

Hap1_WT_10kb_juicer <- load_contacts(signal_path = ‘/data/zhangyong/yyp/yypold/hic/hicdata/outputrna/hicfile_WT_DM_BL-HiC_Rep1/WT_DM_BL-HiC_Rep1.hic’,
sample_name = “WT”,
resolution = 10e3,
balancing = ‘KR’, # this is the default
colour = “black”)

Hap1_WAPL_10kb_juicer <- load_contacts(signal_path = ‘/data/zhangyong/yyp/yypold/hic/hicdata/outputrna/hicfile_WT_DM_BL-HiC_Rep2/WT_DM_BL-HiC_Rep2.hic’,
sample_name = “WT”,
resolution = 10e3,
balancing = ‘KR’, # this is the default
colour = “black”)

RCP_out <- RCP(explist = list(Hap1_WT_10kb_juicer, Hap1_WAPL_10kb_juicer),
chromsToUse = ‘chr1’)

Plot RCP: per-chromosome

visualise(RCP_out,
smooth = T)