SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
DESeq2 Simon Anders Bioinformatics 123 07-06-2015 01:45 AM
DESeq V DESeq2 sebastion RNA Sequencing 35 10-16-2014 06:04 AM
Tenure-track job in Bioinformatics/Computational Biol. at Auburn Univ. KMH Academic/Non-Profit Jobs 0 12-03-2012 11:50 AM
DESeq: question about with replicates and without any replicates. nb509 RNA Sequencing 2 10-25-2011 06:04 AM

Reply
 
Thread Tools
Old 02-24-2016, 10:14 AM   #21
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

There's nothing to normalize before mapping.
dpryan is offline   Reply With Quote
Old 02-23-2019, 08:40 AM   #22
FJwlf
Junior Member
 
Location: Spain

Join Date: Feb 2019
Posts: 2
Default

Quote:
Originally Posted by Simon Anders View Post
To be honest, we couldn't yet be bothered to explain how to analyse such data in DESeq2. It's tricky to write up because too many people will misinterpret whatever I write as if it were actually possible to conduct a meaningful statistical analysis when comparing just two samples.

So, if you promise to not use any such comparisons for actual science, here is how you do it:

Start as above:

Code:
library(DESeq2)
library(pasilla)
data("pasillaGenes")
countData <- counts(pasillaGenes)
countData<-countData[,c("treated1fb","untreated1fb")]
colData <- pData(pasillaGenes)[c("treated1fb","untreated1fb"),c("condition","type")]
dds <- DESeqDataSetFromMatrix(
       countData = countData,
       colData = colData,
       design = ~ condition)
Now use DESeq2's new "rlog transformation". This replaced the VST we had before. It transforms the average of the genes across samples to a log2 scale but "pulls in" those genes for which the evidence for strong fold changes is weak due to low counts.

Code:
rld <- rlogTransformation( dds )
As this is a logarithm-like scale, the differences between the samples can be considered as a "regularized log2 fold change". Let's make a result data frame:
Code:
res <- data.frame(
   assay(rld), 
   avgLogExpr = ( assay(rld)[,2] + assay(rld)[,1] ) / 2,
   rLogFC = assay(rld)[,2] - assay(rld)[,1] )
And now we have a ranking of genes by regularized fold changes:

Code:
> head( res[ order(res$rLogFC), ] )
            treated1fb untreated1fb avgLogExpr    rLogFC
FBgn0011260   7.830359     6.627326   7.228842 -1.203033
FBgn0001226  10.128636     8.929985   9.529311 -1.198652
FBgn0034718   8.503006     7.315640   7.909323 -1.187366
FBgn0003501   7.927864     6.743974   7.335919 -1.183889
FBgn0033635  11.126300     9.973979  10.550139 -1.152321
FBgn0033367  13.411814    12.269436  12.840625 -1.142378
This ranking put ones genes on top which are strongly downregulated in the second sample compared to the first one. If you do this with normal log expressions, the weak genes will appear at the top because they are noisiest and hence tend to have exaggerated fold changes.

The advantage of this procedure is that it does not produce any p values (which would be misleading anyway).
Hi everyone!
just a silly question... why is a (-) and not a division (/) to calculate the rlogFC;

res <- data.frame(
assay(rld),
avgLogExpr = ( assay(rld)[,2] + assay(rld)[,1] ) / 2,
rLogFC = assay(rld)[,2] / assay(rld)[,1] )

Thank you!!
FJwlf is offline   Reply With Quote
Old 02-23-2019, 10:17 AM   #23
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

Quote:
Originally Posted by FJwlf View Post
Hi everyone!
just a silly question... why is a (-) and not a division (/) to calculate the rlogFC;

res <- data.frame(
assay(rld),
avgLogExpr = ( assay(rld)[,2] + assay(rld)[,1] ) / 2,
rLogFC = assay(rld)[,2] / assay(rld)[,1] )

Thank you!!
The values are already log scale and log(a) - log(b) is the same as log(a/b).
dpryan is offline   Reply With Quote
Old 02-24-2019, 08:01 AM   #24
FJwlf
Junior Member
 
Location: Spain

Join Date: Feb 2019
Posts: 2
Default

Quote:
Originally Posted by dpryan View Post
The values are already log scale and log(a) - log(b) is the same as log(a/b).
Thank you!
FJwlf is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 12:15 PM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2019, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO