I am wondering if there legit or if we can find a legit way to adapt RUVSeq -- RUVg approach to SNAP-ChIP assays.

As described on the page listed below, spike-ins (nucleosomes with a synthetic sequence) are typically added to the samples before immuno-precipitation step in such a way that global changes in the histone modifications can be captured by measuring the spike-in intensity.

In our case, we do use SNAP-ChIP spikes on 4 samples : 2 samples that are treated with DMSO, and 2 samples treated with ZT.

The R code that I am using is :

#################################################################################
#################################################################################

library("EDASeq")
library("RUVSeq")
library("RColorBrewer")

spikeins = data.frame(
DMSO_1 = 56955 ,
DMSO_7  = 56765 ,
ZT_2 = 24788 ,
ZT_8 = 47588 
)

rownames(spikeins) = "spikeins"


COUNTS = read.table("INPUT_FILE.txt", 
                         header=T, sep="\t", stringsAsFactors=F) 

COUNTS_dse = rbind(spikeins, COUNTS)

# after the filtering step that follows, and : we work with the data called "filtered"

group <- factor(c("DMSO", "DMSO", "ZT", "ZT"))
group <- relevel(group, ref="DMSO")

#################################################################################
#################################################################################

set <- newSeqExpressionSet(as.matrix(filtered),
                           phenoData = data.frame(group, row.names=colnames(filtered)))

phenoData(set)
featureData(set)

spikes <- rownames(filtered)[grep("^spike", rownames(filtered))]
spikes

loci <- rownames(filtered)[grep("-", rownames(filtered))]
loci

phenoData(set)
featureData(set)

# I understand that I should skip betweenLaneNormalization ()
# set <- betweenLaneNormalization(set, which="upper") 

#################################################################################
#################################################################################

set1 <- RUVg(set, spikes, k=1)
pData(set1)

              x        W_1
DMSO_1 DMSO -0.6887355
DMSO_7 DMSO -0.2536409
TAZ_2   ZT  0.3782554
TAZ_8   ZT  0.5641210

#################################################################################
#################################################################################
####### computing the differential expression using edgeR
# design <- model.matrix(~group + W_1, data=pData(set1))

design <- model.matrix((~0 + group + W_1), data=pData(set1))

contrast.matrix <- makeContrasts(ZT_DMSO=ZT-DMSO, levels=design)

y <- DGEList(counts=counts(set1))
y <- calcNormFactors(y, method="upperquartile").  #  shall we keep the computation of the normalization factors ?
y <- estimateGLMCommonDisp(y, design)
y <- estimateGLMTagwiseDisp(y, design)

fit <- glmFit(y, design)

lrt <- glmLRT(fit, contrast = contrast.matrix[,"ZT_DMSO"]) 

topTags(lrt)
dim(subset(as.data.frame(lrt), (logFC < -1) & (PValue < 0.01 )))[1]

However, the results that I am getting show the same number of up-reg or down-reg regions (approx 10 000 up, and approx 10 000 down), when I should have expected to see a global decrease in the amount of the histone mark (100 000 loci) (as confirmed by regular ChIP)

Thanks,

Bogdan

ps : the description of the spike-in protocol

https://www.thermofisher.com/us/en/home/references/newsletters-and-journals/bioprobes-journal-of-cell-biology-applications/bioprobes-79/snap-chip-histone-PTM-antibody-specificity.html

```