DESeq question - Page 2

avm970 · 07-31-2010, 02:52 PM

Thank you Wolfgang Huber for you help, now the plot command is working but still I didnt get a proper plot.
Now I got a blank plot, I am not understanding y I got a blank plot instead of MvA plot
IS it because of the presence of "-Inf" in log2FoldChange column?

dgu · 08-11-2010, 03:32 PM

Hello Simon and Wolfgang -

Could any of you elaborate on the presence of variance outliers? I've done several DE tests and I've sometimes noticed values in resVarA and resVarB differ a lot for one specific gene (8x, 20x). The manual (pg 10) suggests to excercise caution during interpretation, but is there any rule-of-thumb to determine one of these values is large compared to the other?

Thanks in advance,
Daniel

afadda · 02-09-2011, 05:55 AM

Hi,
I have installed DESeq and all seems right:

> biocLite("DESeq")
Using R version 2.11.1, biocinstall version 2.6.9.
Installing Bioconductor version 2.6 packages:
[1] "DESeq"
Please wait...

Warning: unable to access index for repository http://brainarray.mbni.med.umich.edu...d/contrib/2.11
trying URL 'http://www.bioconductor.org/packages/2.6/bioc/bin/macosx/leopard/contrib/2.11/DESeq_1.0.6.tgz'
Content type 'application/x-gzip' length 1533404 bytes (1.5 Mb)
opened URL
==================================================
downloaded 1.5 Mb

but when I try to use
> cds <- newCountDataSet( countsTable, conds )
Error: could not find function "newCountDataSet"

what seems to be the problem?

Wolfgang Huber · 02-09-2011, 06:01 AM

Afadda,

did you type
library("DESeq")
to load the library?

Also, you're using an old version of R and - what's more of a problem - a rather old version of DESeq. I'd recommend to update.

Wolfgang

afadda · 02-09-2011, 10:17 AM

ok. I'm ignorant in R...
I don't seem to be able to get the newCountDataSet to work. I get this error:
Error in round(countData) : Non-numeric argument to mathematical function
I checked, and non of the cells in the table contain alphabets. But, the class of one of the columns is "factor". I don't know why is it so and how to get around this.
thanks,

afadda · 02-10-2011, 07:51 AM

how about this error

> cds= newCountDataSet( table, conds )
Error in if (any(round(countData) != countData)) stop("The countData is not integer.") :
missing value where TRUE/FALSE needed

here all the table columns are of class integer. if anyone would just give me a clue on what's wrong with my data

Simon Anders · 02-10-2011, 07:58 AM

Please send the output of the commands

Code:

head(table)

and

Code:

str(table)

. They will show you what exactly is in your "table" variable.

(By the way, "table" is not a good name for a variable, as it shadows the function "table" that you may need later.)

On the long run, it may also pay to learn a bit about R. Maybe read a few pages into the Introduction to R.

afadda · 02-10-2011, 09:35 AM

> head(table)
WT1 WT2 T
Tb927.1.10 0 0 0
Tb927.1.30 0 0 0
Tb927.1.40 0 0 0
Tb927.1.110 2260 2914 4825
Tb927.1.270 214 278 359
Tb927.1.280 205 331 589

> str(table)
'data.frame': 7882 obs. of 3 variables:
$ WT1: int 0 0 0 2260 214 205 628 26 379 142 ...
$ WT2: int 0 0 0 2914 278 331 724 58 591 261 ...
$ T : int 0 0 0 4825 359 589 1027 27 757 1380 ...

Simon Anders · 02-11-2011, 02:53 AM

Sorry, I didn't notice that you have missing values, according to your error message.

So, use the command

Code:

which( is.na(table), arr.ind=TRUE )

to find out which rows contain missing values ("NA"). To see, e.g., the content of rows 7, 12, and 17, write

Code:

table[ c(7, 12, 17), ]

(Note the extra comma!) This should give you a hint what went wrong.

Of course, you could simply remove the offending rows with:

Code:

cleanTable <- table[ -which( is.na(table), arr.ind=TRUE )[,"row"], ]

Simon

afadda · 02-11-2011, 07:03 AM

thanks Simon. it's working now!

masterpiece · 03-07-2011, 09:08 PM

Hi Simon,

I'm new to DE analysis. Just tried DESeq for my dataset. FYI, the dataset don't have replicate. When running libsizes

Code:

> libsizes <- c(C=6821199 , E=6037299)

when I checked the size factor, it show

Code:

> sizeFactors (cds)
C             E
1             1

