Hi Michael and Simon,

Here is my complete script for DESeq2:

#open DESeq
library("DESeq2")
setwd("~/Dropbox/TESIS/DGEA_2.0/Chapter_4_thesis/DE_analysis/DESeq2")

#Count matrix input
matrix = read.table (file.choose(), header=TRUE, row.names=1)


CeleDesign <- data.frame(
  row.names = colnames(matrix),
  condition = factor(c("SP", "SP", "EB", "EB")))
dds <- DESeqDataSetFromMatrix(countData = matrix,colData = CeleDesign,design
= ~ condition)
dds

#Est size factor = normalize for library size
dds<- estimateSizeFactors(dds)
dds <- estimateDispersions(dds)
dds <- nbinomWaldTest(dds)

#plot dispersion
plotDispEsts(ddsLocal, ylim = c(1e-6, 1e1))
plotDispEsts(dds, ylim = c(1e-6, 1e1))

#Differential expression analysis
resultsNames(dds)
res <- results(dds)
res <- res[order(res$padj),]
head(res)
plotMA(dds,ylim=c(-3,3),main="DESeq2")
sum(res$padj < .1, na.rm=TRUE)
mcols(res, use.names=TRUE) ##To know the level of conditions

#filter for significant genes
resSig <- res[ which(res$padj < .1), ]
head(resSig)

#filter for upregulated and downregulated genes
##up-regulated
tail(resSig[order(resSig$log2FoldChange), ])
##down-regulated
head(resSig[order(resSig$log2FoldChange), ])
##list most sig genes
head( resSig[ order(resSig$padj), ] )

write.table(as.data.frame(resSig [ order(resSig$padj), ]),
file="MostSigGenes.txt")
write.table(as.data.frame(resSig[ order(
resSig$log2FoldChange,-resSig$baseMean),]), file="DownRegulated.txt")
write.table(as.data.frame(resSig[ order( -resSig$log2FoldChange,
-resSig$baseMean ),]), file="UpRegulated.txt")


   1. MostSigGenes.txt, DownRegulated.txt and UpRegulated.txt always have
   the same number of genes. This is confusing for me.
   2. I understand then that maybe is better to only generate
   MostSigGenes.txt and then divide the result in 2 txt files: one with
   positive log2FC = upregulated and the other with negative log2FC =
   downregulated, right?
   3. I don't understand with if I am entering the conditions like "condition
   = factor(c("SP", "SP", "EB", "EB")))" then is comparing EB vs. SP. I
   would have thought the opposite.


Thanks,
Alicia



--
Alicia R. Pérez-Porro
PhD candidate

Giribet lab
Department of Organismic and Evolutionary Biology
MCZ labs
Harvard University
26 Oxford St, Cambridge MA 02138
phone: +1 617-496-5308
fax: +1 617-495-5667
www.oeb.harvard.edu/faculty/giribet/

Department of Marine Ecology
Center for Advanced Studies of Blanes (CEAB-CSIC)
C/Accés Cala St. Francesc 14
17300 Blanes, Girona, SPAIN
phone: +34 972 336 101
fax: +34 972 337 806
www.ceab.csic.es


On Mon, May 26, 2014 at 5:01 AM, Simon Anders <anders@embl.de> wrote:

> Hi Alicia
>
>
> On 26/05/14 01:39, Alicia R. Pérez-Porro wrote:
>
>> mcols(res, use.names=TRUE) ##To know the level of conditions
>>>
>> DataFrame with 6 rows and 2 columns
>>                         type                                description
>>                  <character>                                <character>
>> baseMean       intermediate                the base mean over all rows
>> log2FoldChange      results log2 fold change (MAP): condition SP vs EB
>> lfcSE               results         standard error: condition SP vs EB
>> stat                results         Wald statistic: condition SP vs EB
>> pvalue              results      Wald test p-value: condition SP vs EB
>> padj                results                       BH adjusted p-values
>>
>> By exporting my up-regulated genes doing:
>>
>> write.table(as.data.frame(resSig[ order( -resSig$log2FoldChange,
>> -resSig$baseMean ),]),
>>            file="UpRegulated.txt")
>>
>> Are my genes up-regulated in the EB condition? or in the SP condition?
>>
>
> Have a look at the column description: It says "log2 fold change (MAP):
> condition SP vs EB".
>
> With "SP vs EB", we mean that SP was compared against EB. Hence: positive
> sign means stronger in SP than in EB.
>
>
> > What about this example: (generated with
> > write.table(as.data.frame(resSig[
> > order( -resSig$log2FoldChange, -resSig$baseMean ),]),
> >          file="UpRegulated.txt")
>
> Please show us how you generated 'resSig' from 'res', and how tyou got
> 'res' from your DESeqDataset object.
>
> BTW, if you can your e-mail program to not wrap the lines, this would make
> it much more readable.
>
>   Simon
>
>
> _______________________________________________
> Bioconductor mailing list
> Bioconductor@r-project.org
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives: http://news.gmane.org/gmane.
> science.biology.informatics.conductor
>

	[[alternative HTML version deleted]]

