Hi, I was trying to investigate the changes of H3K27ac upon conditional KO of gene X.

After mapping and counting the data using "featurecount" from "Rsubread" packages, I used the DEseq2 to perform differential analysis. 

Below is the scatter plot (WT vs CKO), I could obtain genomic region with the significantly changed H3K27ac (as highlighted in "red" in the plot), 

but I also noticed that there is global shift of histone marks.

One possibility I could think of is that this shift might just reflect the fact that the sequence depth is different between WT and CKO samples. 

I checked the "regularized log transformation" part of the vigentte as well as the original paper of "DESeq2",

To be honest, I was not able to fully understand the underlying mathematics, but it seems me that the "rlog" transformation has factored in "size factor", or  the sequence depth . Even though the main rationale of performing "rlog" transformation is to stabilize the variance across gene with difference means. 

So does it mean the shift did reflect the global reduction of H3K27ac upon gene X KO,  given that the wet experiments and other data process performed properly?  ( I guess even it is a correct interpretation, other independent method might be needed to verify the conclusion.) 

Or did I just misunderstand of "rlog" methodology? 

Thanks in advance.

I asked the question at biostar, in a more general sense.

https://www.biostars.org/p/255052/

 Devon Ryan  kindly suggested to use DiffBind to "normalize only to non-peak regions", I checked "DiffBind", could not find the options he referred to. But maybe I missed some key points. So if "Diffbind" could help solve my issue, please kindly let me know. 

if there is global up or down changes you need to specify the controlGenes for the size factor estimation to make sense.

Yes there can be changes when you update versions. If you need the results not to change, best to stay within a release and not change the version of R/Bioconductor.

I don't follow what is being plotted. Is that rlog-ed values? The rlog transformation does correct for sequencing depth, by using already calculated size factors or estimating them if they have not been calculated. This means you can provide whatever size factors you like (e.g. estimated from subsets of regions) before calling rlog().

I'm not sure about the normalization options within DiffBind, but I do know that you'll need to provide your code for someone to offer useful suggestions on that part.

@Michael Love, I have posted the question related to DiffBind, has yet to receive any feedback. but I have found some related posts, 

in replying " DiffBind dba.report() output: Is "concentration" normalized?"

Rory explained, "the counts are normalized. Concentrations are reported as log2 values."

in term of "normalization factor",

by default (bFullLibrarySize=TRUE), 

"the factors are set to:

> DESeq2::sizeFactors(DESeqDataSeq) <- libsize/min(libsize)"

the sizeFactors ill set as library size (total reads numbers in the bam file) 

otherwise(bFullLibrarySize=FALSE), 

then "calls DESeq2::estimateSizeFactors() to calculate the normalization factors. 

I changed my argument when apply DiffBind analysis, basically use the  "size factor" estimated by DESeq2

for normalization.

H3K27ac_CellA_CellD <- dba(sampleSheet="sampleconfiguration.csv",)
H3K27ac_CellA_CellD_Count <- dba.count(H3K27ac_CellA_CellD, score= DBA_SCORE_TMM_READS_EFFECTIVE, mapQCth=0, bRemoveDuplicates=FALSE) 
H3K27ac_CellA_CellD_Diff <- dba.analyze(H3K27ac_CellA_CellD_Count,method=DBA_DESEQ2,　bFullLibrarySize=FALSE)
H3K27ac_CellA_CellD_Diff_All <- dba.report(H3K27ac_CellA_CellD_Diff,method=DBA_DESEQ2, th=1) 

H3K27ac_All <- data.frame (Conc_WT=H3K27ac_CellA_CellD_Diff_All @elementMetadata$Conc_WT,
                           Conc_GeneXcKO= H3K27ac_CellA_CellD_Diff_All @elementMetadata$Conc_GeneXcKO)


p <- ggplot(H3K27ac_All, aes(Conc_WT, Conc_GeneXcKO))
p + geom_point() +
    geom_abline(intercept = 0, slope = 1, linetype=2) +
    scale_x_continuous(name="Epifeatures\n in WT Cell ", limits = c(-1, 11)) +
    scale_y_continuous(name="Epifeatures \n in KO Cell ",limits = c(-1, 11)) 

The plot looks right, in the code above, raw read counts were normalized against size factor calculated by "DESeq2", then log transformed by the base of 2.  It seems to me that the value is arithmatically equivalent to rlog() in DESeq2. ( I am not 100% sure, the original DEseq2 paper explained in statistic notation with the statistics model in mind, which I am not sure that I understood it completely.)

If it is "regularized log transformation" result, I wonder is there any possibility "rlog" might mask the global differences among the samples.

by that, I do not mean rlog() itself are flawed, but is there any occasion "rlog" (or normalize by "size factor" ) should not be applied? 

in the original "DEseq" paper it explained why total read counts might not be a good measure of sequence depth: "a few highly and differentially expressed genes may have strong influence on the total read count, causing the ratio of total read counts not to be a good estimate for the ratio of expected counts.",   it seems to me that it implicitly assumes that most or majority of genes, expression of them are largely comparable (if plotted, distributed alone the line with slope of 1) , does this means it should be cautious when there is possibility that all the genes are disregulated in one direction (up or down) ?  I guess such scenario should be quite rare.  I guess this might be beyond the scope of DESeq2, but if suspect such a scenario (global increase or decrease), which might be the recommended way to normalize? 

PS: I forgot to mention that when I re ran my previous code, I got the different report of differential analysis,  as far as I am aware, I used the same codes. The only thing I could think of is that I reinstalled DESeq2 this week, so the versions might be different, could that contribute to the differences here?

1 Contrast:
   Group1 Members1 Group2 Members2 DB.DESeq2
1 KO        4     WT        4       133

1 Contrast:
   Group1 Members1 Group2 Members2 DB.DESeq2
1 KO        4     WT        4        64