Library and data import
Building objects for future work
Loops from Bonev et al
Analysis of transcriptome with and without CTCF (ES and NS cells)
- ES–>NS validation
- CTCF targets cell type specific genes
- Aldh1a3 expression plot
- Aldh1a3 CTCF-BS KOs
Expression of up and down regulated genes
Chromatin signature with and without CTCF - ATAC-seq and ChIP-seq data
- Peak number changes and deregulation upon CTCF loss
Peaks in ES and NS cells
CTCF, gene expression, and loop gain
Insulators in the NS cells
PLA
PLA : Ctcf-Nono
PLA : FUS-Ddx5
PANTR1 PLA Ddx5-CTCF
PANTR1 PLA FUS-CTCF
qRT-PCR
CTCF clustering - AiryScan on pre-extracted nuclei
CTCF peaks
- CTCF AUC
- Motif search
- Factors enriched at CTCF sites losing CTCF signal in Ddx5-/- NS cells
- Factors enriched at CTCF sites gaining CTCF signal in Ddx5-/- NS cells
- CTCF motif strength
- G4
Load loops
CTCF and G4 at loop anchors
Ddx5 and CTCF at loop anchors
Deep coverage data preparations
Distance decline plot
Loop analysis
Loop signal
Lost and gained loops in the wt and Ddx5-/- cells
IPF loop signal in the intervals of genomic distances
Examples of loops lost in the NS cells
- APA - general
4F11 - global analysis of architectural loops
- Loops
Can we link loop loss to Ddx5 effect on CTCF?
- dTAG13 and chromatin interactions
Pantr1 - general Hi-C
CTCF in Pantr1-/-
APA in Pantr1 on loops lost in Ddx5
Hi-C: Insulation
Published data - insulation
Pantr1 Insulation
SICAP-ChIP
Session Info

Library and data import

data_directory='/Volumes/T7/Dehingia_et_al_2025/data/'
chipseq_directory='/Volumes/T7/Dehingia_et_al_2025/data/ChIP-and-ATAC-seq/'
rnaseq_directory='/Volumes/T7/Dehingia_et_al_2025/data/RNA-seq/'
hic_directory='/Volumes/T7/Dehingia_et_al_2025/data/Hi-C/'
outputs_directory = '/Volumes/T7/Dehingia_et_al_2025/data/outputs/'
objects_directory = '/Volumes/T7/Dehingia_et_al_2025/data/objects/'
scripts_directory = '/Volumes/T7/Dehingia_et_al_2025/scripts/'
source_data_directory = '/Volumes/T7/Dehingia_et_al_2025/source_files/'

library(rtracklayer)

## Warning: pakiet 'rtracklayer' został zbudowany w wersji R 4.1.1

## Warning: pakiet 'S4Vectors' został zbudowany w wersji R 4.1.1

## Warning: pakiet 'GenomeInfoDb' został zbudowany w wersji R 4.1.1

library(sf)

## Warning: pakiet 'sf' został zbudowany w wersji R 4.1.2

library(DESeq2)

## Warning: pakiet 'MatrixGenerics' został zbudowany w wersji R 4.1.1

library(org.Hs.eg.db)
library(LSD)
library(geneplotter)
library(biomaRt)

## Warning: pakiet 'biomaRt' został zbudowany w wersji R 4.1.1

library(gplots)
library(VennDiagram)

## Warning: pakiet 'VennDiagram' został zbudowany w wersji R 4.1.2

library(vsn)
library(ggVennDiagram)
library(GenomicRanges)
library(gwasrapidd)
library(fgsea)
library(goseq)
library(ggplot2)
library(RColorBrewer)

## Warning: pakiet 'RColorBrewer' został zbudowany w wersji R 4.1.2

library(pheatmap)
library(TxDb.Mmusculus.UCSC.mm10.knownGene)

## Warning: pakiet 'GenomicFeatures' został zbudowany w wersji R 4.1.1

library(dplyr)
library(tidyverse)

## Warning: pakiet 'tidyverse' został zbudowany w wersji R 4.1.2

## Warning: pakiet 'tibble' został zbudowany w wersji R 4.1.2

## Warning: pakiet 'forcats' został zbudowany w wersji R 4.1.2

library(scales)
library(smoothmest)

## Warning: pakiet 'smoothmest' został zbudowany w wersji R 4.1.2

## Warning: pakiet 'MASS' został zbudowany w wersji R 4.1.2

library(BSgenome.Mmusculus.UCSC.mm10)
library(Matrix)
library(ggpubr)

## Warning: pakiet 'ggpubr' został zbudowany w wersji R 4.1.2

set.seed(22)

ms=12000*1024^2 
options(future.globals.maxSize=ms)

buildDfForGGBARPLOT = function( normCounts, metaTable, metaColumn, geneSymbol, mappings4genes ){
  # normCounts=ddsFdf; metaTable=as.data.frame(coldata);metaColumn="cell_type";geneSymbol="Aldh1a3";mappings4genes=genemap_unique
  theCol=which(colnames(metaTable)==metaColumn)
  data.frame( norm_counts = normCounts[rownames(normCounts) == unique(mappings4genes$gene_id[mappings4genes$geneName==geneSymbol])],
              annot = metaTable[ match(colnames(normCounts),metaTable$sample), theCol ] )}


source( paste0(scripts_directory,'functions.R' ) )

Building objects for future work

gtf = import.gff(paste0(data_directory,'Mus_musculus.GRCm38.101.gtf'))

txdb = TxDb.Mmusculus.UCSC.mm10.knownGene
columns <- c("tx_name", "gene_id","GENEID")
promGR = promoters(txdb, upstream=200, downstream=200,columns=columns)
TSSGR = promoters(txdb, upstream=0, downstream=0,columns=columns)
promGRExt = promoters(txdb, upstream=2000, downstream=2000,columns=columns)

## Warning in valid.GenomicRanges.seqinfo(x, suggest.trim = TRUE): GRanges object contains 5 out-of-bound ranges located on sequences
##   chr4_GL456350_random, chr4_JH584293_random, chr4_JH584295_random,
##   chr5_JH584296_random, and chrUn_GL456239. Note that ranges located on a
##   sequence whose length is unknown (NA) or on a circular sequence are not
##   considered out-of-bound (use seqlengths() and isCircular() to get the
##   lengths and circularity flags of the underlying sequences). You can use
##   trim() to trim these ranges. See ?`trim,GenomicRanges-method` for more
##   information.

promGRExt=trim(promGRExt)

transcrpit2ensg = data.frame(tid=gtf$transcript_id,gid=gtf$gene_id,stringsAsFactors = FALSE)
promGRExt$transcipt1 = unlist(strsplit(promGRExt$tx_name,'\\.'))[seq(1,2*length(promGRExt),by=2)]
promGRExt$ensgeneid = transcrpit2ensg$gid[match(promGRExt$transcipt1,promGRExt$transcipt1)]
seqlevelsStyle(promGR)='ucsc'
genemap_unique = data.frame(gene_id=gtf$gene_id,geneName=gtf$gene_name,stringsAsFactors = FALSE)
genemap_unique = genemap_unique[!duplicated(genemap_unique$gene_id),]

promoters_ap = data.frame( chr=as.character(chrom(gtf)),
                           start=as.numeric(start(gtf)),
                           end=as.numeric(end(gtf)),
                           strand=as.character(strand(gtf)),
                           transcript_id=as.character(gtf$transcript_id),
                           gene_id=as.character(gtf$gene_id),
                           gene_name = as.character(gtf$gene_name),
                           gene_biotype = as.character(gtf$gene_biotype),
                           type = gtf$type,
                           stringsAsFactors = FALSE )
length(unique(promoters_ap$gene_id))

## [1] 55487

promoters_ap = promoters_ap[promoters_ap$type == "transcript",]
promoters_sp = split(promoters_ap,promoters_ap$gene_id)

## for each gene find the 5-prime most TSS
promoters_tss = do.call('rbind', lapply( promoters_sp, function(x){
    tss = ifelse( as.character(unique(x$strand))=="+", 
                  x[which.min(x$start),'start'], 
                  x[which.max(x$end),'end'] ) 
    if(  as.character(unique(x$strand))=="+") tp = x[which.min(x$start),] else tp = x[which.max(x$end),]
    tp$tss = tss
  return(tp)
  }))

## 
promoters_tss_gr = GRanges(seqnames = promoters_tss$chr,
                           ranges = IRanges(as.numeric(promoters_tss$tss)-500 ,
                                            end=as.numeric(promoters_tss$tss) + 500,
                                            names=promoters_tss$transcript_id),
                           strand = promoters_tss$strand,
                           gene_id = promoters_tss$gene_id,
                           gene_name = promoters_tss$gene_name,
                           gene_biotype = promoters_tss$gene_biotype,
                           tss = promoters_tss$tss)
seqlevelsStyle(promoters_tss_gr) = 'ucsc'

ctcf_motif = readBed_filterChromsStraded( paste0(data_directory,'CTCF_HUMAN.H11MO.0.A.bed'),
                                          chroms=paste0('chr',c(1:19,'X','Y')), 5 )

Loops from Bonev et al

es_loops = importBEDPE_Filter_Loops( paste0(hic_directory, 'Bonev/hiccups_output_esc/' ),
                                     max_size=5000000, extend_anchors_by=0 )
ns_loops = importBEDPE_Filter_Loops( paste0(hic_directory, 'Bonev/hiccups_output_npc/' ),
                                     max_size=5000000, extend_anchors_by=0 )

es_loops_size=es_loops$y2-es_loops$x1
ns_loops_size=ns_loops$y2-ns_loops$x1

es_loops_left_anchor = getGR( es_loops$X.chr1, es_loops$x1, es_loops$x2,0 )
es_loops_right_anchor = getGR( es_loops$X.chr1, es_loops$y1, es_loops$y2,0 )

ns_loops_left_anchor = getGR( ns_loops$X.chr1, ns_loops$x1, ns_loops$x2,0 )
ns_loops_right_anchor = getGR( ns_loops$X.chr1, ns_loops$y1, ns_loops$y2,0 )

seqlevelsStyle(es_loops_left_anchor) = 'ucsc'
seqlevelsStyle(es_loops_right_anchor) = 'ucsc'
seqlevelsStyle(ns_loops_left_anchor) = 'ucsc'
seqlevelsStyle(ns_loops_right_anchor) = 'ucsc'

es_anchors_GR = c( es_loops_left_anchor, es_loops_right_anchor )
ns_anchors_GR = c( ns_loops_left_anchor, ns_loops_right_anchor )

es_ns_loops = compareLoops( es_loops, ns_loops, offset=5000, maxLoopSize=5000000 )
es_spe_loops = es_ns_loops$First_set_specific
ns_spe_loops = es_ns_loops$Second_set_specific
common_loops = es_ns_loops$First_set_common
es_spe_loops_size=es_spe_loops$y2-es_spe_loops$x1
ns_spe_loops_size=ns_spe_loops$y2-ns_spe_loops$x1


es_spe_loops_left_anchor = getGR( es_spe_loops$X.chr1, es_spe_loops$x1, es_spe_loops$x2,0 )
es_spe_loops_right_anchor = getGR( es_spe_loops$X.chr1, es_spe_loops$y1, es_spe_loops$y2,0 )

ns_spe_loops_left_anchor = getGR( ns_spe_loops$X.chr1, ns_spe_loops$x1, ns_spe_loops$x2,0 )
ns_spe_loops_right_anchor = getGR( ns_spe_loops$X.chr1, ns_spe_loops$y1, ns_spe_loops$y2,0 )

en_com_loops_left_anchor = getGR( common_loops$X.chr1, common_loops$x1, common_loops$x2,0 )
en_com_loops_right_anchor = getGR( common_loops$X.chr1, common_loops$y1, common_loops$y2,0 )

es_spe_loops_gr = getGR( es_spe_loops$X.chr1, es_spe_loops$x1, es_spe_loops$y2,0 )
ns_spe_loops_gr = getGR( ns_spe_loops$X.chr1, ns_spe_loops$x1, ns_spe_loops$y2,0 )
cm_loops_gr = getGR( common_loops$X.chr1, common_loops$x1, common_loops$y2,0 )

es_spe_loops_size=width(es_spe_loops_gr)
ns_spe_loops_size=width(ns_spe_loops_gr)

es_loops_gr = getGR( es_loops$X.chr1,es_loops$x1, es_loops$y2, 0 )
ns_loops_gr = getGR( ns_loops$X.chr1,ns_loops$x1, ns_loops$y2, 0 )

es_spe_loops_domain_gr = getGR( es_spe_loops$X.chr1, es_spe_loops$x1+25000, es_spe_loops$y2-25000,0 )
ns_spe_loops_domain_gr = getGR( ns_spe_loops$X.chr1, ns_spe_loops$x1+25000, ns_spe_loops$y2-25000,0 )
common_loops_domain_gr = getGR( common_loops$X.chr1, common_loops$x1+25000, common_loops$y2-25000,0 )

seqlevelsStyle(es_spe_loops_left_anchor) = 'ucsc'
seqlevelsStyle(es_spe_loops_right_anchor) = 'ucsc'
seqlevelsStyle(ns_spe_loops_left_anchor) = 'ucsc'
seqlevelsStyle(ns_spe_loops_right_anchor) = 'ucsc'
seqlevelsStyle(en_com_loops_left_anchor) = 'ucsc'
seqlevelsStyle(en_com_loops_right_anchor) = 'ucsc'
seqlevelsStyle(es_spe_loops_gr) = 'ucsc'
seqlevelsStyle(ns_spe_loops_gr) = 'ucsc'
seqlevelsStyle(es_loops_gr) = 'ucsc'
seqlevelsStyle(ns_loops_gr) = 'ucsc'
seqlevelsStyle(es_spe_loops_domain_gr) = 'ucsc'
seqlevelsStyle(ns_spe_loops_domain_gr) = 'ucsc'
seqlevelsStyle(common_loops_domain_gr) = 'ucsc'
seqlevelsStyle(cm_loops_gr) = 'ucsc'


es_spe_loops_filt = es_spe_loops[-queryHits(findOverlaps(es_spe_loops_gr,c(en_com_loops_left_anchor,en_com_loops_right_anchor))),]
ns_spe_loops_filt = ns_spe_loops[-queryHits(findOverlaps(ns_spe_loops_gr,c(en_com_loops_left_anchor,en_com_loops_right_anchor))),]
es_spe_loops_filt_gr = getGR(es_spe_loops_filt$X.chr1,es_spe_loops_filt$x1,es_spe_loops_filt$y2,0)
ns_spe_loops_filt_gr = getGR(ns_spe_loops_filt$X.chr1,ns_spe_loops_filt$x1,ns_spe_loops_filt$y2,0)

