This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Using the VST normalization rather than the DESeqDataSetFromMatrix function for DESeq2

Hi,

I do have a drosophila 2x150 illumina RNAseq data. It is a ph KD condition, involved in the Polycomb way, which is a transcriptional repressor. Then, when it is knock down, it will deregulate the genes expression.

To show that, I do a DESeq2 comparison, using a ph control. But there is an issue concerning the concordance between the bigwig coverage files and my log2FC. To show that, I will focus on a specific gene (Ago2) :

enter image description here

-> there are 3 controls (upper panel) and 3 KD (bottom panel), I did an autoscale of the coverage : we can see an increase of the coverage for the kd samples (5x vs 35x -> ratio 6).

Then, I run DESeq2 with a custom script, using the classic normalization with the DESeq() function . The problem I have concerns the size factors, which are not centered around. In a typical scenario, only a few genes are up- or down-regulated (~5–10%), so the size factors (i.e., the median of the ratios for each gene) are centered around 1. Normalization (raw count / size factor) has a negligible impact on the counts, and the log2FC of the up- or down-regulated genes is estimated correctly. In my experiment, several thousand genes are upregulated, which pulls the median of the ratios far from 1, with a size factor of ~1.4 for Ph-KD and ~0.7 for the control. The normalized counts are therefore artificially inflated for the control and underestimated for Ph-KD, which compresses the log2FC relative to reality. The log2FC of Ago2 is +1.28 with a padj of 2.25-e16 (x=2.42 , underestimated compared to ~6)

I do a new DESeq2 run , using a selection of the genes which are stable :

 # VST blind=TRUE 
  vsd     <- varianceStabilizingTransformation(dds, blind=TRUE)
  vst_mat <- assay(vsd)

  # variation coefficient on the VST values
  cv_vst <- apply(vst_mat, 1, function(x) sd(x) / mean(x))
  bm_vst <- rowMeans(counts_filtered)

  # raw log2FC on the raw counts (without normalization)
  mean_ctrl_vst <- rowMeans(counts_filtered[, ctrl_samples, drop=FALSE])
  mean_test_vst <- rowMeans(counts_filtered[, test_samples, drop=FALSE])
  lfc_brut_vst  <- log2((mean_test_vst + 1) / (mean_ctrl_vst + 1))

  # Double filters
  stable_cv  <- names(cv_vst[cv_vst < opt$cv_threshold & bm_vst >= opt$base_mean_min])
  stable_vst <- names(cv_vst[
    cv_vst  <  opt$cv_threshold  &
    bm_vst  >= opt$base_mean_min &
    abs(lfc_brut_vst) < opt$lfc_threshold
  ])

 # Save the stable features
   stable_df <- data.frame(
  feature_id = clean_feature_id(stable_vst),
  CV_VST     = round(cv_vst[stable_vst], 4),
  lfc_brut   = round(lfc_brut_vst[stable_vst], 3),
  baseMean   = round(bm_vst[stable_vst], 1)
 )

# Compute the size factors 
stable_counts_mat <- as.matrix(counts_filtered[stable_vst, ])
sf_stable <- estimateSizeFactorsForMatrix(stable_counts_mat)

 # DESeq2 with the new size factors
dds <- DESeqDataSetFromMatrix(countData = counts_filtered,
                              colData   = meta_subset,
                              design    = ~condition)
sizeFactors(dds) <- sf_stable
dds <- estimateDispersions(dds)
dds <- nbinomWaldTest(dds)

-> That way will correct the log2FC for AGO2 (and others similar genes) : 2,63 with a padj of 2.34e-52.

I don't know if it is a good way to do?

I know one solution would be to use some spike in , but for that sequencing, we did not.

Best

normalization deseq2 rnaseq
• 306 views • link
written by pailloufat • pailloufat

I don't think that the size factors need to be symmetric around one. In DESeq2 they scale with library size, so maybe one groups is just sequenced a bit deeper than the other. Rather than all this code and coverage plots, just run default analysis and show output of plotMA(). This will tell whether normalization is ok or not.

• 0 views • link
written by ATpoint ATpoint

Thanks for the reply! Here is the depth sequecing :

file                                             format  type    num_seqs        sum_len  min_len  avg_len  max_len
rawdata/No-ph-KD-rep1/No-ph-KD-rep1_R1.fastq.gz  FASTQ   DNA   19,110,016  2,866,502,400      150      150      150
rawdata/No-ph-KD-rep1/No-ph-KD-rep1_R2.fastq.gz  FASTQ   DNA   19,110,016  2,866,502,400      150      150      150
rawdata/No-ph-KD-rep2/No-ph-KD-rep2_R1.fastq.gz  FASTQ   DNA   20,444,393  3,066,658,950      150      150      150
rawdata/No-ph-KD-rep2/No-ph-KD-rep2_R2.fastq.gz  FASTQ   DNA   20,444,393  3,066,658,950      150      150      150
rawdata/No-ph-KD-rep3/No-ph-KD-rep3_R1.fastq.gz  FASTQ   DNA   22,128,177  3,319,226,550      150      150      150
rawdata/No-ph-KD-rep3/No-ph-KD-rep3_R2.fastq.gz  FASTQ   DNA   22,128,177  3,319,226,550      150      150      150
rawdata/Ph-KD-rep1/Ph-KD-rep1_R1.fastq.gz        FASTQ   DNA   22,091,755  3,313,763,250      150      150      150
rawdata/Ph-KD-rep1/Ph-KD-rep1_R2.fastq.gz        FASTQ   DNA   22,091,755  3,313,763,250      150      150      150
rawdata/Ph-KD-rep2/Ph-KD-rep2_R1.fastq.gz        FASTQ   DNA   21,164,573  3,174,685,950      150      150      150
rawdata/Ph-KD-rep2/Ph-KD-rep2_R2.fastq.gz        FASTQ   DNA   21,164,573  3,174,685,950      150      150      150
rawdata/Ph-KD-rep3/Ph-KD-rep3_R1.fastq.gz        FASTQ   DNA   22,942,903  3,441,435,450      150      150      150
rawdata/Ph-KD-rep3/Ph-KD-rep3_R2.fastq.gz        FASTQ   DNA   22,942,903  3,441,435,450      150      150      150

-> the coverage is quite homogeneous on the samples.

Here is the MA-plot. As you can see, the kd will affect massively the transcriptome :

maplot

• 0 views • link
written by pailloufat • pailloufat

0 answers

No answers yet.

Log in to answer this question.

More posts like this