Entering edit mode
green.rebeccam
•
0
@greenrebeccam-9713
Last seen 9.0 years ago
Hi All, I am working with an RNAseq data set of 20 samples across 4 geneotypes with increasing severity of the phenotype across the genotypes. One of the questions we are trying to answer is does the variance of gene expression (level of disregulation) change across the groups. I know that DeSeq2 calculates a gene-wise diserpersion estimate. Does anyone have any ideas to use this to calculate a sample dispersion level and to compare that to other samples? - Thanks!!!!
So would this be a correct set of code? I have run my original anaylsis in Deseq2. counts should have been generated from Deseq2 summarizeOverlaps. Correct? I see a lot of variance, but it isn't correlated with my genotypes and I am trying to make sure i didn't screw something up.
Then run this through limma
#generate overlaps
se <- summarizeOverlaps(features=ebg, reads=bamfiles, mode="Union", singleEnd=TRUE, ignore.strand=FALSE, preprocess.reads=complementStrand)
meta <- read.csv(“sample_meta_litter.tab")
(colData(se) <- DataFrame(meta))
dds<- DESeq(dds)
design <- model.matrix(se$Category)
voom <- voomWithQualityWeights(counts, design=design, plot=TRUE)
summary (lm(voom$weight ~ se$Category))
No, the code isn't correct. Lots of issues. The model.matrix() command will give an error, and the last line probably will also. What are you trying to achieve by regressing voom weights on Category? That doesn't seem a sensible thing to do, and the y and x vectors are different lengths in your command. Where does 'counts' come from? What is the DESeq() command for -- it doesn't seem to do anything here because you don't use the output? If you want to model different dispersions for different groups, you need to specify the var.design argument to voomWithQualityWeights.
I can't quite tell what you're trying to do. Why don't you just follow the usual sort of voom RNA-seq analysis?
I think this might have been intended as a reply to my answer, where I mentioned that I plotted and regressed the weights against covariates to see whether the heteroskedasticity seems related to groups/batches or whether it is more random.
Oh I see. I think this would be a better way to visualize it:
Is voomaByGroup different than using voomWithQualityWeights with the var.design argument?
The main difference in this context is that it will plot a separate voom curve for every category value, so you can see how different the curves are.
voomWithQualityWeights() can't plot different curves because it handles cases that a more general and can't be summarized as distinct curves. There are other differences, but they are perhaps minor in this context.
sorry, but with is the function cpm?
never mind, its part of edgeR... sorrt
Yes - I was trying to see if the heteroskedasticity had anything to do with groups. You are correct that I am not using the deseq here for this output, its used other places in my code. I am however, using the count table generated by the summarizeOverlaps.
You are right on the model.matrix, that was a type and missing the ~, but the matrix looks as expected.