Analysis of transcriptome with and without CTCF (ES and NS cells)

Check differentiation of ES 46 and Nora ES cells

rna_46C = data.frame(es1=read.delim(paste0(rnaseq_directory,'gene_counts_FEB_2022.txt'),skip=1)[,13],
                     es2=read.delim(paste0(rnaseq_directory,'gene_counts_FEB_2022.txt'),skip=1)[,14],
                     ns1=read.delim(paste0(rnaseq_directory,'gene_counts_FEB_2022.txt'),skip=1)[,8],
                     ns2=read.delim(paste0(rnaseq_directory,'gene_counts_FEB_2022.txt'),skip=1)[,10],
                     row.names=read.delim(paste0(rnaseq_directory,'gene_counts_FEB_2022.txt'),skip=1)[,1])
coldata46c = data.frame(cell_type=c('ES','ES','NS','NS'),genotype='46C',row.names=colnames(rna_46C))

load( paste0(rnaseq_directory,"CTCF_degron_countdata.RData") )
count_df_filt = requested_countdata[,1:14]
coldata = requested_metadata[1:14,]
coldata$culture[9:14]='NS'
coldata$cell_type = paste(coldata$culture,coldata$condition,sep="_")
ddsF = DESeqDataSetFromMatrix(countData = count_df_filt,
                              colData = coldata,
                              design= ~cell_type )

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors

ddsF = DESeq(ddsF)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

ddsFdf = counts(ddsF,normalized=TRUE )

ES–>NS validation

Diff_Nora__DDS = DESeqDataSetFromMatrix(countData = count_df_filt[,c(3,4,7,8,10,11,14)],
                                        colData = coldata[c(3,4,7,8,10,11,14),],
                                        design= ~0+culture )

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors

es2i_ns_comp = results(DESeq(Diff_Nora__DDS), contrast=c("culture", "2i", "NS") )

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

es2i_ns_comp$padj_squished=es2i_ns_comp$padj
es2i_ns_comp$padj_squished[es2i_ns_comp$padj_squished<10^(-20)]=10^(-20)

par(pty="s")
plot(x=es2i_ns_comp$log2FoldChange,y=-log10(es2i_ns_comp$padj_squished),
     col=ifelse(es2i_ns_comp$padj<0.01,"green4","gray"),
     ylim=c(0,20),pch=19,cex=0.5)
axis(1,lwd=2)
axis(2,lwd=2)
box(col="black",lwd=2)

df = data.frame(LFC=es2i_ns_comp$log2FoldChange,log10pval=-log10(es2i_ns_comp$padj_squished))
write.table(df,file=paste0(source_data_directory,"volcano_differentiationNcells.txt"),row.names = FALSE, sep="\t",quote=FALSE)

Diff_Nora__DDS = DESeqDataSetFromMatrix(countData = count_df_filt[,c(3,4,7,8,10,11,14)],
                                        colData = coldata[c(3,4,7,8,10,11,14),],
                                        design= ~0+culture )

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors

es2i_ns_comp_counts = counts(DESeq(Diff_Nora__DDS),normalized=TRUE )

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

es2i_ns_comp = es2i_ns_comp[!is.na(es2i_ns_comp$padj),]
thisColV = rep("gray",nrow(es2i_ns_comp_counts))
thisColV[rownames(es2i_ns_comp_counts) %in% rownames(es2i_ns_comp[es2i_ns_comp$padj<0.01 & es2i_ns_comp$log2FoldChange>0,])] = "red"
thisColV[rownames(es2i_ns_comp_counts) %in% rownames(es2i_ns_comp[es2i_ns_comp$padj<0.01 & es2i_ns_comp$log2FoldChange<0,])] = "blue"

par(pty="s",bty="O")
plot(x=log2(0.5+rowMeans(es2i_ns_comp_counts[,1:4])),
     y=log2(0.5+rowMeans(es2i_ns_comp_counts[,5:7])),
     col=thisColV,xlab="ES", ylab="NS",
     ylim=c(-2,20),xlim=c(-2,20),pch=19,cex=0.25)
axis(1,lwd=2)
axis(2,lwd=2)
box(col="black",lwd=2)

df = data.frame(ES=log2(0.5+rowMeans(es2i_ns_comp_counts[,1:4])),NS=log2(0.5+rowMeans(es2i_ns_comp_counts[,5:7])))
write.table(df,file=paste0(source_data_directory,"scatterplot_Differentiation_Ncells.txt"),row.names = FALSE, sep="\t",quote=FALSE)

CTCF targets cell type specific genes

LFC_thr=0;pval_thr=0.1
ti_DDS = DESeqDataSetFromMatrix(countData = count_df_filt[,1:8],
                                colData = coldata[1:8,],
                                design= ~cell_type )

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors

ti_DDS=DESeq(ti_DDS)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

ti = results(ti_DDS, contrast=c("cell_type", "2i_WT", "2i_IAA") )


ns_DDS = DESeqDataSetFromMatrix(countData = count_df_filt[,9:14],
                                colData = coldata[9:14,],
                                design= ~0+cell_type )

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors

ns_DDS = DESeq(ns_DDS)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

ns = results(ns_DDS, contrast=c("cell_type", "NS_WT", "NS_IAA") )

es_ns_IAA_merged = merge( as.data.frame(ti), as.data.frame(ns), by='row.names')


ti = ti[!is.na(ti$padj),]
ns = ns[!is.na(ns$padj),]

TI_sig_genes = ( ti[ti$padj<pval_thr & abs(ti$log2FoldChange)>LFC_thr,] )
TI_sig_genes_up = rownames( ti[ti$padj<pval_thr & ti$log2FoldChange<(-1*LFC_thr),] )
TI_sig_genes_dn = rownames( ti[ti$padj<pval_thr & ti$log2FoldChange>LFC_thr,] )
TI_common = rownames( ti[ti$padj>0.1 & abs(ti$log2FoldChange)<log2(1.25),] )

NS_sig_genes = ( ns[ns$padj<pval_thr & abs(ns$log2FoldChange)>LFC_thr,] )
NS_sig_genes_up = rownames( ns[ns$padj<pval_thr & ns$log2FoldChange<(-1 *LFC_thr),] )
NS_sig_genes_dn = rownames( ns[ns$padj<pval_thr & ns$log2FoldChange>LFC_thr,] )
NS_common = rownames( ns[ns$padj>0.1 & abs(ns$log2FoldChange)<log2(1.25),] )

write.table( TI_sig_genes_up, file=paste0(outputs_directory,'TI_sig_genes_up.txt'),
             quote=FALSE, row.names=FALSE, col.names=FALSE, sep='\t' )
write.table( TI_sig_genes_dn, file=paste0(outputs_directory,'TI_sig_genes_dn.txt'),
             quote=FALSE, row.names=FALSE, col.names=FALSE, sep='\t' )
write.table( NS_sig_genes_up, file=paste0(outputs_directory,'NS_sig_genes_up.txt'),
             quote=FALSE, row.names=FALSE, col.names=FALSE, sep='\t' )
write.table( NS_sig_genes_dn, file=paste0(outputs_directory,'NS_sig_genes_dn.txt'),
             quote=FALSE, row.names=FALSE, col.names=FALSE, sep='\t' )

length(unique(c(TI_sig_genes_up,TI_sig_genes_dn,NS_sig_genes_up,NS_sig_genes_dn)))

## [1] 1250

Volcano plots for the extended figure

TI = data.frame(Untreated=rowMeans(ddsFdf[,c('ES_2i_1','ES_2i_2','ES_2i_3','ES_2i_4')]),
                IAA = rowMeans(ddsFdf[,c('ES_2i_IAA_1','ES_2i_IAA_2','ES_2i_IAA_3','ES_2i_IAA_4')]),
                stringsAsFactors=FALSE )

par(mar=c(5,5,3,1))
plot(ti$log2FoldChange,-log10(ti$pvalue), 
     col=ifelse(abs(ti$log2FoldChange)>0 & ti$padj<pval_thr,'red3','gray'), 
     pch=19, cex=0.25,
     xlab=expression('log'[2]*'(Untr/IAA)'),
     ylab=expression('-log'[10]*'(P-val)'),
     ylim=c(0,30),xlim=c(-3,3),
     main='ES 2i untr. vs. IAA')
axis(1,lwd=3)
axis(2,lwd=3)
box(lwd=3,col='black')

df = data.frame(LFC=ti$log2FoldChange,log10pval=-log10(ti$pvalue),col=ifelse(abs(ti$log2FoldChange)>0 & ti$padj<pval_thr,'red3','gray'))
write.table(df,file=paste0(source_data_directory,"volcano_ES_IAA.txt"),row.names = FALSE, sep="\t",quote=FALSE)

NS = data.frame(Untreated=rowMeans(ddsFdf[,c('NS_FBS_3','NS_FBS_5','NS_FBS_4')]),
                IAA = rowMeans(ddsFdf[,c('NS_FBS_IAA_3','NS_FBS_IAA_5','NS_FBS_IAA_4')]),
                stringsAsFactors=FALSE )

par(mar=c(5,5,3,1))
plot(ns$log2FoldChange,-log10(ns$pvalue), 
     col=ifelse(abs(ns$log2FoldChange)>0 & ns$padj<pval_thr,'red3','gray'), 
     pch=19, cex=0.25,
     ylim=c(0,30),xlim=c(-3,3),
     xlab=expression('log'[2]*'(Untr/IAA)'),
     ylab=expression('-log'[10]*'(P-val)'),
     main='NS')
axis(1,lwd=3)
axis(2,lwd=3)
box(lwd=3,col='black')

df = data.frame(LFC=ns$log2FoldChange,log10pval=-log10(ns$pvalue),col=ifelse(abs(ns$log2FoldChange)>0 & ns$padj<pval_thr,'red3','gray'))
write.table(df,file=paste0(source_data_directory,"volcano_NS_IAA.txt"),row.names = FALSE, sep="\t",quote=FALSE)

More up-regulated genes in the NS cells

pval_thr_plot=0.1
x = ti[ti$padj<pval_thr_plot,]
y = ns[ns$padj<pval_thr_plot,]
m=rbind(ES_2i=table( x$log2FoldChange>0 ),
        NS=table( y$log2FoldChange>0 ) )
m

##       FALSE TRUE
## ES_2i   357  418
## NS      358  198

fisher.test(m)

## 
##  Fisher's Exact Test for Count Data
## 
## data:  m
## p-value = 4.223e-11
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.3751452 0.5944608
## sample estimates:
## odds ratio 
##  0.4726291

par(mar=c(5,10,5,10),mfrow=c(1,1))
barplot(100*(m/rowSums(m)), 
        beside=T, ylim=c(0,80),
        col=c('red','blue'), ylab='%',
        names=c('Up','Down'),lwd=1.5)
fisher.test(matrix(c(49,51,65,35),2,2))

## 
##  Fisher's Exact Test for Count Data
## 
## data:  matrix(c(49, 51, 65, 35), 2, 2)
## p-value = 0.03188
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.2815690 0.9484736
## sample estimates:
## odds ratio 
##   0.519077

axis(1,lwd=3,at=c(2,5),c('Up','Down'))
axis(2,lwd=3)

Direction of change

es_ns_IAA_merged_filt = es_ns_IAA_merged[! is.na(es_ns_IAA_merged$padj.x),]
es_ns_IAA_merged_filt = es_ns_IAA_merged_filt[! is.na(es_ns_IAA_merged_filt$padj.y),]
es_ns_IAA_merged_filt = es_ns_IAA_merged_filt[rowSums(es_ns_IAA_merged_filt[,c('padj.x','padj.y')]<0.1)>0,]
es_ns_IAA_merged_filt$geneName = genemap_unique$geneName[match( es_ns_IAA_merged_filt$Row.names,genemap_unique$gene_id ) ]

es_ns_IAA_merged_filt$col = 'black'
es_ns_IAA_merged_filt$col[es_ns_IAA_merged_filt$padj.x<0.1 & es_ns_IAA_merged_filt$padj.y<0.1] = 'green4'
es_ns_IAA_merged_filt$col[es_ns_IAA_merged_filt$padj.x<0.1 & es_ns_IAA_merged_filt$padj.y>0.1] = 'red'
es_ns_IAA_merged_filt$col[es_ns_IAA_merged_filt$padj.x>0.1 & es_ns_IAA_merged_filt$padj.y<0.1] = 'blue'


es_ns_IAA_merged_filt$size = 0.5
es_ns_IAA_merged_filt$size[es_ns_IAA_merged_filt$padj.x<0.1 & es_ns_IAA_merged_filt$padj.y<0.1] = 0.5

par(mar=c(5,5,3,1))
plot(-1*es_ns_IAA_merged_filt$log2FoldChange.x,
     -1*es_ns_IAA_merged_filt$log2FoldChange.y, 
     col=es_ns_IAA_merged_filt$col, 
     pch=19, cex=es_ns_IAA_merged_filt$size,
     ylim=c(-4,4),xlim=c(-4,4),
     xlab=expression('ES [log'[2]*'(IAA/Untreated)]'),
     ylab=expression('NS [log'[2]*'(IAA/Untreated)]'),
     main='Without versus with CTCF')
abline( v=0,col='gray',lwd=3)
abline( h=0,col='gray',lwd=3)
axis(1,lwd=3)
axis(2,lwd=3)
box(col='black',lwd=3)

points( es_ns_IAA_merged_filt$log2FoldChange.x[es_ns_IAA_merged_filt$Row.names == 'ENSMUSG00000024268'],
        es_ns_IAA_merged_filt$log2FoldChange.y[es_ns_IAA_merged_filt$Row.names == 'ENSMUSG00000024268'], pch=19, cex=0.5, col='black')

cor.test(es_ns_IAA_merged_filt$log2FoldChange.x, es_ns_IAA_merged_filt$log2FoldChange.y)

## 
##  Pearson's product-moment correlation
## 
## data:  es_ns_IAA_merged_filt$log2FoldChange.x and es_ns_IAA_merged_filt$log2FoldChange.y
## t = 13.692, df = 979, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.3469996 0.4521298
## sample estimates:
##       cor 
## 0.4008836

es_ns_IAA_merged_filt_genes = es_ns_IAA_merged_filt[!is.na(es_ns_IAA_merged_filt$geneName),]