Is the size factor for both sample '1' is normal. Is it due to working with no replicate dataset?

NicoBxl · 03-08-2011, 12:46 AM

Doing DE analysis without any replicate in any of the condition is a none sense. You can not have significant results without replicates.

For your question :

I tried with artificial data and it works well.

A <- data.frame(x=c(1000,10,10,10,0),y=c(100,2,2,2,1))
conds<-c("a","b")
cds <- newCountDataSet(A,conds)
cds <- estimateSizeFactors(cds)
sizeFactors(cds)
x y
2.2360680 0.4472136

masterpiece · 03-08-2011, 02:51 AM

Hi NicoBxl

I've tried plying with the example dataset. I treated only two sample (T1b and T2) trying to mimic my dataset as there is no replicate. I still got the sizeFactors not to '1' (as shown below). So I believed non replicate dataset is not an issue here.

Code:

> sizeFactors(cds)
     T1b       T2
0.604218 1.655032

Yes, now understand that without any replicate, its hard to get concrete conclusion on DE analysis. That is our mistakes here. Any how we try to get as much as possible input from our dataset.

Still not sure what is the prob, why my sizeFactor = 1 on both sample.

Simon Anders · 03-08-2011, 06:17 AM

Hi

Quote:

Originally Posted by masterpiece

Code:

> libsizes <- c(C=6821199 , E=6037299)

The vignette is a bit badly worded here. Ignore the part about the 'libsizes', you do not need this, and the information is overwritten, once you call 'estimateLibrarySizes'.

Have you done this and the got this here?

Quote:

Code:

> sizeFactors (cds)
C             E
1             1

Make a plot with

Code:

plot( counts(cds)[,1] + .01, counts(cds)[,2] + .01, log="xy" )

to see how your data looks.If you are still confused, show us the plot.

Simon

masterpiece · 03-08-2011, 06:52 AM

Hi Simon,

Thanks for the reply. Hope I got you right. Using my dataset I've follow exactly as instructed in the tutorial until I got

Code:

> sizeFactors (cds)
C             E
1             1

If i do not need to run libsizes, can I normalize my sample? About the 'estimateLibrarySizes'. Do i need to run them? because i can't see it in the DESeq

I've create the plot. Attached is the plot

Simon Anders · 03-08-2011, 07:11 AM

Quote:

Originally Posted by masterpiece

If i do not need to run libsizes, can I normalize my sample? About the 'estimateLibrarySizes'. Do i need to run them? because i can't see it in the DESeq

I don't quite get your last sentence, but anyway. Them minimal set of commands you always need are:

Code:

cds <- newCountDataSet( countsTable, conds )
cds <- estimateSizeFactors( cds )
cds <- estimateVarianceFunctions( cds )
res <- nbinomTest( cds, "T", "N")

Quote:

I've create the plot. Attached is the plot

You could have maybe made it a little larger... Nevertheless, with quite some squinting, it seems to me that your data really scatters quite exactly around the x=y diagonal, i.e., it might be close to equal sequencing depth. If so, 1 and 1 as size factors should be fine.

Simon

masterpiece · 03-08-2011, 07:54 AM

Hi Simon,

About the attachment, yup, tried to attach larger plot, but the maximum size attachment allowed is only 680x280. So the plot size has to be reduced. Sorry about that.

Great. At least now I have some confident to work with the data. Thanks a lot for your explanation.

Simon Anders · 03-11-2011, 05:10 AM

Hi masterpiece

I've got another report from somebody who got exactly 1 as size factors, and so had a closer look. The problem is that the size factor estimation tends to round values when you have many genes with only very few counts. So, I've just improved the method slightly and the newest version of DESeq, version 1.3.3, now offers a new optional argument to 'estimateSizeFactors'.

So, please try again, and now call estimateSizeFactors as follows:

Code:

library( genefilter )
cds <- estimateSizeFactors( cds, locfunc=shorth )

See help page to 'estimateSizeFactorsForMatrix' for details.

Version 1.3.3 is available via SVN now, and for download in a day or so.

Simon

masterpiece · 03-14-2011, 03:21 AM

Hi Simon,

It does make sense. Our data is from DGE sample and we dont do any count filtering on the data. I can say most of the sequence have very few count. I'll try redo the analysis.

Thanks again

KellerMac · 06-03-2011, 11:41 AM

