Hi there,

I wish to compare differential enrichment of a histone modification (H3K9me2) across 5 replicates of control and treated samples, for each of two cell types. This modification is very broad and not amenable to peak analysis. Therefore I have used featurecounts to count up reads falling in gene-bodies for each gene and sample, and now have a counts matrix of K9me2 enrichment at each gene.

I would like to use this counts matrix in DESeq2 for differential enrichment analysis, as if treating it as RNA-seq in essence. However one issue I have is that for Cut&Run and for analysis of global changes in enrichment, I first normalize the data to the number of mapped reads generated from S.cerevisiae spike-in, within each sample. This means my data table now contains spike-in normalized values and not raw counts, which I know DESeq2 normally requires for DEseq(). If I round my values so that no error is thrown, is it OK to continue to use DEseq() to find genes with differential enrichment of the modification (and what would you recommend to set for the betaPrior?) If not advisable, is there a better way to allow normalization to spike-in values plus differential testing with DESeq2?

Many thanks for any help.

Hi,

I was recently very happy to find RUVSeq (https://bioconductor.org/packages/release/bioc/html/RUVSeq.html) for normalisation prior differential analysis. They have a method specifically developed for spike-in normalisation (RUVg), and provide further advice how to proceed with DEseq2.

Maybe this helps you. 

Best,

Alex

Do not normalize the data for sequencing depth and provide to DESeq2. This point is made frequently in our documentation.

You can instead do:

dds <- estimateSizeFactors(dds, controlGenes=...)
dds <- DESeq(dds)