df = data.frame(ES=-1*es_ns_IAA_merged_filt$log2FoldChange.x,NS=-1*es_ns_IAA_merged_filt$log2FoldChange.y,col=es_ns_IAA_merged_filt$col)
write.table(df,file=paste0(source_data_directory,"scatterplot_IAAimpact_Ncells.txt"),row.names = FALSE, sep="\t",quote=FALSE)

Aldh1a3 expression plot

Plots for Figure 10

tp=buildDfForGGBARPLOT( ddsFdf,as.data.frame(coldata),
                     "cell_type",
                     "Sox9",
                     genemap_unique )
tp$annot = factor(tp$annot,levels = c("NS_WT","NS_IAA","2i_WT","2i_IAA"))
ggbarplot(tp, x = "annot", y = "norm_counts", 
          color="annot",
          add = c("mean_se", "jitter"), size=1.25,
          palette=c('blue','blue4','red','coral3'),
          title="Sox9") + theme( axis.line = element_line(colour = 'black', size = 0.75))

## Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

tp

##    norm_counts  annot
## 1     5.448699 2i_IAA
## 2     1.689030 2i_IAA
## 3     3.136609  2i_WT
## 4     4.230981  2i_WT
## 5     7.634598 2i_IAA
## 6     5.358430 2i_IAA
## 7     4.846634  2i_WT
## 8     5.989234  2i_WT
## 9   978.291257 NS_IAA
## 10 1530.401815  NS_WT
## 11 1868.185338  NS_WT
## 12  890.618161 NS_IAA
## 13  889.752300 NS_IAA
## 14 1371.565757  NS_WT

tp=buildDfForGGBARPLOT( ddsFdf,as.data.frame(coldata),
                     "cell_type",
                     "Ngfr",
                     genemap_unique )
tp$annot = factor(tp$annot,levels = c("NS_WT","NS_IAA","2i_WT","2i_IAA"))
ggbarplot(tp, x = "annot", y = "norm_counts", 
          color="annot",
          add = c("mean_se", "jitter"), size=1.25,
          palette=c('blue','blue4','red','coral3'),
          title="Ngfr") + theme( axis.line = element_line(colour = 'black', size = 0.75)) + ylim(c(0,125))

## Warning: Removed 8 rows containing non-finite outside the scale range
## (`stat_summary()`).

## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_bar()`).

## Warning: Removed 8 rows containing missing values or values outside the scale range
## (`geom_point()`).

tp

##    norm_counts  annot
## 1    530.34002 2i_IAA
## 2    913.76539 2i_IAA
## 3    564.58968  2i_WT
## 4    829.27229  2i_WT
## 5    329.67584 2i_IAA
## 6    306.62127 2i_IAA
## 7    428.44241  2i_WT
## 8    387.05422  2i_WT
## 9    101.76855 NS_IAA
## 10    47.31081  NS_WT
## 11    44.81172  NS_WT
## 12   110.09226 NS_IAA
## 13    94.72400 NS_IAA
## 14    55.29920  NS_WT

tp=buildDfForGGBARPLOT( ddsFdf,as.data.frame(coldata),
                     "cell_type",
                     "Camk2a",
                     genemap_unique )
tp$annot = factor(tp$annot,levels = c("NS_WT","NS_IAA","2i_WT","2i_IAA"))
ggbarplot(tp, x = "annot", y = "norm_counts", 
          color="annot",
          add = c("mean_se", "jitter"), size=1.25,
          palette=c('blue','blue4','red','coral3'),
          title="Camk2a") + theme( axis.line = element_line(colour = 'black', size = 0.75)) + ylim(c(0,600))

tp

##    norm_counts  annot
## 1     187.0720 2i_IAA
## 2     167.2140 2i_IAA
## 3     196.0381  2i_WT
## 4     205.2026  2i_WT
## 5     263.7407 2i_IAA
## 6     270.8984 2i_IAA
## 7     268.5035  2i_WT
## 8     243.3126  2i_WT
## 9     591.7349 NS_IAA
## 10    165.2450  NS_WT
## 11    285.8679  NS_WT
## 12    338.0397 NS_IAA
## 13    595.7813 NS_IAA
## 14    344.1648  NS_WT

Aldh1a3 CTCF-BS KOs

library(dplyr) 
library(ggplot2)
library(ggpubr)
theme_set(theme_bw())

mm = read.delim(paste0(data_directory,"Aldh1a3_ES-NS qPCR.txt"),header=TRUE)
mm = mm[-1,]
mm$Genotype=factor(mm$Genotype,levels=c("Wild_type","KO_1","KO_2","KO_3"))

mmES = mm[mm$Cell_type=="ES",]
m = ggplot(mmES, aes(x = Genotype, y = Expression))

m + geom_jitter(
  aes(shape = Genotype, color = Genotype), 
  position = position_jitter(0.2),
  size = 5 ) +
  stat_summary(
    aes(color = Genotype),
    fun.data="mean_sdl",  fun.args = list(mult=1), 
    geom = "pointrange",  size = 0.7) + scale_color_manual(values =  c("red","coral3", "coral3", "coral3")) + ylim(0,1) + theme( axis.line = element_line(colour = "black", 
                      linewidth = 1, linetype = "solid"))

mmNS = mm[mm$Cell_type=="NS",]
m = ggplot(mmNS, aes(x = Genotype, y = Expression))

m + geom_jitter(
  aes(shape = Genotype, color = Genotype), 
  position = position_jitter(0.1),
  size = 5 ) +
  stat_summary(
    aes(color = Genotype),
    fun.data="mean_sdl",  fun.args = list(mult=1), 
    geom = "pointrange",  size = 0.7 ) + scale_color_manual(values =  c("blue","steelblue3", "steelblue3", "steelblue3")) + ylim(0,0.5) + theme( axis.line = element_line(colour = "black", 
                      linewidth = 1, linetype = "solid"))

t.test( mmNS$Expression[mmNS$Genotype=="Wild_type"],mmNS$Expression[mmNS$Genotype=="KO_1"])

## 
##  Welch Two Sample t-test
## 
## data:  mmNS$Expression[mmNS$Genotype == "Wild_type"] and mmNS$Expression[mmNS$Genotype == "KO_1"]
## t = -4.7588, df = 3.6289, p-value = 0.0113
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.35052505 -0.08552556
## sample estimates:
##  mean of x  mean of y 
## 0.07367433 0.29169964

Check if the differentiation between ES2i and NS worked in the Nora and the 46C cells.

marker_genes=c("Prdm14","Nanog","Fgf4","Nodal","Lin28a","Zfp42","Dppa2","Esrrb","Tbx3","Klf4","Pou5f1","Phc1","Tet1","Dll1","Pax6","Notch1","Rxra","Hes6","Wnt5a","Fzd2","Fzd1","Nes","Bmi1","Hes5","Ascl1","Olig1","Olig2")

ens2gene = data.frame(gtf$gene_name,gtf$gene_id)
ens2gene = ens2gene[!duplicated(ens2gene$gtf.gene_id),]
marker_genes = data.frame(gene_names=marker_genes,
                          gene_id = ens2gene$gtf.gene_id[match(marker_genes,ens2gene$gtf.gene_name)])

marker_genes$lfc = -1*es2i_ns_comp$log2FoldChange[match(marker_genes$gene_id,rownames(es2i_ns_comp))]
marker_genes$padj = es2i_ns_comp$padj[match(marker_genes$gene_id,rownames(es2i_ns_comp))]
marker_genes = marker_genes[order(marker_genes$lfc,decreasing=FALSE),]

barplot(marker_genes$lfc,col=ifelse(marker_genes$lfc<0,"red","blue"),
        names=marker_genes$gene_names,las=2,ylim=c(-20,20),cex.axis = 1.5)

df = data.frame(LFC=marker_genes$lfc,col=ifelse(marker_genes$lfc<0,"red","blue"))
write.table(df,file=paste0(source_data_directory,"barplot_ES__NS_differentiation.txt"))

Expression of up and down regulated genes

rnL = read.delim(paste0(rnaseq_directory,'gene_counts_FEB_2022.txt'),skip=1)[,c(1,6)]
requested_countdata_tpm = data.frame( ES = rowSums( requested_countdata[,c(3,4,7,8)] ),
                                      NS = rowSums( requested_countdata[,c(10,11,14)] ),
                                      Length = rnL$Length[match(rownames(requested_countdata),rnL$Geneid)] )
requested_countdata_tpm = GetTPM(requested_countdata_tpm,1:2,rownames(requested_countdata_tpm))

par(mfrow=c(1,2),bty='n')
boxplot( requested_countdata_tpm[rownames(requested_countdata_tpm) %in% TI_sig_genes_up, "ES"],
         requested_countdata_tpm[rownames(requested_countdata_tpm) %in% TI_sig_genes_dn, "ES"],
         outline=FALSE,col='white', ylab="Expression (TPM)",
         border=c('red','red4'),lwd=3,ylim=c(0,100))
axis(1,lwd=3,at=c(1,2),c("Up","Down"))
axis(2,lwd=3)
boxplot( requested_countdata_tpm[rownames(requested_countdata_tpm) %in% NS_sig_genes_up, "NS"],
         requested_countdata_tpm[rownames(requested_countdata_tpm) %in% NS_sig_genes_dn, "NS"],
         outline=FALSE,col='white', border=c('blue','steelblue3'),
         ylab="Expression (TPM)",
         lwd=3,ylim=c(0,100))
axis(1,lwd=3,at=c(1,2),c("Up","Down"))
axis(2,lwd=3)

t.test(requested_countdata_tpm[rownames(requested_countdata_tpm) %in% TI_sig_genes_up, "ES"],
       requested_countdata_tpm[rownames(requested_countdata_tpm) %in% TI_sig_genes_dn, "ES"])

## 
##  Welch Two Sample t-test
## 
## data:  requested_countdata_tpm[rownames(requested_countdata_tpm) %in% TI_sig_genes_up, "ES"] and requested_countdata_tpm[rownames(requested_countdata_tpm) %in% TI_sig_genes_dn, "ES"]
## t = -5.2704, df = 441.46, p-value = 2.134e-07
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -102.59209  -46.86077
## sample estimates:
## mean of x mean of y 
##  12.29998  87.02641

t.test(requested_countdata_tpm[rownames(requested_countdata_tpm) %in% NS_sig_genes_up, "NS"],
         requested_countdata_tpm[rownames(requested_countdata_tpm) %in% NS_sig_genes_dn, "NS"])

## 
##  Welch Two Sample t-test
## 
## data:  requested_countdata_tpm[rownames(requested_countdata_tpm) %in% NS_sig_genes_up, "NS"] and requested_countdata_tpm[rownames(requested_countdata_tpm) %in% NS_sig_genes_dn, "NS"]
## t = -2.5177, df = 210.1, p-value = 0.01256
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -83.12087 -10.11752
## sample estimates:
## mean of x mean of y 
##  21.75334  68.37254

theNtimes = max(c(length(TI_sig_genes_up),length(TI_sig_genes_dn),length(NS_sig_genes_up),length(NS_sig_genes_dn)))
df = data.frame( ES_up_by_IAA = rep(NA,theNtimes), ES_dn_by_IAA = rep(NA,theNtimes),
                 NS_up_by_IAA = rep(NA,theNtimes), NS_dn_by_IAA = rep(NA,theNtimes) )
df$ES_up_by_IAA[1:length(TI_sig_genes_up)] = requested_countdata_tpm[rownames(requested_countdata_tpm) %in% TI_sig_genes_up, "ES"]
df$ES_dn_by_IAA[1:length(TI_sig_genes_dn)] = requested_countdata_tpm[rownames(requested_countdata_tpm) %in% TI_sig_genes_dn, "ES"]
df$NS_up_by_IAA[1:length(NS_sig_genes_up)] = requested_countdata_tpm[rownames(requested_countdata_tpm) %in% NS_sig_genes_up, "ES"]
df$NS_dn_by_IAA[1:length(NS_sig_genes_dn)] = requested_countdata_tpm[rownames(requested_countdata_tpm) %in% NS_sig_genes_dn, "ES"]
write.table(df,file=paste0(source_data_directory,"boxplots_TPM.txt"),row.names = FALSE, sep="\t",quote=FALSE)

Chromatin signature with and without CTCF - ATAC-seq and ChIP-seq data

ATAC peaks +/- CTCF

ti_atac_peaks_Nora_unt = readBed_filterChroms(paste0(chipseq_directory,'ES_2i_PLNOV_merged_filtered.bam_peaks.narrowPeak'),
                                     chroms=paste0('chr',c(1:19,'X','Y')), 7)
ti_atac_peaks_Nora_iaa = readBed_filterChroms(paste0(chipseq_directory,'ES_2i_IAA_PLNOV_merged_filtered.bam_peaks.narrowPeak'),
                                     chroms=paste0('chr',c(1:19,'X','Y')), 7)

ns_atac_peaks_Nora_unt = readBed_filterChroms(paste0(chipseq_directory,'NPC_PLDEC_merged_filtered.bam_peaks.narrowPeak'),
                                     chroms=paste0('chr',c(1:19,'X','Y')), 7)
ns_atac_peaks_Nora_iaa = readBed_filterChroms(paste0(chipseq_directory,'NPC_IAA_PLDEC_merged_filtered.bam_peaks.narrowPeak'),
                                     chroms=paste0('chr',c(1:19,'X','Y')), 7)

ti_atac_peaks_Nora_unt_summit = readNarrowPeak2getSummit(paste0(chipseq_directory,'ES_2i_PLNOV_merged_filtered.bam_peaks.narrowPeak'),
                                     chroms=paste0('chr',c(1:19,'X','Y')), 5)
ti_atac_peaks_Nora_iaa_summit = readNarrowPeak2getSummit(paste0(chipseq_directory,'ES_2i_IAA_PLNOV_merged_filtered.bam_peaks.narrowPeak'),
                                     chroms=paste0('chr',c(1:19,'X','Y')), 5)

ns_atac_peaks_Nora_unt_summit = readNarrowPeak2getSummit(paste0(chipseq_directory,'NPC_PLDEC_merged_filtered.bam_peaks.narrowPeak'),
                                     chroms=paste0('chr',c(1:19,'X','Y')), 5)
ns_atac_peaks_Nora_iaa_summit = readNarrowPeak2getSummit(paste0(chipseq_directory,'NPC_IAA_PLDEC_merged_filtered.bam_peaks.narrowPeak'),
                                     chroms=paste0('chr',c(1:19,'X','Y')), 5)

