size factors in DESeq

[highlei]

My case is:
We have a deep-sequencing for RNA: two samples (WT and mutation) and 3 replicates for each one.
These small RNA sequences were map to ref-genome.
The ref-genome was seperated into 500bp windows along the whole genome.
We wanna identify the windows with decreased express from WT to Mutation.

Firstly, I try to use two different methods to estimate the size factors.
1. I try to use total counts of reads to estimate the size factors.
p_adj < 0.05 53,866 windows were got.
2. I try to use the estimateSizeFactors to estimate the size factors.
p_adj < 0.05 4,555 windows were got.

Secondly, since the genome is long, I make a filter: the windows with total raw counts in all 6 libraries < 50 were discarded.
1. I try to use total counts of reads to estimate the size factors.
p_adj < 0.05 48,005 windows were got.
2. I try to use the estimateSizeFactors to estimate the size factors.
p_adj < 0.05 464 windows were got.
I donot know I should use the total counts or use estimateSizeFactors to estimate the size factors.

[Simon Anders]

To check whether your normalization is good, compare two samples with an MvA plot, as follows:

First, get the normalized counts:

ncounts <- counts( cds, normalized=TRUE )
# or, for the release version:
ncounts <- t( t( counts(cds) ) / sizeFactors(cds) )

The, plot the mean of, say samples 1 and 2, versus their ratio:

plot( 
   ( ncounts[,2] + ncounts[,1] )/2,
   ncounts[,2] / ncounts[,1],
   log="xy", pch="." )
abline( h = 1)

If the normalization is good, the bulk of the counts should lie on the horizontal line marking equal value in both samples.

You can also plot each sample against the mean of all the other samples:

means <- rowMeans( ncounts )
plot( 
   means,
   ncounts[,2] / means,
   log="xy", pch="." )
abline( h = 1)

Further remarks:

- "We have a deep-sequencing for RNA: two samples (WT and mutation) and 3 replicates for each one." I hope you mean that you have six samples, 3 WT and 3 mutation. if you really just have sequenced the same sample three times (as you wrote), you better start over.

- "We wanna identify the windows with decreased express from WT to Mutation." Do you have reason to believe that there is much more down-regulation than up-regulation? If so, you may need a custom normalization procedure (even though the new option "locfunc=shorth" in the devel version of DESeq might already help).

[highlei]

Thanks for your reply.

- "We have a deep-sequencing for RNA: two samples (WT and mutation) and 3 replicates for each one." I hope you mean that you have six samples, 3 WT and 3 mutation.
Yes. I have six samples.
- "We wanna identify the windows with decreased express from WT to Mutation." Do you have reason to believe that there is much more down-regulation than up-regulation?
Yes. We mutate a key protein in the pathway.

I have used your code to plot figures. But I don't know how to evaluate the figures.

6.ge50.1.jpg sizeFactors with the total counts.
6.ge50.2.jpg sizeFactors with the total counts.
6.ge50.SF.1.jpg sizeFactors with estimateSizeFactor.
6.ge50.SF.2.jpg sizeFactors with estimateSizeFactor.

[alittleboy]

To check whether your normalization is good, compare two samples with an MvA plot, as follows:

First, get the normalized counts:

ncounts <- counts( cds, normalized=TRUE )
# or, for the release version:
ncounts <- t( t( counts(cds) ) / sizeFactors(cds) )

The, plot the mean of, say samples 1 and 2, versus their ratio:

plot( 
   ( ncounts[,2] + ncounts[,1] )/2,
   ncounts[,2] / ncounts[,1],
   log="xy", pch="." )
abline( h = 1)

If the normalization is good, the bulk of the counts should lie on the horizontal line marking equal value in both samples.

You can also plot each sample against the mean of all the other samples:

means <- rowMeans( ncounts )
plot( 
   means,
   ncounts[,2] / means,
   log="xy", pch="." )
abline( h = 1)

Further remarks:

- "We have a deep-sequencing for RNA: two samples (WT and mutation) and 3 replicates for each one." I hope you mean that you have six samples, 3 WT and 3 mutation. if you really just have sequenced the same sample three times (as you wrote), you better start over.

- "We wanna identify the windows with decreased express from WT to Mutation." Do you have reason to believe that there is much more down-regulation than up-regulation? If so, you may need a custom normalization procedure (even though the new option "locfunc=shorth" in the devel version of DESeq might already help).

So just curious: when having 6 samples, to check the performance of normalization, shall I plot/check all pairwise MA-plot (e.g. 1 vs. 2, 1 vs. 3, etc.)? I am a little confused as I think normalization is done for different samples (liver vs. kidney, say) -- but in the code, ncounts[ ,1] refers to replicate 1 of the first sample?