Question: Filter by Fold change in DESeq2
5
gravatar for EpiExplorer
3.4 years ago by
EpiExplorer90
New Zealand
EpiExplorer90 wrote:

Hi,

I am using DESeq2 for RNA-Seq data. I got 557 genes deferentially expressed with 5% FDR and Log Fold change >= +/- 0.

I would like to choose the top genes that have a fold change >=2. How do I do that in DESeq2? Can I just export the list of 557 genes and filter them by fold change >=2 or do I have to use any argument in the results function to statistically calculate that?

The functions that I used in DESeq2 to get the 557 deferentially expressed genes is below:

res <- results(dds, alpha=0.05)

Thanks.

rna-seq • 4.8k views
ADD COMMENTlink modified 3.4 years ago by firatuyulur300 • written 3.4 years ago by EpiExplorer90

Hi, I need urgent help in this please. I have read the vignette but I am not clear how to get the genes that are significant with 5% FDR and also greater than or equal to fold change of -2/+2.

I have the list of genes which was extracted by results function with alpha value as 0.05. Now if I want to further filter this list based on the fold change, can I use the logfold change column in the file to filter the significant genes with -1/+1( Log(2 fold up or down) or do I need to necessarily rerun the statistics on fold change by specifying the lfsthreshold argument in the results function along with the alpha function set to 0.05? When I run results with only alpha, I get 557 genes significant res <- results(dds, alpha=0.05)

But when I run results with foldchange and alpha, i get 91 genes: res1 <- results(dds,alpha=0.05, lfcThreshold=1, altHypothesis="greaterAbs") By comparing the results with 557 genes and 91 genes, I found that the in list of 557 genes about only 180 genes are below +1/-1 logfold change. So I am wandering why I get only 91 genes when filtered by alpha and fold change both. I am sure DESeq2 is doing some internal statistics for p-value adjustment and fold change significance but not sure how and what it is doing.

Can somebody please explain me if I am running a correct function with right arguments for what I want as output?

ADD REPLYlink written 3.4 years ago by EpiExplorer90
1

Just a comment: in my opinion, the right approach is to run res <- results(dds, alpha=0.05) and then filter based on log2FoldChange (eventually remove the rows with NAs, which can appear if the number of reads in one condition equals zero). Thus I would assume that 180 genes are what you are looking for.

Now by specifying the altHypothesis parameter in res1 <- results(dds,alpha=0.05, lfcThreshold=1, altHypothesis="greaterAbs"), you calculate DEGs based on a different null hypothesis:

|beta| <= lfcThreshold

Per default, the null hypothesis is

beta = 0.

This explains, in my opinion, the difference in the number of obtained genes (91 vs 180), since the hypothesis

|beta| <= lfcThreshold

is more conservative.

ADD REPLYlink written 3.4 years ago by e.rempel980
1
gravatar for e.rempel
3.4 years ago by
e.rempel980
Germany, Heidelberg
e.rempel980 wrote:

Hi,

you could filter the res object directly in R

res <- res[res$log2FoldChange >= 1, ]
ADD COMMENTlink written 3.4 years ago by e.rempel980

Thanks for the reply. I tried this but I get an error below:

Error in NSBS(i, x, exact = exact, upperBoundIsStrict = !allow.append) : subscript contains NAs

ADD REPLYlink written 3.4 years ago by EpiExplorer90
0
gravatar for firatuyulur
3.4 years ago by
firatuyulur300
firatuyulur300 wrote:

There is a function, complete.cases(), in R that helps you to get rid of NAs.

mydata[complete.cases(mydata),]

will give you your dataframe with no NAs. Then you can apply any numeric filtering you like.

ADD COMMENTlink modified 3.4 years ago • written 3.4 years ago by firatuyulur300

Thank you for your reply. But this is not I want.

ADD REPLYlink written 3.4 years ago by EpiExplorer90

Explain how this is not what you're looking for - a previous answer had you showing an error owing to NAs, so how is removing NAs not the way forward?

ADD REPLYlink written 3.4 years ago by _r_am32k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1455 users visited in the last hour
_