overlaps=findOverlaps(ti_atac_peaks_Nora_unt,ti_atac_peaks_Nora_iaa)
draw.pairwise.venn(area1=queryLength(overlaps),
                   area2=subjectLength(overlaps),
                   cross.area=length(overlaps),
                   category=c('',''), col=c('red','coral4'),
                   fill=rep('white',2),cat.cex=1.2,lwd=8)

## (polygon[GRID.polygon.374], polygon[GRID.polygon.375], polygon[GRID.polygon.376], polygon[GRID.polygon.377], text[GRID.text.378], text[GRID.text.379], text[GRID.text.380], text[GRID.text.381], text[GRID.text.382])

overlaps=findOverlaps(ns_atac_peaks_Nora_unt,ns_atac_peaks_Nora_iaa)
draw.pairwise.venn(area1=queryLength(overlaps),
                   area2=subjectLength(overlaps),
                   cross.area=length(overlaps),
                   category=c('',''), col=c('blue','blue4'),
                   fill=rep('white',2),cat.cex=1.2,lwd=8)

## (polygon[GRID.polygon.383], polygon[GRID.polygon.384], polygon[GRID.polygon.385], polygon[GRID.polygon.386], text[GRID.text.387], text[GRID.text.388], text[GRID.text.389], text[GRID.text.390], text[GRID.text.391])

H3K27ac peaks

ti_k27ac_peaks_Nora_unt = readBed_filterChroms(paste0(chipseq_directory,'ChIP_Seq_H3K27ac_05-22_MusMus_ESC_AID_KI_CTCF-AID-GFP_OsTIR_TIGRE_Nora_2i_Rep_1_peaks.narrowPeak'),
                                     chroms=paste0('chr',c(1:19,'X','Y')), 7)
ti_k27ac_peaks_Nora_iaa = readBed_filterChroms(paste0(chipseq_directory,'ChIP_Seq_H3K27ac_05-22_MusMus_ESC_AID_KI_CTCF-AID-GFP_OsTIR_TIGRE_Nora_2i_IAA_Rep_1_peaks.narrowPeak'),
                                     chroms=paste0('chr',c(1:19,'X','Y')), 7)

ns_k27ac_peaks_Nora_unt = readBed_filterChroms(paste0(chipseq_directory,'ChIP_Seq_H3K27ac_05-22_MusMus_es-NPC_AID_KI_CTCF-AID-GFP_OsTIR_TIGRE_Nora_Rep_1_peaks.narrowPeak'),
                                     chroms=paste0('chr',c(1:19,'X','Y')), 7)
ns_k27ac_peaks_Nora_iaa = readBed_filterChroms(paste0(chipseq_directory,'ChIP_Seq_H3K27ac_05-22_MusMus_es-NPC_AID_KI_CTCF-AID-GFP_OsTIR_TIGRE_Nora_IAA_Rep_1_peaks.narrowPeak'),
                                     chroms=paste0('chr',c(1:19,'X','Y')), 7)

overlaps=findOverlaps(ti_k27ac_peaks_Nora_unt,ti_k27ac_peaks_Nora_iaa)
draw.pairwise.venn(area1=queryLength(overlaps),
                   area2=subjectLength(overlaps),
                   cross.area=length(overlaps),
                   category=c('',''), col=c('red','coral4'),
                   fill=rep('white',2),cat.cex=1.2,lwd=8)

## (polygon[GRID.polygon.392], polygon[GRID.polygon.393], polygon[GRID.polygon.394], polygon[GRID.polygon.395], text[GRID.text.396], text[GRID.text.397], text[GRID.text.398], text[GRID.text.399], text[GRID.text.400])

overlaps=findOverlaps(ns_k27ac_peaks_Nora_unt,ns_k27ac_peaks_Nora_iaa)
draw.pairwise.venn(area1=queryLength(overlaps),
                   area2=subjectLength(overlaps),
                   cross.area=length(overlaps),
                   category=c('',''), col=c('blue','blue4'),
                   fill=rep('white',2),cat.cex=1.2,lwd=8)

## (polygon[GRID.polygon.401], polygon[GRID.polygon.402], polygon[GRID.polygon.403], polygon[GRID.polygon.404], text[GRID.text.405], text[GRID.text.406], text[GRID.text.407], text[GRID.text.408], text[GRID.text.409])

Peak number changes and deregulation upon CTCF loss

es_up_proms = GenomicRanges::resize(promoters_tss_gr[promoters_tss_gr$gene_id %in% TI_sig_genes_up],500000,fix='center')
es_dn_proms = GenomicRanges::resize(promoters_tss_gr[promoters_tss_gr$gene_id %in% TI_sig_genes_dn],500000,fix='center')

ns_up_proms = GenomicRanges::resize(promoters_tss_gr[promoters_tss_gr$gene_id %in% NS_sig_genes_up],500000,fix='center')
ns_dn_proms = GenomicRanges::resize(promoters_tss_gr[promoters_tss_gr$gene_id %in% NS_sig_genes_dn],500000,fix='center')

par(mfrow=c(1,4),mar=c(6,5,1,1),pty="m")
es_up_promc = data.frame(unt=countOverlaps(es_up_proms,ti_k27ac_peaks_Nora_unt),
                         iaa=countOverlaps(es_up_proms,ti_k27ac_peaks_Nora_iaa))
es_dn_promc = data.frame(unt=countOverlaps(es_dn_proms,ti_k27ac_peaks_Nora_unt),
                         iaa=countOverlaps(es_dn_proms,ti_k27ac_peaks_Nora_iaa))
boxplot( es_dn_promc$iaa-es_dn_promc$unt, 
         es_up_promc$iaa-es_up_promc$unt, 
         border=c("red4","coral4"),
         outline=FALSE, col="white",names=c("Down","Up"),
         ylim=c(-10,10),lwd=2, ylab="IAA-veh")
write.table(es_up_promc,file=paste0(source_data_directory,"es_up_promc.txt"),quote=FALSE,row.names = FALSE, sep="\t")
write.table(es_dn_promc,file=paste0(source_data_directory,"es_dn_promc.txt"),quote=FALSE,row.names = FALSE, sep="\t")


es_up_promc = data.frame(unt=countOverlaps(es_up_proms,ti_atac_peaks_Nora_unt),
                         iaa=countOverlaps(es_up_proms,ti_atac_peaks_Nora_iaa))
es_dn_promc = data.frame(unt=countOverlaps(es_dn_proms,ti_atac_peaks_Nora_unt),
                         iaa=countOverlaps(es_dn_proms,ti_atac_peaks_Nora_iaa))

boxplot( es_dn_promc$iaa-es_dn_promc$unt, 
         es_up_promc$iaa-es_up_promc$unt, 
         border=c("red4","coral4"),
         outline=FALSE, col="white", names=c("Down","Up"),
         ylim=c(-10,10), lwd=2, ylab="IAA-veh" )
write.table(es_up_promc,file=paste0(source_data_directory,"es_up_promc.txt"),quote=FALSE,row.names = FALSE, sep="\t")
write.table(es_dn_promc,file=paste0(source_data_directory,"es_dn_promc.txt"),quote=FALSE,row.names = FALSE, sep="\t")

ns_dn_promc = data.frame(unt=countOverlaps(ns_dn_proms,ns_k27ac_peaks_Nora_unt),
                         iaa=countOverlaps(ns_dn_proms,ns_k27ac_peaks_Nora_iaa))
ns_up_promc = data.frame(unt=countOverlaps(ns_up_proms,ns_k27ac_peaks_Nora_unt),
                         iaa=countOverlaps(ns_up_proms,ns_k27ac_peaks_Nora_iaa))
boxplot( ns_dn_promc$unt-ns_dn_promc$iaa, 
         ns_up_promc$unt-ns_up_promc$iaa, 
         border=c("blue4","steelblue3"),
         outline=FALSE, col="white",names=c("Down","Up"),
         ylim=c(-10,10), lwd=2, ylab="veh-IAA" )
write.table(ns_dn_promc,file=paste0(source_data_directory,"ns_dn_promc.txt"),quote=FALSE,row.names = FALSE, sep="\t")
write.table(ns_up_promc,file=paste0(source_data_directory,"ns_up_promc.txt"),quote=FALSE,row.names = FALSE, sep="\t")

ns_dn_promc = data.frame(unt=countOverlaps(ns_dn_proms,ns_atac_peaks_Nora_unt),
                         iaa=countOverlaps(ns_dn_proms,ns_atac_peaks_Nora_iaa))
ns_up_promc = data.frame(unt=countOverlaps(ns_up_proms,ns_atac_peaks_Nora_unt),
                         iaa=countOverlaps(ns_up_proms,ns_atac_peaks_Nora_iaa))
boxplot( ns_dn_promc$unt-ns_dn_promc$iaa, 
         ns_up_promc$unt-ns_up_promc$iaa, 
         border=c("blue4","steelblue3"),
         outline=FALSE, col="white",names=c("Down","Up"),
         ylim=c(-10,10), lwd=2, ylab="veh-IAA" )

write.table(ns_dn_promc,file=paste0(source_data_directory,"ns_dn_promc.txt"),quote=FALSE,row.names = FALSE, sep="\t")
write.table(ns_up_promc,file=paste0(source_data_directory,"ns_up_promc.txt"),quote=FALSE,row.names = FALSE, sep="\t")

CTCF controls - we define peaks in untreated cells and look at CTCF signal in untreated and treated cells

tictmm10Nora_peak = readNarrowPeak2getSummit( paste0(chipseq_directory, 'ChIP_Seq_CTCF_07-22_MusMus_ESC_AID_KI_CTCF-AID-GFP_OsTIR_TIGRE_2i_Rep_1_peaks.narrowPeak'),
                                     chroms=paste0('chr',c(1:19,'X','Y')), 5)

nsctmm10Nora_peak = readNarrowPeak2getSummit( paste0(chipseq_directory, 'ChIP_Seq_CTCF_05-22_MusMus_es-NPC_AID_KI_CTCF-AID-GFP_OsTIR_TIGRE_Nora_Rep_1_peaks.narrowPeak'),
                                     chroms=paste0('chr',c(1:19,'X','Y')), 5)

ti_nora_ctcf_untreated_bp = import.bw(paste0(chipseq_directory,'ChIP_Seq_CTCF_07-22_MusMus_ESC_AID_KI_CTCF-AID-GFP_OsTIR_TIGRE_2i_Rep_1_RPGC.bw'))
ti_nora_ctcf_iaa_bp = import.bw(paste0(chipseq_directory,'ChIP_Seq_CTCF_07-22_MusMus_ESC_AID_KI_CTCF-AID-GFP_OsTIR_TIGRE_2i-IAA_Rep_1_RPGC.bw'))
ns_nora_ctcf_untreated_bp = import.bw(paste0(chipseq_directory,'ChIP_Seq_CTCF_05-22_MusMus_es-NPC_AID_KI_CTCF-AID-GFP_OsTIR_TIGRE_Nora_Rep_1_RPGC.bw'))
ns_nora_ctcf_iaa_bp = import.bw(paste0(chipseq_directory,'ChIP_Seq_CTCF_05-22_MusMus_es-NPC_AID_KI_CTCF-AID-GFP_OsTIR_TIGRE_Nora_IAA_Rep_1_RPGC.bw'))

seqlevelsStyle(ti_nora_ctcf_untreated_bp)='ucsc'
seqlevelsStyle(ns_nora_ctcf_untreated_bp)='ucsc'
seqlevelsStyle(ti_nora_ctcf_iaa_bp)='ucsc'
seqlevelsStyle(ns_nora_ctcf_iaa_bp)='ucsc'

Ctcf_nora_summit_ranges_ES = do.call('rbind', lapply( as.list(seq(1:length(tictmm10Nora_peak))), function(i){
  g = tictmm10Nora_peak[i]
  tss = start(g)
  chr = as.character( chrom(g) )
  peak=i
      return( data.frame( chr = rep( (chr), 200),
                  starts = seq( tss-1000, tss+990, length.out=200),
                  ends =  seq( tss-1000, tss+990, length.out=200)+9,
                  peak = rep( peak, 200) ) )
} ) )
Ctcf_nora_summit_ranges_NS = do.call('rbind', lapply( as.list(seq(1:length(nsctmm10Nora_peak))), function(i){
  g = nsctmm10Nora_peak[i]
  tss = start(g)
  chr = as.character( chrom(g) )
  peak=i
      return( data.frame( chr = rep( (chr), 200),
                  starts = seq( tss-1000, tss+990, length.out=200),
                  ends =  seq( tss-1000, tss+990, length.out=200)+9,
                  peak = rep( peak, 200) ) )
} ) )

Ctcf_nora_summit_ranges_ES = GRanges(seqnames = Rle(Ctcf_nora_summit_ranges_ES$chr),
                               ranges = IRanges(as.numeric(Ctcf_nora_summit_ranges_ES$starts),
                                                end = as.numeric(Ctcf_nora_summit_ranges_ES$ends),
                                                names = seq(1, nrow(Ctcf_nora_summit_ranges_ES))),
                               strand = Rle(rep("*", nrow(Ctcf_nora_summit_ranges_ES))),
                               peak = Ctcf_nora_summit_ranges_ES$peak )

Ctcf_nora_summit_ranges_NS = GRanges(seqnames = Rle(Ctcf_nora_summit_ranges_NS$chr),
                               ranges = IRanges(as.numeric(Ctcf_nora_summit_ranges_NS$starts),
                                                end = as.numeric(Ctcf_nora_summit_ranges_NS$ends),
                                                names = seq(1, nrow(Ctcf_nora_summit_ranges_NS))),
                               strand = Rle(rep("*", nrow(Ctcf_nora_summit_ranges_NS))),
                               peak = Ctcf_nora_summit_ranges_NS$peak )

save(Ctcf_nora_summit_ranges_ES,file=paste0(objects_directory,'Ctcf_nora_summit_ranges_ES.RData'))
save(Ctcf_nora_summit_ranges_NS,file=paste0(objects_directory,'Ctcf_nora_summit_ranges_NS.RData'))

Ctcf_nora_summit_ranges_ES_untreated = getSignalInBins( Ctcf_nora_summit_ranges_ES, ti_nora_ctcf_untreated_bp, 1 )
Ctcf_nora_summit_ranges_ES_iaa = getSignalInBins( Ctcf_nora_summit_ranges_ES, ti_nora_ctcf_iaa_bp, 1 )