Great program for processing data. I used it to find fold change and which genes had the greatest fold change between groups. I want to know if its possible to sort the data according to fold change and if I could view more than just the top six in the "head" display. Thank you much.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
DESeq: question about with replicates and without any replicates.	nb509	RNA Sequencing	2	10-25-2011 07:04 AM
DESeq: question about baseMean	Azazel	Bioinformatics	2	05-22-2011 06:40 PM
DESeq: question about baseMean. Also, replicates.	Azazel	Bioinformatics	5	05-18-2011 11:51 PM
DESeq question	gfmgfm	Bioinformatics	2	04-18-2011 04:15 AM
Another DESeq question	shurjo	Bioinformatics	2	05-15-2010 10:25 PM

08-11-2010, 03:32 PM	#22
dgu Junior Member Location: Seattle Join Date: Jun 2010 Posts: 5	Variance outliers? Columns resVarA or resVarB Hello Simon and Wolfgang - Could any of you elaborate on the presence of variance outliers? I've done several DE tests and I've sometimes noticed values in resVarA and resVarB differ a lot for one specific gene (8x, 20x). The manual (pg 10) suggests to excercise caution during interpretation, but is there any rule-of-thumb to determine one of these values is large compared to the other? Thanks in advance, Daniel

02-09-2011, 06:01 AM	#24
Wolfgang Huber Senior Member Location: Heidelberg, Germany Join Date: Aug 2009 Posts: 109	Afadda, did you type library("DESeq") to load the library? Also, you're using an old version of R and - what's more of a problem - a rather old version of DESeq. I'd recommend to update. Wolfgang __________________ Wolfgang Huber EMBL

02-10-2011, 07:58 AM	#27
Simon Anders Senior Member Location: Heidelberg, Germany Join Date: Feb 2010 Posts: 994	Please send the output of the commands Code: head(table) and Code: str(table) . They will show you what exactly is in your "table" variable. (By the way, "table" is not a good name for a variable, as it shadows the function "table" that you may need later.) On the long run, it may also pay to learn a bit about R. Maybe read a few pages into the Introduction to R.

02-11-2011, 02:53 AM	#29
Simon Anders Senior Member Location: Heidelberg, Germany Join Date: Feb 2010 Posts: 994	Sorry, I didn't notice that you have missing values, according to your error message. So, use the command Code: which( is.na(table), arr.ind=TRUE ) to find out which rows contain missing values ("NA"). To see, e.g., the content of rows 7, 12, and 17, write Code: table[ c(7, 12, 17), ] (Note the extra comma!) This should give you a hint what went wrong. Of course, you could simply remove the offending rows with: Code: cleanTable <- table[ -which( is.na(table), arr.ind=TRUE )[,"row"], ] Simon

07-31-2010, 02:52 PM	#21
avm970 Junior Member Location: Birmingham, UK Join Date: Jul 2010 Posts: 5	Thank you Wolfgang Huber for you help, now the plot command is working but still I didnt get a proper plot. Now I got a blank plot, I am not understanding y I got a blank plot instead of MvA plot IS it because of the presence of "-Inf" in log2FoldChange column?

02-09-2011, 05:55 AM	#23
afadda Member Location: Qatar Join Date: Jun 2010 Posts: 15	Hi, I have installed DESeq and all seems right: > biocLite("DESeq") Using R version 2.11.1, biocinstall version 2.6.9. Installing Bioconductor version 2.6 packages: [1] "DESeq" Please wait... Warning: unable to access index for repository http://brainarray.mbni.med.umich.edu...d/contrib/2.11 trying URL 'http://www.bioconductor.org/packages/2.6/bioc/bin/macosx/leopard/contrib/2.11/DESeq_1.0.6.tgz' Content type 'application/x-gzip' length 1533404 bytes (1.5 Mb) opened URL ================================================== downloaded 1.5 Mb but when I try to use > cds <- newCountDataSet( countsTable, conds ) Error: could not find function "newCountDataSet" what seems to be the problem?

02-09-2011, 10:17 AM	#25
afadda Member Location: Qatar Join Date: Jun 2010 Posts: 15	ok. I'm ignorant in R... I don't seem to be able to get the newCountDataSet to work. I get this error: Error in round(countData) : Non-numeric argument to mathematical function I checked, and non of the cells in the table contain alphabets. But, the class of one of the columns is "factor". I don't know why is it so and how to get around this. thanks,

02-10-2011, 07:51 AM	#26
afadda Member Location: Qatar Join Date: Jun 2010 Posts: 15	how about this error > cds= newCountDataSet( table, conds ) Error in if (any(round(countData) != countData)) stop("The countData is not integer.") : missing value where TRUE/FALSE needed here all the table columns are of class integer. if anyone would just give me a clue on what's wrong with my data