Ctcf_nora_summit_ranges_NS_untreated = getSignalInBins( Ctcf_nora_summit_ranges_NS, ns_nora_ctcf_untreated_bp, 1 )
Ctcf_nora_summit_ranges_NS_iaa = getSignalInBins( Ctcf_nora_summit_ranges_NS, ns_nora_ctcf_iaa_bp, 1 )

save( Ctcf_nora_summit_ranges_ES_untreated, 
      Ctcf_nora_summit_ranges_ES_iaa, 
      Ctcf_nora_summit_ranges_NS_untreated, 
      Ctcf_nora_summit_ranges_NS_iaa,
      file=paste0(objects_directory,'Ctcf_signal_Nora_CTCF_peaks.RData'))

Display the effect of IAA on CTCF level.

load(paste0(objects_directory,'Ctcf_signal_Nora_CTCF_peaks.RData'))
x=log2((0.1+Ctcf_nora_summit_ranges_ES_iaa)/(0.1+Ctcf_nora_summit_ranges_ES_untreated))
y=log2((0.1+Ctcf_nora_summit_ranges_NS_iaa)/(0.1+Ctcf_nora_summit_ranges_NS_untreated))

par(mfrow=c(1,1),mar=c(5,5,1,1))
plot(density(x[,100]),col='red',lwd=3,ty='l',
     xlim=c(-10,10),main='',axes=F,ylim=c(0,0.5))
axis(1,lwd=3)
axis(2,lwd=3)
lines(density(y[,100]),col='blue',lwd=3 )
abline(v=0,lwd=3,col='black',lty=2)

t.test(x[,100])

## 
##  One Sample t-test
## 
## data:  x[, 100]
## t = -185.78, df = 13849, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -1.985657 -1.944194
## sample estimates:
## mean of x 
## -1.964926

t.test(y[,100])

## 
##  One Sample t-test
## 
## data:  y[, 100]
## t = -213.78, df = 18626, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -2.215984 -2.175719
## sample estimates:
## mean of x 
## -2.195851

df = data.frame(ES=density(x[,100])$y,NS=density(y[,100])$y)
write.table(df,file=paste0(source_data_directory,"ChIP-seq_FC_IAA.txt"),row.names = FALSE, sep="\t",quote=FALSE)

Peaks in ES and NS cells

ti_atac_peaks = readBed_filterChroms(paste0(chipseq_directory,'ATAC_Seq_03-22_MusMus_ESC_MOD_SOX1-GFP_46C_2i_Rep_1_peaks.narrowPeak'),
                                     chroms=paste0('chr',c(1:19,'X','Y')), 7)
ns_atac_peaks = readBed_filterChroms(paste0(chipseq_directory,'ATAC_Seq_05-22_MusMus_es-NPC_MOD_SOX1-GFP_46C_Rep_1_peaks.narrowPeak'),
                                     chroms=paste0('chr',c(1:19,'X','Y')), 7)

ti_atac_summits = readNarrowPeak2getSummit(paste0(chipseq_directory,'ATAC_Seq_03-22_MusMus_ESC_MOD_SOX1-GFP_46C_2i_Rep_1_peaks.narrowPeak'),
                                     chroms=paste0('chr',c(1:19,'X','Y')), 5)
ns_atac_summits = readNarrowPeak2getSummit(paste0(chipseq_directory,'ATAC_Seq_05-22_MusMus_es-NPC_MOD_SOX1-GFP_46C_Rep_1_peaks.narrowPeak'),
                                     chroms=paste0('chr',c(1:19,'X','Y')), 5)


ctcf_motif$presence = 1

K27peaks_ES = readBed_filterChroms(paste0(chipseq_directory,'ChIP_Seq_12-21_MusMus_es-ESC_H3K27ac_MOD_SOX1-GFP_2i_Rep_1.narrowPeak'),
                                     chroms=paste0('chr',c(1:19,'X','Y')), 7)
K27peaks_NS = readBed_filterChroms( paste0(chipseq_directory, 'ChIP_Seq_12-21_MusMus_es-NPC_H3K27ac_MOD_SOX1-GFP_Rep_1.narrowPeak'),
                                     chroms=paste0('chr',c(1:19,'X','Y')), 7)

K27peaks_summit_ES  = readNarrowPeak2getSummit(paste0(chipseq_directory, 'ChIP_Seq_12-21_MusMus_es-ESC_H3K27ac_MOD_SOX1-GFP_2i_Rep_1.narrowPeak'),
                                         chroms=paste0('chr',c(1:19,'X','Y')), 5)

K27peaks_summit_NS = readNarrowPeak2getSummit(paste0(chipseq_directory, 'ChIP_Seq_12-21_MusMus_es-NPC_H3K27ac_MOD_SOX1-GFP_Rep_1.narrowPeak'),
                                         chroms=paste0('chr',c(1:19,'X','Y')), 5)


ti_enhancers = K27peaks_ES[-queryHits(findOverlaps(K27peaks_ES,promGRExt))]
ns_enhancers = K27peaks_NS[-queryHits(findOverlaps(K27peaks_NS,promGRExt))]

ti_enhancers_atac = ti_atac_peaks[queryHits(findOverlaps(ti_atac_peaks,ti_enhancers))]
ns_enhancers_atac = ns_atac_peaks[queryHits(findOverlaps(ns_atac_peaks,ns_enhancers))]

CTCF, gene expression, and loop gain

Induced genes are closer to each other than reduced genes

par(mfrow=c(1,2),mar=c(5,5,1,1),cex.axis=1)
boxplot(elementMetadata(distanceToNearest(promoters_tss_gr[which(promoters_tss_gr$gene_id %in% NS_sig_genes_up)]))$distance/1000,
        elementMetadata(distanceToNearest(promoters_tss_gr[which(promoters_tss_gr$gene_id %in% NS_sig_genes_dn)]))$distance/1000,
        outline=FALSE,border=c('green4','orange3'),col='white',names=c("Up","Down"),ylab="distance (kbp)")
ks.test(elementMetadata(distanceToNearest(promoters_tss_gr[which(promoters_tss_gr$gene_id %in% NS_sig_genes_up)]))$distance,
        elementMetadata(distanceToNearest(promoters_tss_gr[which(promoters_tss_gr$gene_id %in% NS_sig_genes_dn)]))$distance)

## Warning in
## ks.test(elementMetadata(distanceToNearest(promoters_tss_gr[which(promoters_tss_gr$gene_id
## %in% : p-value will be approximate in the presence of ties

## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  elementMetadata(distanceToNearest(promoters_tss_gr[which(promoters_tss_gr$gene_id %in% NS_sig_genes_up)]))$distance and elementMetadata(distanceToNearest(promoters_tss_gr[which(promoters_tss_gr$gene_id %in% NS_sig_genes_dn)]))$distance
## D = 0.24865, p-value = 3.431e-07
## alternative hypothesis: two-sided

## distance to nearest enhancer
boxplot(elementMetadata(distanceToNearest(promoters_tss_gr[which(promoters_tss_gr$gene_id %in% NS_sig_genes_up)],ns_enhancers_atac))$distance/1000,
        elementMetadata(distanceToNearest(promoters_tss_gr[which(promoters_tss_gr$gene_id %in% NS_sig_genes_dn)],ns_enhancers_atac))$distance/1000,
        outline=FALSE,border=c('green4','orange3'),col='white',names=c("Up","Down"))

Boxplot data

theNtimes = max( c(length(NS_sig_genes_up),length(NS_sig_genes_dn)))
df = data.frame( Prom_prom_up = rep(NA,theNtimes), 
                 Prom_prom_dn = rep(NA,theNtimes),
                 Prom_enh_up = rep(NA,theNtimes), 
                 Prom_enh_dn = rep(NA,theNtimes) )
df$Prom_prom_up[1:length(elementMetadata(distanceToNearest(promoters_tss_gr[which(promoters_tss_gr$gene_id %in% NS_sig_genes_up)]))$distance/1000)] = elementMetadata(distanceToNearest(promoters_tss_gr[which(promoters_tss_gr$gene_id %in% NS_sig_genes_up)]))$distance/1000
df$Prom_prom_dn[1:length(elementMetadata(distanceToNearest(promoters_tss_gr[which(promoters_tss_gr$gene_id %in% NS_sig_genes_dn)]))$distance/1000)] = elementMetadata(distanceToNearest(promoters_tss_gr[which(promoters_tss_gr$gene_id %in% NS_sig_genes_dn)]))$distance/1000

df$Prom_enh_up[1:length(elementMetadata(distanceToNearest(promoters_tss_gr[which(promoters_tss_gr$gene_id %in% NS_sig_genes_up)],ns_enhancers_atac))$distance/1000)] = elementMetadata(distanceToNearest(promoters_tss_gr[which(promoters_tss_gr$gene_id %in% NS_sig_genes_up)],ns_enhancers_atac))$distance/1000
df$Prom_enh_dn[1:length(elementMetadata(distanceToNearest(promoters_tss_gr[which(promoters_tss_gr$gene_id %in% NS_sig_genes_dn)],ns_enhancers_atac))$distance/1000)] = elementMetadata(distanceToNearest(promoters_tss_gr[which(promoters_tss_gr$gene_id %in% NS_sig_genes_dn)],ns_enhancers_atac))$distance/1000

write.table(df,file=paste0(source_data_directory,"boxplots_Prom_enh_Prom_Prom_distance.txt"),row.names = FALSE, sep="\t",quote=FALSE)

Insulators in the NS cells

Do we see that architectural CTCF peaks separate the induced genes from enhancers more often than reduced genes?

Find the smallest loop domain for each up and down gene, identify regions that flank the insulators (500 upstream), find the region with the highest number of enhancers and compare for induced and reduced genes. We will also check the number of CTCF peaks.

nsctmm10peak = readBed_filterChroms( paste0(chipseq_directory, 'NS_CTCF_ChIP_NIAMS_merged_filtered_peaks.narrowPeak'),
                                     chroms=paste0('chr',c(1:19,'X','Y')), 7)

nsctmm10 = readNarrowPeak2getSummit( paste0(chipseq_directory,'NS_CTCF_ChIP_NIAMS_merged_filtered_peaks.narrowPeak'),
                                     chroms=paste0('chr',c(1:19,'X','Y')), 5)
nsctmm10$score = nsctmm10peak$score
nsctmm10peak = appendMotifInformation( nsctmm10peak, ctcf_motif,summitFile=nsctmm10 )


ctcf_peaks_ns_mm10_strand = nsctmm10peak[which(nsctmm10peak$motif_strand != "*")]
strand(ctcf_peaks_ns_mm10_strand) = ctcf_peaks_ns_mm10_strand$motif_strand


sum( elementMetadata(distanceToNearest(promoters_tss_gr[which(promoters_tss_gr$gene_id %in% NS_sig_genes_up)],ctcf_peaks_ns_mm10_strand,ignore.strand=TRUE))$distance == 0)/length(NS_sig_genes_up)

## [1] 0.1731844

sum( elementMetadata(distanceToNearest(promoters_tss_gr[which(promoters_tss_gr$gene_id %in% NS_sig_genes_dn)],ctcf_peaks_ns_mm10_strand,ignore.strand=TRUE))$distance == 0)/length(NS_sig_genes_dn)

## [1] 0.4545455

m=matrix(c(sum( elementMetadata(distanceToNearest(promoters_tss_gr[which(promoters_tss_gr$gene_id %in% NS_sig_genes_up)],ctcf_peaks_ns_mm10_strand,ignore.strand=TRUE))$distance == 0),
                     length(NS_sig_genes_up),
                     sum( elementMetadata(distanceToNearest(promoters_tss_gr[which(promoters_tss_gr$gene_id %in% NS_sig_genes_dn)],ctcf_peaks_ns_mm10_strand,ignore.strand=TRUE))$distance == 0),
                     length(NS_sig_genes_dn)),2,2)
fisher.test(m)

## 
##  Fisher's Exact Test for Count Data
## 
## data:  m
## p-value = 2.369e-07
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  0.2591761 0.5589507
## sample estimates:
## odds ratio 
##  0.3815533

par(mfrow=c(1,1))
barplot(100*m[1,]/m[2,],col=c('green4','orange3'),
        ylab="%",ylim=c(0,50),names=c("up","down"))
axis(2,lwd=2)

##      [,1] [,2]
## [1,]   62   90
## [2,]  358  198

How many enhancers flank the loop domains that contain different classes of genes?

dis2anchor_induced = c()
induced_intervals_left = data.frame()
induced_intervals_right = data.frame()
induced_smallest_loop = c()
induced_largest_loop = c()


for( gene in NS_sig_genes_up ){
  thisGene=promoters_tss_gr[which(promoters_tss_gr$gene_id %in% gene)]
  tp = findOverlaps(ns_loops_gr,thisGene)
  if(length(tp)>0){
  loops = ns_loops_gr[queryHits(tp)]
  loops2 = loops[which.max(width(loops))]
  loops = loops[which.min(width(loops))]
  
  induced_smallest_loop = rbind(induced_smallest_loop,
                                data.frame(chrom(loops),start(loops),end(loops)))
  induced_largest_loop = rbind(induced_largest_loop, data.frame(chrom(loops2),start(loops2),end(loops2)))
  
  dis2anchor_induced = c(dis2anchor_induced,elementMetadata(distanceToNearest(thisGene,c(ns_loops_left_anchor[loops$loop],ns_loops_right_anchor[loops$loop])))$distance)
  induced_intervals_left = rbind(induced_intervals_left,data.frame(chrom(loops),start(loops)-500000,start(loops)))
  induced_intervals_right = rbind(induced_intervals_right,data.frame(chrom(loops),end(loops),end(loops)+500000))}
}
induced_smallest_loop_gr = getGR(induced_smallest_loop[,1],induced_smallest_loop[,2],induced_smallest_loop[,3],0)

induced_intervals_left_gr = getGR(induced_intervals_left[,1],induced_intervals_left[,2],induced_intervals_left[,3],0)
induced_intervals_right_gr = getGR(induced_intervals_right[,1],induced_intervals_right[,2],induced_intervals_right[,3],0)
induced_smallest_loop_gr = getGR(induced_smallest_loop[,1],induced_smallest_loop[,2],induced_smallest_loop[,3],0)
induced_largest_loop_gr = getGR(induced_largest_loop[,1],induced_largest_loop[,2],induced_largest_loop[,3],0)

## reduced genes
dis2anchor_reduced = c()
reduced_intervals_left = data.frame()
reduced_intervals_right = data.frame()
reduced_smallest_loop = c()
reduced_largest_loop = c()

for( gene in NS_sig_genes_dn ){
  thisGene=promoters_tss_gr[which(promoters_tss_gr$gene_id %in% gene)]
  tp = findOverlaps(ns_loops_gr,thisGene)
  if(length(tp)>0){
  loops = ns_loops_gr[queryHits(tp)]
  loops2 = loops[which.max(width(loops))]
  loops = loops[which.min(width(loops))]
  reduced_smallest_loop = rbind(reduced_smallest_loop,
                                data.frame(chrom(loops),start(loops),end(loops)))
  reduced_largest_loop = rbind(reduced_largest_loop,
                                data.frame(chrom(loops2),start(loops2),end(loops2)))
  dis2anchor_reduced = c(dis2anchor_reduced,elementMetadata(distanceToNearest(thisGene,c(ns_loops_left_anchor[loops$loop],ns_loops_right_anchor[loops$loop])))$distance)
  reduced_intervals_left = rbind(reduced_intervals_left,data.frame(chrom(loops),start(loops)-500000,start(loops)))
  reduced_intervals_right = rbind(reduced_intervals_right,data.frame(chrom(loops),end(loops),end(loops)+500000))}
}


reduced_intervals_left_gr = getGR(reduced_intervals_left[,1],reduced_intervals_left[,2],reduced_intervals_left[,3],0)
reduced_intervals_right_gr = getGR(reduced_intervals_right[,1],reduced_intervals_right[,2],reduced_intervals_right[,3],0)
reduced_smallest_loop_gr = getGR(reduced_smallest_loop[,1],reduced_smallest_loop[,2],reduced_smallest_loop[,3],0)
reduced_largest_loop_gr = getGR(reduced_largest_loop[,1],reduced_largest_loop[,2],reduced_largest_loop[,3],0)

## ------------
dis2anchor_random = c()
random_intervals_left = data.frame()
random_intervals_right = data.frame()
random_smallest_loop = c()
random_largest_loop = c()

for( gene in promoters_tss_gr[sample(seq(1,length(promoters_tss_gr)),2000,replace = FALSE)]$gene_id ){
  thisGene=promoters_tss_gr[which(promoters_tss_gr$gene_id %in% gene)]
  tp = findOverlaps(ns_loops_gr,thisGene)
  if(length(tp)>0){
  loops = ns_loops_gr[queryHits(tp)]
  loops2 = loops[which.max(width(loops))]
  loops = loops[which.min(width(loops))]
  random_smallest_loop = rbind(random_smallest_loop,
                               data.frame(chrom(loops),start(loops),end(loops)))
  random_largest_loop = rbind(random_largest_loop,
                               data.frame(chrom(loops2),start(loops2),end(loops2)))
  dis2anchor_random = c(dis2anchor_random,elementMetadata(distanceToNearest(thisGene,c(ns_loops_left_anchor[loops$loop],ns_loops_right_anchor[loops$loop])))$distance)
  random_intervals_left = rbind(random_intervals_left,data.frame(chrom(loops),start(loops)-500000,start(loops)))
  random_intervals_right = rbind(random_intervals_right,data.frame(chrom(loops),end(loops),end(loops)+500000))}
}


random_intervals_left_gr = getGR(random_intervals_left[,1],random_intervals_left[,2],random_intervals_left[,3],0)
random_intervals_right_gr = getGR(random_intervals_right[,1],random_intervals_right[,2],random_intervals_right[,3],0)
random_smallest_loop_gr = getGR(random_smallest_loop[,1],random_smallest_loop[,2],random_smallest_loop[,3],0)
random_largest_loop_gr = getGR(random_largest_loop[,1],random_largest_loop[,2],random_largest_loop[,3],0)


## ----------
induced_left_right = data.frame(left=countOverlaps(induced_intervals_left_gr,
                                                   ns_enhancers_atac),

                                right = countOverlaps(induced_intervals_right_gr,
                                                      ns_enhancers_atac))

reduced_left_right = data.frame(left=countOverlaps(reduced_intervals_left_gr,
                                                   ns_enhancers_atac),

                                right = countOverlaps(reduced_intervals_right_gr,
                                                      ns_enhancers_atac))

random_left_right = data.frame(left=countOverlaps(random_intervals_left_gr,
                                                   ns_enhancers_atac),

                                right = countOverlaps(random_intervals_right_gr,
                                                      ns_enhancers_atac))

Check the landscape in the loops

induced_largest_loop_gr$loop = paste(chrom(induced_largest_loop_gr),start(induced_largest_loop_gr),sep="_")
reduced_largest_loop_gr$loop = paste(chrom(reduced_largest_loop_gr),start(reduced_largest_loop_gr),sep="_")
induced_largest_loop_gr_filt = induced_largest_loop_gr[which(! duplicated(induced_largest_loop_gr$loop))]
reduced_largest_loop_gr_filt = reduced_largest_loop_gr[which(! duplicated(reduced_largest_loop_gr$loop))]

boxplot(rowSums(induced_left_right[which(! duplicated(induced_largest_loop_gr$loop)),]),
        rowSums(reduced_left_right[which(! duplicated(reduced_largest_loop_gr$loop)),]),
        rowSums(random_left_right),col="white",border=c("green4","orange3","gray60"),
        names=c("up","down","random"),ylab="number of enhancers in flanks")
axis(1,lwd=2, at=c(1,2,3),c("up","down","random"))
axis(2,lwd=2)
box(col='black',lwd=2)

t.test(rowSums(induced_left_right[which(! duplicated(induced_largest_loop_gr$loop)),]),
       rowSums(reduced_left_right[which(! duplicated(reduced_largest_loop_gr$loop)),]))

## 
##  Welch Two Sample t-test
## 
## data:  rowSums(induced_left_right[which(!duplicated(induced_largest_loop_gr$loop)), ]) and rowSums(reduced_left_right[which(!duplicated(reduced_largest_loop_gr$loop)), ])
## t = 4.6081, df = 313.74, p-value = 5.919e-06
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  2.053708 5.114292
## sample estimates:
## mean of x mean of y 
##    12.000     8.416

t.test(rowSums(random_left_right[which(! duplicated(reduced_largest_loop_gr$loop)),]),rowSums(reduced_left_right))

## 
##  Welch Two Sample t-test
## 
## data:  rowSums(random_left_right[which(!duplicated(reduced_largest_loop_gr$loop)), ]) and rowSums(reduced_left_right)
## t = -2.8088, df = 243.77, p-value = 0.005376
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -3.6673615 -0.6439448
## sample estimates:
## mean of x mean of y 
##  7.064000  9.219653

theNtimes = max( length(rowSums(induced_left_right[which(! duplicated(induced_largest_loop_gr$loop)),])),
                 length(rowSums(reduced_left_right[which(! duplicated(reduced_largest_loop_gr$loop)),])),
                 length(rowSums(random_left_right)) )
df = data.frame( induced = rep(NA,theNtimes),
                 reduced = rep(NA,theNtimes),
                 random = rep(NA,theNtimes) )

df$induced[1:length(rowSums(induced_left_right[which(! duplicated(induced_largest_loop_gr$loop)),]))] = rowSums(induced_left_right[which(! duplicated(induced_largest_loop_gr$loop)),])
df$reduced[1:length(rowSums(reduced_left_right[which(! duplicated(reduced_largest_loop_gr$loop)),]))] = rowSums(reduced_left_right[which(! duplicated(reduced_largest_loop_gr$loop)),])
df$random[1:length(rowSums(random_left_right))] = rowSums(random_left_right)

write.table(df,file=paste0(source_data_directory,"boxplots_flanks.txt"),row.names = FALSE, sep="\t",quote=FALSE)

PLA

rnase = read.delim(paste0(data_directory,"PLA_microscopy/Fig-RNaseA-PLA-Ddx5-CTCF_Fus-CTCF-WT-ES_NS.txt"),
                   header=TRUE)
rnase$Cond=factor(paste(rnase$Cell.type, rnase$Treatment),
                  levels=c("ESC Untreated","ESC RNaseA", "NSC Untreated", "NSC RNaseA"))

par(bty="n",mfrow=c(2,1),mar=c(3,3,2,1))
boxplot( split(rnase[rnase$PLA=="Ddx5-Ctcf",2],
               rnase$Cond[rnase$PLA=="Ddx5-Ctcf"]),
         col="white",border=c("red","red","blue","blue"),
         lwd=2,ylim=c(0,20),axes=FALSE,main="Ddx5-Ctcf")
axis(1,c(1,2,3,4),
     c("","", "", ""),
     lwd=2)
axis(2,lwd=2)

boxplot( split(rnase[rnase$PLA=="Fus-Ctcf",2],
               rnase$Cond[rnase$PLA=="Fus-Ctcf" ]),
         col="white",border=c("red","red","blue","blue"),
         lwd=2,ylim=c(0,15),main="Fus-Ctcf",axes=FALSE)
axis(1,c(1,2,3,4),
     c("ESC Unt.","ESC RNaseA", "NSC Unt.", "NSC RNaseA"),
     lwd=2)
axis(2,lwd=2)

PLA : Ctcf-Nono

Ctcf_nono = read.delim(paste0(data_directory,"PLA_microscopy/Fig-PLA-Nono-Ctcf-WT-ES_NS.txt"),
                   header=TRUE)
x = split(Ctcf_nono[,2],Ctcf_nono$Cell.type)

par(bty="n",mfrow=c(1,1),mar=c(3,3,2,1),pty="m")
boxplot( x,
         col="white",border=c("red","blue"),
         lwd=2,ylim=c(0,40),axes=FALSE,main="Ctcf-Nono")
axis(1,c(1,2), c("ES","NS"), lwd=2)
axis(2,lwd=2)

t.test(x[[1]],x[[2]])

## 
##  Welch Two Sample t-test
## 
## data:  x[[1]] and x[[2]]
## t = -3.0945, df = 33.582, p-value = 0.003958
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -13.552177  -2.804965
## sample estimates:
## mean of x mean of y 
##  13.00000  21.17857

PLA : FUS-Ddx5

ddx_fus = read.delim(paste0(data_directory,"PLA_microscopy/Fig-PLA-Ddx5-Fus-WT-ES_NS.txt"),
                   header=TRUE)
par(bty="n",mfrow=c(1,1),mar=c(3,3,2,1))
boxplot( split(ddx_fus[,2],
               ddx_fus$Cell.type),
         col="white",border=c("red","blue"),
         lwd=2,ylim=c(0,40),axes=FALSE,main="Ddx5-Fus")
axis(1,c(1,2),
     c("ES","NS"),
     lwd=2)
axis(2,lwd=2)

PANTR1 PLA Ddx5-CTCF

pantr1 = read.delim(paste0(data_directory,"PLA_microscopy/PLA_Ddx5-CTCF_end_Dec_2023.txt"),
                   header=TRUE)
pantr1$Genotype = factor(pantr1$Genotype,levels=c("Wildtype-A2","Wildtype-A3",
                                                  "Pantr1-KO-PB6","Pantr1-KO-PE3"))


pantr1s = split(pantr1[,2],pantr1$Genotype)
par(bty="n",mfrow=c(1,1),mar=c(3,3,2,1))
boxplot( pantr1s[c(1,2,4,3)],
         col="white",border=c("blue","blue","gray60","gray60"),
         lwd=2,ylim=c(0,100),axes=FALSE,main="CTCF-Ddx5")
axis(1,c(1,2,3,4),
     c("Wt1","Wt2",
       "PE3","PB6"),
     lwd=2)
axis(2,lwd=2)

t.test(c(pantr1s[[1]],pantr1s[[2]]), c(pantr1s[[3]],pantr1s[[4]]))

## 
##  Welch Two Sample t-test
## 
## data:  c(pantr1s[[1]], pantr1s[[2]]) and c(pantr1s[[3]], pantr1s[[4]])
## t = 18.201, df = 129.76, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  38.40084 47.76715
## sample estimates:
## mean of x mean of y 
##  66.28947  23.20548

PANTR1 PLA FUS-CTCF

pantr1_fus = read.delim(paste0(data_directory,"PLA_microscopy/PLA-Fus-CTCF-Pantr-KO-new.txt"),
                   header=TRUE)
pantr1_fus$Genotype = factor(pantr1_fus$Genotype,
                             levels=c("Wildtype-A2","Wildtype-A3",
                                      "Pantr1-KO-PB6","Pantr1-KO-PE3"))
pantr1_fus = pantr1_fus[pantr1_fus$Genotype %in% c("Wildtype-A2","Wildtype-A3","Pantr1-KO-PB6","Pantr1-KO-PE3"),]
pantr1_fus$Genotype = droplevels(pantr1_fus$Genotype)

pantr1s = split(pantr1_fus[,2],pantr1_fus$Genotype)
par(bty="n",mfrow=c(1,1),mar=c(3,3,2,1))
boxplot( pantr1s,
         col="white",border=c("blue","blue","gray60","gray60"),
         lwd=2,ylim=c(0,100),axes=FALSE,main="CTCF-FUS")
axis(1,c(1,2,3,4),
     c("Wt1","Wt2",
       "PE3","PB6"),
     lwd=2)
axis(2,lwd=2)

t.test(c(unlist(pantr1s[c(1,2)])),c(unlist(pantr1s[c(3,4)])))

## 
##  Welch Two Sample t-test
## 
## data:  c(unlist(pantr1s[c(1, 2)])) and c(unlist(pantr1s[c(3, 4)]))
## t = 6.2177, df = 120.78, p-value = 7.483e-09
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  10.27250 19.87033
## sample estimates:
## mean of x mean of y 
##  58.63768  43.56627

qRT-PCR

pantr1_qPCR = c(0.34,0,0.01,0,0.01,0)
names(pantr1_qPCR) = c("Wt +","Wt -","PB6+","PB6-","PE3+","PE3-")
barplot(pantr1_qPCR,col=c("blue","blue",rep("gray",4)),ylim=c(0,0.4) )
axis(2,lwd=2)

CTCF clustering - AiryScan on pre-extracted nuclei

ctcf_clusters = read.delim(paste0(data_directory,"PLA_microscopy/Figure_CTCF-Clusters-WT_ES_NS.txt"),header=TRUE)

par(bty="n",mfrow=c(1,1),mar=c(3,3,2,1))
x=split(ctcf_clusters[,1],ctcf_clusters$Cell.type)
boxplot( x, col="white",border=c("red","blue"),
         lwd=2,ylim=c(0,200),axes=FALSE,main="CTCF cluster size")
axis(1,c(1,2),
     c('ES','NS'),
     lwd=2)