02-10-2011, 09:35 AM	#28
afadda Member Location: Qatar Join Date: Jun 2010 Posts: 15	> head(table) WT1 WT2 T Tb927.1.10 0 0 0 Tb927.1.30 0 0 0 Tb927.1.40 0 0 0 Tb927.1.110 2260 2914 4825 Tb927.1.270 214 278 359 Tb927.1.280 205 331 589 > str(table) 'data.frame': 7882 obs. of 3 variables: $ WT1: int 0 0 0 2260 214 205 628 26 379 142 ... $ WT2: int 0 0 0 2914 278 331 724 58 591 261 ... $ T : int 0 0 0 4825 359 589 1027 27 757 1380 ...

02-11-2011, 07:03 AM	#30
afadda Member Location: Qatar Join Date: Jun 2010 Posts: 15	thanks Simon. it's working now!

03-07-2011, 09:08 PM	#31
masterpiece Member Location: malaysia Join Date: Mar 2009 Posts: 40	Hi Simon, I'm new to DE analysis. Just tried DESeq for my dataset. FYI, the dataset don't have replicate. When running libsizes Code: > libsizes <- c(C=6821199 , E=6037299) when I checked the size factor, it show Code: > sizeFactors (cds) C E 1 1 Is the size factor for both sample '1' is normal. Is it due to working with no replicate dataset?

03-08-2011, 12:46 AM	#32
NicoBxl not just another member Location: Belgium Join Date: Aug 2010 Posts: 264	Doing DE analysis without any replicate in any of the condition is a none sense. You can not have significant results without replicates. For your question : I tried with artificial data and it works well. A <- data.frame(x=c(1000,10,10,10,0),y=c(100,2,2,2,1)) conds<-c("a","b") cds <- newCountDataSet(A,conds) cds <- estimateSizeFactors(cds) sizeFactors(cds) x y 2.2360680 0.4472136

03-08-2011, 02:51 AM	#33
masterpiece Member Location: malaysia Join Date: Mar 2009 Posts: 40	Hi NicoBxl I've tried plying with the example dataset. I treated only two sample (T1b and T2) trying to mimic my dataset as there is no replicate. I still got the sizeFactors not to '1' (as shown below). So I believed non replicate dataset is not an issue here. Code: > sizeFactors(cds) T1b T2 0.604218 1.655032 Yes, now understand that without any replicate, its hard to get concrete conclusion on DE analysis. That is our mistakes here. Any how we try to get as much as possible input from our dataset. Still not sure what is the prob, why my sizeFactor = 1 on both sample.

03-08-2011, 07:54 AM	#37
masterpiece Member Location: malaysia Join Date: Mar 2009 Posts: 40	Hi Simon, About the attachment, yup, tried to attach larger plot, but the maximum size attachment allowed is only 680x280. So the plot size has to be reduced. Sorry about that. Great. At least now I have some confident to work with the data. Thanks a lot for your explanation.

03-11-2011, 05:10 AM	#38
Simon Anders Senior Member Location: Heidelberg, Germany Join Date: Feb 2010 Posts: 994	Hi masterpiece I've got another report from somebody who got exactly 1 as size factors, and so had a closer look. The problem is that the size factor estimation tends to round values when you have many genes with only very few counts. So, I've just improved the method slightly and the newest version of DESeq, version 1.3.3, now offers a new optional argument to 'estimateSizeFactors'. So, please try again, and now call estimateSizeFactors as follows: Code: library( genefilter ) cds <- estimateSizeFactors( cds, locfunc=shorth ) See help page to 'estimateSizeFactorsForMatrix' for details. Version 1.3.3 is available via SVN now, and for download in a day or so. Simon

03-14-2011, 03:21 AM	#39
masterpiece Member Location: malaysia Join Date: Mar 2009 Posts: 40	Hi Simon, It does make sense. Our data is from DGE sample and we dont do any count filtering on the data. I can say most of the sequence have very few count. I'll try redo the analysis. Thanks again

06-03-2011, 11:41 AM	#40
KellerMac Member Location: Baton Rouge, Louisiana Join Date: Jun 2011 Posts: 11	Viewing processed data Great program for processing data. I used it to find fold change and which genes had the greatest fold change between groups. I want to know if its possible to sort the data according to fold change and if I could view more than just the top six in the "head" display. Thank you much.