axis(2,lwd=2)

t.test(x[[1]],x[[2]])

## 
##  Welch Two Sample t-test
## 
## data:  x[[1]] and x[[2]]
## t = -5.9133, df = 64.346, p-value = 1.403e-07
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -67.54590 -33.43454
## sample estimates:
## mean of x mean of y 
##  19.67186  70.16207

CTCF peaks

Overwhelming majority of CTCF peaks intersects a motif.

ctcf_motif = readBed_filterChromsStraded( paste0(data_directory,'CTCF_HUMAN.H11MO.0.A.bed'),
                                          chroms=paste0('chr',c(1:19,'X','Y')), 5 )
ctcf_wt_NS_peaks = import.bed( paste0(chipseq_directory,'ChIP_SEQ_CTCF_Merged_rpgc_narrow_summits.bed'))
seqlevelsStyle(ctcf_wt_NS_peaks)="ucsc"
ctcf_wt_NS_peaks = ctcf_wt_NS_peaks[which(chrom(ctcf_wt_NS_peaks) %in% paste0("chr",c(1:19,"x")))]

ctcf_wt_NS_peaks = GenomicRanges::resize(ctcf_wt_NS_peaks,width=200,fix="center")
ctcf_wt_NS_peaks = appendMotifInformation( ctcf_wt_NS_peaks, ctcf_motif,summitFile=ctcf_wt_NS_peaks )
table(ctcf_wt_NS_peaks$motif_strand)

## 
##     *     +     - 
##  5038 18472 18722

ctcf_wt_NS_peaks_summit = GenomicRanges::resize(ctcf_wt_NS_peaks,fix="center",1)
names(ctcf_wt_NS_peaks_summit) = seq(1:length(ctcf_wt_NS_peaks_summit))
ctcf_wt_NS_peaks_resized = GenomicRanges::resize( ctcf_wt_NS_peaks_summit, fix="center",width=500 )

write.table( ctcf_wt_NS_peaks,
             file=paste0(data_directory,"ctcf_wt_NS_peaks.txt"),
             sep='\t',quote=FALSE,row.names=FALSE,col.names=TRUE )

CTCF AUC

a3_1_ctcf_rpgc = import.bw(paste0(chipseq_directory,"ChIP_Seq_CTCF_06-23_MusMus_es-NPC_MOD_CTCF-Cterm_HALO_A3_Control_Rep_1_RPGC.bw"))
a3_2_ctcf_rpgc = import.bw(paste0(chipseq_directory,"ChIP_Seq_CTCF_08-23_MusMus_es-NPC_MOD_CTCF-Cterm_HALO_A3_Control_Rep_2_RPGC.bw"))
cb1_ctcf_rpgc = import.bw(paste0(chipseq_directory,"ChIP_Seq_CTCF_03-23_MusMus_es-NPC_DDX5_KO_CTCF-Cterm_HALO_CB1_Rep_1_RPGC.bw"))
ce10_ctcf_rpgc = import.bw(paste0(chipseq_directory,"ChIP_Seq_CTCF_03-23_MusMus_es-NPC_DDX5_KO_CTCF-Cterm_HALO_CE10_Rep_1_RPGC.bw"))
seqlevelsStyle(a3_1_ctcf_rpgc) = "ucsc"
seqlevelsStyle(a3_2_ctcf_rpgc) = "ucsc"
seqlevelsStyle(cb1_ctcf_rpgc) = "ucsc"
seqlevelsStyle(ce10_ctcf_rpgc) = "ucsc"

a3_1_ctcf_raw = import.bw(paste0(chipseq_directory,"ChIP_Seq_CTCF_06-23_MusMus_es-NPC_MOD_CTCF-Cterm_HALO_A3_Control_Rep_1_filtered.bw"))
a3_2_ctcf_raw = import.bw(paste0(chipseq_directory,"ChIP_Seq_CTCF_08-23_MusMus_es-NPC_MOD_CTCF-Cterm_HALO_A3_Control_Rep_2_filtered.bw"))
cb1_ctcf_raw = import.bw(paste0(chipseq_directory,"ChIP_Seq_CTCF_03-23_MusMus_es-NPC_DDX5_KO_CTCF-Cterm_HALO_CB1_Rep_1_filtered.bw"))
ce10_ctcf_raw = import.bw(paste0(chipseq_directory,"ChIP_Seq_CTCF_03-23_MusMus_es-NPC_DDX5_KO_CTCF-Cterm_HALO_CE10_Rep_1_filtered.bw"))
seqlevelsStyle(a3_1_ctcf_raw) = "ucsc"
seqlevelsStyle(a3_2_ctcf_raw) = "ucsc"
seqlevelsStyle(cb1_ctcf_raw) = "ucsc"
seqlevelsStyle(ce10_ctcf_raw) = "ucsc"

ctcf_NS_ranges = do.call('rbind', lapply( as.list(seq(1:length(ctcf_wt_NS_peaks_summit))), 
                                          function(i){
  g = ctcf_wt_NS_peaks_summit[i]
  tss = start(g)
  chr = as.character( chrom(g) )
  peak=names(ctcf_wt_NS_peaks_summit)[i]
      return( data.frame( chr = rep( (chr), 200),
                          starts = seq( tss-1000, tss+990, length.out=200),
                          ends =  seq( tss-1000, tss+990, length.out=200)+9,
                          peak = rep( peak, 200) ) ) } ) )


ctcf_NS_ranges = GRanges(seqnames = Rle(ctcf_NS_ranges$chr),
                         ranges = IRanges(as.numeric(ctcf_NS_ranges$starts),
                                          end = as.numeric(ctcf_NS_ranges$ends),
                                          names = seq(1, nrow(ctcf_NS_ranges))),
                         strand = Rle(rep("*", nrow(ctcf_NS_ranges))),
                         peak = ctcf_NS_ranges$peak )

save(ctcf_NS_ranges,file=paste0(objects_directory,'ctcf_NS_ranges.RData'))

a3_1_ctcf_rpgc_ctcf_AP = getSignalInBins(ctcf_NS_ranges,a3_1_ctcf_rpgc,1)
a3_2_ctcf_rpgc_ctcf_AP = getSignalInBins(ctcf_NS_ranges,a3_2_ctcf_rpgc,1)
cb1_ctcf_rpgc_ctcf_AP = getSignalInBins(ctcf_NS_ranges,cb1_ctcf_rpgc,1)
ce10_ctcf_rpgc_ctcf_AP = getSignalInBins(ctcf_NS_ranges,ce10_ctcf_rpgc,1)
rownames(a3_1_ctcf_rpgc_ctcf_AP) = names(ctcf_wt_NS_peaks_summit)
rownames(a3_2_ctcf_rpgc_ctcf_AP) = names(ctcf_wt_NS_peaks_summit)
rownames(cb1_ctcf_rpgc_ctcf_AP) = names(ctcf_wt_NS_peaks_summit)
rownames(ce10_ctcf_rpgc_ctcf_AP) = names(ctcf_wt_NS_peaks_summit)

save(a3_1_ctcf_rpgc_ctcf_AP,
     a3_2_ctcf_rpgc_ctcf_AP,
     cb1_ctcf_rpgc_ctcf_AP,
     ce10_ctcf_rpgc_ctcf_AP,
     file=paste0(objects_directory,"CTCF_Ddx5ko_NS_A3_RPGC.RData"))

a3_1_ctcf_raw_ctcf_AP = getSignalInBins(ctcf_NS_ranges,a3_1_ctcf_raw,1)
a3_2_ctcf_raw_ctcf_AP = getSignalInBins(ctcf_NS_ranges,a3_2_ctcf_raw,1)
cb1_ctcf_raw_ctcf_AP = getSignalInBins(ctcf_NS_ranges,cb1_ctcf_raw,1)
ce10_ctcf_raw_ctcf_AP = getSignalInBins(ctcf_NS_ranges,ce10_ctcf_raw,1)
rownames(a3_1_ctcf_raw_ctcf_AP) = names(ctcf_wt_NS_peaks_summit)
rownames(a3_2_ctcf_raw_ctcf_AP) = names(ctcf_wt_NS_peaks_summit)
rownames(cb1_ctcf_raw_ctcf_AP) = names(ctcf_wt_NS_peaks_summit)
rownames(ce10_ctcf_raw_ctcf_AP) = names(ctcf_wt_NS_peaks_summit)

save(a3_1_ctcf_raw_ctcf_AP,
     a3_2_ctcf_raw_ctcf_AP,
     cb1_ctcf_raw_ctcf_AP,
     ce10_ctcf_raw_ctcf_AP,
     file=paste0(objects_directory,"CTCF_Ddx5ko_NS_A3_RAW.RData"))

load(paste0(objects_directory,"CTCF_Ddx5ko_NS_A3_RAW.RData"))
load(paste0(objects_directory,'ctcf_NS_ranges.RData'))
a=90; b=110

elbat_ctcf = data.frame( a3_1 = ( rowSums(a3_1_ctcf_raw_ctcf_AP[,a:b])),
                         a3_2 = ( rowSums(a3_2_ctcf_raw_ctcf_AP[,a:b])),
                         cb1 = ( rowSums(cb1_ctcf_raw_ctcf_AP[,a:b])),
                         ce10 = ( rowSums(ce10_ctcf_raw_ctcf_AP[,a:b])) )
metadata_CTCF_peaks = data.frame( genotype = c("wt","wt","Ddx5","Ddx5"),
                                  clone = c("a3_1","a3_2","cb1","ce10"),
                                  row.names = colnames(elbat_ctcf))
Ddx5_ctcf_Deseq2_RPGC = DESeqDataSetFromMatrix(
 countData = elbat_ctcf,
 colData = metadata_CTCF_peaks,
 design = ~ genotype )

## converting counts to integer mode

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors

Ddx5_ctcf_Deseq2_RPGC = DESeq(Ddx5_ctcf_Deseq2_RPGC, fitType = 'local')

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

res_Ddx5_ctcf_Deseq2_RPGC = results(Ddx5_ctcf_Deseq2_RPGC, 
                              contrast = c("genotype", "Ddx5", "wt"))
summary( res_Ddx5_ctcf_Deseq2_RPGC )

## 
## out of 42232 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0 (up)       : 343, 0.81%
## LFC < 0 (down)     : 398, 0.94%
## outliers [1]       : 0, 0%
## low counts [2]     : 0, 0%
## (mean count < 7)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

ctcf_normalized = counts( Ddx5_ctcf_Deseq2_RPGC,normalized=TRUE )

heatscatter( x=res_Ddx5_ctcf_Deseq2_RPGC$baseMean,
             y=res_Ddx5_ctcf_Deseq2_RPGC$log2FoldChange, 
             pch=19, cex=0.1, add.contour = T, 
             xlim=c(0,1200),ylim=c(-2,2),
             xlab=expression(Log[2] ~ (Signal)),
             ylab=expression(Log[2] ~ (LFC)), 
             main="WT mutant")
abline( h=0 )
axis(1,lwd=2)
axis(2,lwd=2)
box(col="black",lwd=2)

df = data.frame( M=res_Ddx5_ctcf_Deseq2_RPGC$baseMean,A=res_Ddx5_ctcf_Deseq2_RPGC$log2FoldChange)
write.table(df,file=paste0(source_data_directory,"MA_ddx5KO.txt"),row.names = FALSE, sep="\t",quote=FALSE)

res_Ddx5_ctcf_Deseq2_RPGC = res_Ddx5_ctcf_Deseq2_RPGC[!is.na(res_Ddx5_ctcf_Deseq2_RPGC$pvalue),]
res_Ddx5_ctcf_Deseq2_RPGC = res_Ddx5_ctcf_Deseq2_RPGC[!is.na(res_Ddx5_ctcf_Deseq2_RPGC$padj),]

res_Ddx5_ctcf_Deseq2_RPGC = res_Ddx5_ctcf_Deseq2_RPGC[!is.na(res_Ddx5_ctcf_Deseq2_RPGC$padj),]
res_Ddx5_ctcf_Deseq2_RPGC$col = "gray90"

res_Ddx5_ctcf_Deseq2_RPGC$col[res_Ddx5_ctcf_Deseq2_RPGC$padj<0.25 & res_Ddx5_ctcf_Deseq2_RPGC$log2FoldChange>0]="#049494"
res_Ddx5_ctcf_Deseq2_RPGC$col[res_Ddx5_ctcf_Deseq2_RPGC$padj<0.25 & res_Ddx5_ctcf_Deseq2_RPGC$log2FoldChange<0 ]="#305494"
res_Ddx5_ctcf_Deseq2_RPGC$pvalue[res_Ddx5_ctcf_Deseq2_RPGC$pvalue<0.0000000001] = 0.0000000001


par(mfrow=c(1,1),mar=c(5,5,5,5))
plot(res_Ddx5_ctcf_Deseq2_RPGC$log2FoldChange,
     -log10(res_Ddx5_ctcf_Deseq2_RPGC$pvalue),
     col=res_Ddx5_ctcf_Deseq2_RPGC$col,pch=19, cex=0.25,
     ylab=expression(-Log[10] ~ (P-val.)),
     xlab=expression(Log[2] ~ (Ddx5/WT)), 
     axes=FALSE, xlim=c(-3,3),cex.lab=1.5 )
axis(1,lwd=2,cex.axis=1.5)
axis(2,lwd=2,cex.axis=1.5)
box(lwd=2, col='black')
points( x=res_Ddx5_ctcf_Deseq2_RPGC$log2FoldChange[res_Ddx5_ctcf_Deseq2_RPGC$col=="#049494"],
        y=-log10(res_Ddx5_ctcf_Deseq2_RPGC$pvalue[res_Ddx5_ctcf_Deseq2_RPGC$col=="#049494"]),
        col="#049494")
points( x=res_Ddx5_ctcf_Deseq2_RPGC$log2FoldChange[res_Ddx5_ctcf_Deseq2_RPGC$col=="#305494"],
        y=-log10(res_Ddx5_ctcf_Deseq2_RPGC$pvalue[res_Ddx5_ctcf_Deseq2_RPGC$col=="#305494"]),
        col="#305494")

table(res_Ddx5_ctcf_Deseq2_RPGC$col)

## 
## #049494 #305494  gray90 
##     932    1099   40201

df = data.frame( LFC=res_Ddx5_ctcf_Deseq2_RPGC$log2FoldChange,
                 P=res_Ddx5_ctcf_Deseq2_RPGC$log2FoldChange)
write.table(df,file=paste0(source_data_directory,"MA_ddx5KO.txt"),row.names = FALSE, sep="\t",quote=FALSE)

lost_peaks = rownames(res_Ddx5_ctcf_Deseq2_RPGC[res_Ddx5_ctcf_Deseq2_RPGC$padj<0.2 & res_Ddx5_ctcf_Deseq2_RPGC$log2FoldChange<0,])
gained_peaks = rownames(res_Ddx5_ctcf_Deseq2_RPGC[res_Ddx5_ctcf_Deseq2_RPGC$padj<0.2 & res_Ddx5_ctcf_Deseq2_RPGC$log2FoldChange>0,])

lost_Ddx5 = ctcf_wt_NS_peaks_resized[which(names(ctcf_wt_NS_peaks_resized) %in% lost_peaks)]
gained_Ddx5 = ctcf_wt_NS_peaks_resized[which(names(ctcf_wt_NS_peaks_resized) %in% gained_peaks)]

write.table( res_Ddx5_ctcf_Deseq2_RPGC,
             file=paste0(objects_directory,"res_Ddx5_ctcf_Deseq2_RPGC.txt"),
             sep='\t',quote=FALSE,row.names=FALSE,col.names=TRUE )

load(paste0(objects_directory,"CTCF_Ddx5ko_NS_A3_RPGC.RData"))
par(mfrow=c(1,2),mar=c(5,5,5,3))
plot( x=seq(-1000,1000,length.out=200),
      y=colMeans(a3_1_ctcf_rpgc_ctcf_AP[rownames(a3_1_ctcf_rpgc_ctcf_AP) %in% lost_peaks,]),
      ty="l", col="blue3", ylim=c(0,100), 
      ylab="RPGC", xlab="Distance from CTCF peak summit")
lines(x=seq(-1000,1000,length.out=200,lwd=2),
      y=colMeans(a3_2_ctcf_rpgc_ctcf_AP[rownames(a3_2_ctcf_rpgc_ctcf_AP) %in% lost_peaks,]),
      col="blue3",lwd=2 )

## Warning: W poleceniu 'seq.default(-1000, 1000, length.out = 200, lwd = 2)':
##  dodatkowy argument 'lwd' zostanie odrzucony

lines(x=seq(-1000,1000,length.out=200),
      y=colMeans(cb1_ctcf_rpgc_ctcf_AP[rownames(cb1_ctcf_rpgc_ctcf_AP) %in% lost_peaks,]),
      col="steelblue3",lwd=2 )
lines(x=seq(-1000,1000,length.out=200),
      y=colMeans(ce10_ctcf_rpgc_ctcf_AP[rownames(ce10_ctcf_rpgc_ctcf_AP) %in% lost_peaks,]),
      col="steelblue3",lwd=2 )


plot( x=seq(-1000,1000,length.out=200),
      y=colMeans(a3_1_ctcf_rpgc_ctcf_AP[rownames(a3_1_ctcf_rpgc_ctcf_AP) %in% gained_peaks,]),
      ty="l", col="blue3", ylim=c(0,100), ylab="RPGC", xlab="Distance from CTCF peak summit",lwd=2)
lines(x=seq(-1000,1000,length.out=200),
      y=colMeans(a3_2_ctcf_rpgc_ctcf_AP[rownames(a3_2_ctcf_rpgc_ctcf_AP) %in% gained_peaks,]),
      col="blue3",lwd=2 )
lines(x=seq(-1000,1000,length.out=200),
      y=colMeans(cb1_ctcf_rpgc_ctcf_AP[rownames(cb1_ctcf_rpgc_ctcf_AP) %in% gained_peaks,]),
      col="steelblue3",lwd=2 )
lines(x=seq(-1000,1000,length.out=200),
      y=colMeans(ce10_ctcf_rpgc_ctcf_AP[rownames(ce10_ctcf_rpgc_ctcf_AP) %in% gained_peaks,]),
      col="steelblue3",lwd=2 )

df = data.frame( location=seq(-1000,1000,length.out=200),
                 lost_wt1=colMeans(a3_1_ctcf_rpgc_ctcf_AP[rownames(a3_1_ctcf_rpgc_ctcf_AP) %in% lost_peaks,]),
                 lost_wt2=colMeans(a3_2_ctcf_rpgc_ctcf_AP[rownames(a3_2_ctcf_rpgc_ctcf_AP) %in% lost_peaks,]),
                 lost_cb1=colMeans(cb1_ctcf_rpgc_ctcf_AP[rownames(cb1_ctcf_rpgc_ctcf_AP) %in% lost_peaks,]),
                 lost_ce10=colMeans(ce10_ctcf_rpgc_ctcf_AP[rownames(ce10_ctcf_rpgc_ctcf_AP) %in% lost_peaks,]),
                 gained_wt1=colMeans(a3_1_ctcf_rpgc_ctcf_AP[rownames(a3_1_ctcf_rpgc_ctcf_AP) %in% gained_peaks,]),
                 gained_wt2=colMeans(a3_2_ctcf_rpgc_ctcf_AP[rownames(a3_2_ctcf_rpgc_ctcf_AP) %in% gained_peaks,]),
                 gained_cb1=colMeans(cb1_ctcf_rpgc_ctcf_AP[rownames(cb1_ctcf_rpgc_ctcf_AP) %in% gained_peaks,]),
                 gained_ce10=colMeans(ce10_ctcf_rpgc_ctcf_AP[rownames(ce10_ctcf_rpgc_ctcf_AP) %in% gained_peaks,]))
                 
write.table(df,file=paste0(source_data_directory,"AP_CTCF_wt_ddx5KO.txt"),row.names = FALSE, sep="\t",quote=FALSE)

DMSO1_ctcf_rpgc = import.bw(paste0(chipseq_directory,"ChIP_Seq_CTCF_11-24_MusMus_es-NPC_DDX5_FKBP_KI_CTCF-Cterm_HALO+DDX5-FKBP-KI_NPC_Ddx5-FKBP_KI_4F11_DMSO-Rep1_Rep_1_RPGC.bw"))
dTAG1_ctcf_rpgc = import.bw(paste0(chipseq_directory,"ChIP_Seq_CTCF_11-24_MusMus_es-NPC_DDX5_FKBP_KI_CTCF-Cterm_HALO+DDX5-FKBP-KI_NPC_Ddx5-FKBP_KI_4F11_dTAG13-Rep1_Rep_1_RPGC.bw"))
DMSO2_ctcf_rpgc = import.bw(paste0(chipseq_directory,"ChIP_Seq_CTCF_11-24_MusMus_es-NPC_DDX5_FKBP_KI_CTCF-Cterm_HALO+DDX5-FKBP-KI_NPC_Ddx5-FKBP_KI_4F11_DMSO-Rep2_Rep_1_RPGC.bw"))
dTAG2_ctcf_rpgc = import.bw(paste0(chipseq_directory,"ChIP_Seq_11-24_MusMus_es-NPC_DDX5_FKBP_KI_CTCF-Cterm_HALO+DDX5-FKBP-KI_NPC_Ddx5-FKBP_KI_4F11_dTAG13-Rep2_Rep_1_RPGC.bw"))
seqlevelsStyle(DMSO1_ctcf_rpgc) = "ucsc"
seqlevelsStyle(dTAG1_ctcf_rpgc) = "ucsc"
seqlevelsStyle(DMSO2_ctcf_rpgc) = "ucsc"
seqlevelsStyle(dTAG2_ctcf_rpgc) = "ucsc"


DMSO1_ctcf_rpgc_AP = getSignalInBins(ctcf_NS_ranges,DMSO1_ctcf_rpgc,1)
dTAG1_ctcf_rpgc_AP = getSignalInBins(ctcf_NS_ranges,dTAG1_ctcf_rpgc,1)
DMSO2_ctcf_rpgc_AP = getSignalInBins(ctcf_NS_ranges,DMSO2_ctcf_rpgc,1)
dTAG2_ctcf_rpgc_AP = getSignalInBins(ctcf_NS_ranges,dTAG2_ctcf_rpgc,1)
rownames(DMSO1_ctcf_rpgc_AP) = ctcf_wt_NS_peaks_summit$name
rownames(dTAG1_ctcf_rpgc_AP) = ctcf_wt_NS_peaks_summit$name
rownames(DMSO2_ctcf_rpgc_AP) = ctcf_wt_NS_peaks_summit$name
rownames(dTAG2_ctcf_rpgc_AP) = ctcf_wt_NS_peaks_summit$name

save(DMSO1_ctcf_rpgc_AP,
     dTAG1_ctcf_rpgc_AP,
     DMSO2_ctcf_rpgc_AP,
     dTAG2_ctcf_rpgc_AP,
     file=paste0(objects_directory,"CTCF_4F11_NS_RPGC.RData"))

load(paste0(objects_directory,"CTCF_4F11_NS_RPGC.RData"))
a=90;b=110
c4f11_ctcf = data.frame( dmso1 = ( rowSums(DMSO1_ctcf_rpgc_AP[,a:b])),
                         dmso2 = ( rowSums(DMSO2_ctcf_rpgc_AP[,a:b])),
                         dTAG1 = ( rowSums(dTAG1_ctcf_rpgc_AP[,a:b])),
                         dTAG2 = ( rowSums(dTAG2_ctcf_rpgc_AP[,a:b])),
                         row.names = rownames(DMSO1_ctcf_rpgc_AP))

Check with raw signal like Ddx5

load(paste0(objects_directory,'ctcf_NS_ranges.RData'))
DMSO1_ctcf_raw = import.bw(paste0(chipseq_directory,"ChIP_Seq_CTCF_11-24_MusMus_es-NPC_DDX5_FKBP_KI_CTCF-Cterm_HALO+DDX5-FKBP-KI_NPC_Ddx5-FKBP_KI_4F11_DMSO-Rep1_Rep_1_filtered_unnormalized.bw"))
dTAG1_ctcf_raw = import.bw(paste0(chipseq_directory,"ChIP_Seq_CTCF_11-24_MusMus_es-NPC_DDX5_FKBP_KI_CTCF-Cterm_HALO+DDX5-FKBP-KI_NPC_Ddx5-FKBP_KI_4F11_dTAG13-Rep1_Rep_1_filtered_unnormalized.bw"))
DMSO2_ctcf_raw = import.bw(paste0(chipseq_directory,"ChIP_Seq_CTCF_11-24_MusMus_es-NPC_DDX5_FKBP_KI_CTCF-Cterm_HALO+DDX5-FKBP-KI_NPC_Ddx5-FKBP_KI_4F11_DMSO-Rep2_Rep_1_filtered_unnormalized.bw"))
dTAG2_ctcf_raw = import.bw(paste0(chipseq_directory,"ChIP_Seq_11-24_MusMus_es-NPC_DDX5_FKBP_KI_CTCF-Cterm_HALO+DDX5-FKBP-KI_NPC_Ddx5-FKBP_KI_4F11_dTAG13-Rep2_Rep_1_filtered_unnormalized.bw"))
seqlevelsStyle(DMSO1_ctcf_raw) = "ucsc"
seqlevelsStyle(dTAG1_ctcf_raw) = "ucsc"
seqlevelsStyle(DMSO2_ctcf_raw) = "ucsc"
seqlevelsStyle(dTAG2_ctcf_raw) = "ucsc"

DMSO1_ctcf_raw_AP = getSignalInBins(ctcf_NS_ranges,DMSO1_ctcf_raw,1)
dTAG1_ctcf_raw_AP = getSignalInBins(ctcf_NS_ranges,dTAG1_ctcf_raw,1)
DMSO2_ctcf_raw_AP = getSignalInBins(ctcf_NS_ranges,DMSO2_ctcf_raw,1)
dTAG2_ctcf_raw_AP = getSignalInBins(ctcf_NS_ranges,dTAG2_ctcf_raw,1)
rownames(DMSO1_ctcf_raw_AP) = ctcf_wt_NS_peaks_summit$name
rownames(dTAG1_ctcf_raw_AP) = ctcf_wt_NS_peaks_summit$name
rownames(DMSO2_ctcf_raw_AP) = ctcf_wt_NS_peaks_summit$name
rownames(dTAG2_ctcf_raw_AP) = ctcf_wt_NS_peaks_summit$name

save(DMSO1_ctcf_raw_AP,
     dTAG1_ctcf_raw_AP,
     DMSO2_ctcf_raw_AP,
     dTAG2_ctcf_raw_AP,
     file=paste0(objects_directory,"CTCF_4F11_NS_RAW.RData"))

dTAGD_ctcf = data.frame( dmso1 = ( rowSums(DMSO1_ctcf_raw_AP[,a:b])),
                         dmso2 = ( rowSums(DMSO2_ctcf_raw_AP[,a:b])),
                         dtag1 = ( rowSums(dTAG1_ctcf_raw_AP[,a:b])),
                         dtag2 = ( rowSums(dTAG2_ctcf_raw_AP[,a:b])) )
metadata_dTAG = data.frame( treatment = c("dmso","dmso","dtag","dtag"),
                                  row.names = colnames(dTAGD_ctcf))


dTAG_ctcf_Deseq2 = DESeqDataSetFromMatrix(
 countData = dTAGD_ctcf,
 colData = metadata_dTAG,
 design = ~ treatment )

## converting counts to integer mode

## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors

dTAG_ctcf_Deseq2 = DESeq(dTAG_ctcf_Deseq2, fitType = 'local')

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

res_dTAG_ctcf_Deseq2 = results(dTAG_ctcf_Deseq2, 
                              contrast = c("treatment", "dtag", "dmso"))
summary( res_dTAG_ctcf_Deseq2 )

## 
## out of 42232 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0 (up)       : 287, 0.68%
## LFC < 0 (down)     : 256, 0.61%
## outliers [1]       : 0, 0%
## low counts [2]     : 0, 0%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

dtag_normalized = counts( dTAG_ctcf_Deseq2,normalized=TRUE )

heatscatter( x=res_dTAG_ctcf_Deseq2$baseMean,
             y=res_dTAG_ctcf_Deseq2$log2FoldChange, 
             pch=19, cex=0.1, add.contour = T, 
             xlim=c(0,1200),ylim=c(-2,2),
             xlab=expression(Log[2] ~ (Signal)),
             ylab=expression(Log[2] ~ (LFC)), 
             main="dTAG13/DMSO")
abline( h=0 )
axis(1,lwd=2)
axis(2,lwd=2)
box(col="black",lwd=2)

RNA binding proteins mediate the maturation of chromatin topology during differentiation

Aleksandra Pekowska, Misbah Abbas