This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Setting FDR at 0.05

Hi people, thank you for helping me.

I would like to ask for help because i cannot change the FDR from defalult 0.1 , to 0.05.

Briefly, my experiment is Differential Expression between Liver and Lung, and i want ,in the last step of this script, filter by FDR 0.05, but when i print the summary, it says the FDR was 0.1.

########### pulmon vs higado ################ 

library(DESeq2) 
library("gplots") 
library(RColorBrewer) 
library(genefilter)


###set directory##

setwd("/home/isma/Alignment_50bp_TEST/Featurecounts/Counts/merged")


###countdata###

countData <- read.table("countData.txt", header = TRUE, dec = ".", sep
 = "\t", row.names = "GeneID") colData <- read.table("colData.txt", header = TRUE, sep = "\t") head(colData)


##Create matrix
dds <- DESeqDataSetFromMatrix(countData = countData,
                           colData = colData,
                             design = ~treatment)

### Filter low counts ### 
dds <- estimateSizeFactors(dds) idx <- rowSums( counts(dds, normalized=TRUE) >= 50 ) >= 12  dds <- dds[idx,]


#### RUN DEA #### 
dds <- DESeq(dds)
dds

#Extract normalized counts### 
normalized_counts <- counts(dds, normalized=TRUE) View(normalized_counts)
write.table(as.data.frame(normalized_counts), file = "Count_DEseqNormCount.txt", sep = "\t", dec = ".")

### Differential expression genes ### 
res_HIGADOvsPULMON <- results(dds,contrast = c("treatment", "HIGADO", "PULMON")) 
summary(res_HIGADOvsPULMON) 

### Filter significatives genes for FDR ### 
res_HIGADOvsPULMON_Sig <- res_HIGADOvsPULMON[which(res_HIGADOvsPULMON$padj < 0.05),] 
summary(res_HIGADOvsPULMON_Sig) 

#output from summary(res_HIGADOvsPULMON_Sig) 
Out of 46 with nonzero total read count adjusted p-value < 0.1 LFC > 0 (up)       : 35, 76%
> LFC < 0 (down)     : 11, 24% outliers [1]       : 0, 0% low counts [2]
> : 0, 0% (mean count < 64) [1] see 'cooksCutoff' argument of ?results
> [2] see 'independentFiltering' argument of ?results
rna-seq
• 3,323 views • link
updated by Ram Ram • written by ovariohisterectomia • ovariohisterectomia

For your information, I removed the script because it was very malformatted. You can simply copy and paste the script out of Rstudio and then highlight code with the code button 10101 to have it properly displayed. As from what I've seen in it you simply must do:

summary(res, alpha=0.05)

The alpha argument indicates your FDR choice.

• 0 views • link
written by ATpoint ATpoint

Looks better but please for the future focus on the relevant parts. All the lines above this line:

res_HIGADOvsPULMON <- results(dds,contrast = c("treatment", "HIGADO", "PULMON"))

are not relevant for the question as you are interested in the summary function.

• 0 views • link
written by ATpoint ATpoint

2 answers

DESeq2 summary function has two main arguments: object and alpha which means the adjusted p-value cutoff. To change cutoff you need to change alpha argument value. For example: 

summary(object = res, alpha = 0.05)

Or shorter version without writing argument names:

summary(res, 0.05)
• 0 views • link
written by zubenel zubenel

Thank you so much, ive already done that, but the results are the same. Im i doing something wrong before this line? Thank you

res_HIGADOvsPULMON_Sig <- res_HIGADOvsPULMON[which(res_HIGADOvsPULMON$padj < 0.05),] 
summary(res_HIGADOvsPULMON_Sig)

out of 46 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)       : 35, 76%
LFC < 0 (down)     : 11, 24%
outliers [1]       : 0, 0%
low counts [2]     : 0, 0%
(mean count < 64)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

summary(res_HIGADOvsPULMON_Sig, 0.05)

out of 46 with nonzero total read count
adjusted p-value < 0.05
LFC > 0 (up)       : 35, 76%
LFC < 0 (down)     : 11, 24%
outliers [1]       : 0, 0%
low counts [2]     : 0, 0%
(mean count < 64)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results
• 0 views • link
written by ovariohisterectomia • ovariohisterectomia

By applying filter which(res_HIGADOvsPULMON$padj < 0.05) you have removed all the cases where adjusted p-values >= 0.05 and that includes all the cases where p-value is between 0.05 and 0.1.

If you want to compare summary results you have to use unfiltered res_HIGADOvsPULMON object.

• 0 views • link
written by zubenel zubenel

I performed the script once again (without changes) and now it works! thank you all

• 0 views • link
written by ovariohisterectomia • ovariohisterectomia
<h6>##### pulmon vs higado ################</h6> 
library(DESeq2) 
library("gplots") 
library(RColorBrewer) 
library(genefilter)

set directory

   setwd("/home/isma/Alignment_50bp_TEST/Featurecounts/Counts/merged")

countdata

countData <- read.table("countData.txt", header = TRUE, dec = ".", sep
 = "\t", row.names = "GeneID") colData <- read.table("colData.txt", header = TRUE, sep = "\t") head(colData)

Create matrix

dds <- DESeqDataSetFromMatrix(countData = countData,
                           colData = colData,
                             design = ~treatment)

Filter low counts ###

dds <- estimateSizeFactors(dds) idx <- rowSums( counts(dds, normalized=TRUE) >= 50 ) >= 12  dds <- dds[idx,]

RUN DEA ####

dds <- DESeq(dds)
 dds

Extract normalized counts###

normalized_counts <- counts(dds, normalized=TRUE) View(normalized_counts)
write.table(as.data.frame(normalized_counts), file = "Count_DEseqNormCount.txt", sep = "\t", dec = ".")

Differential expression genes ###

res_HIGADOvsPULMON <- results(dds,contrast = c("treatment", "HIGADO", "PULMON")) 
summary(res_HIGADOvsPULMON)

Filter significatives genes for FDR ###

res_HIGADOvsPULMON_Sig <- res_HIGADOvsPULMON[which(res_HIGADOvsPULMON$padj < 0.05),] 
summary(res_HIGADOvsPULMON_Sig)

output from summary(res_HIGADOvsPULMON_Sig)

Out of 46 with nonzero total read count adjusted p-value < 0.1 LFC > 0 (up)       : 35, 76%
> LFC < 0 (down)     : 11, 24% outliers [1]       : 0, 0% low counts [2]
> : 0, 0% (mean count < 64) [1] see 'cooksCutoff' argument of ?results
> [2] see 'independentFiltering' argument of ?results
• 0 views • link
written by ovariohisterectomia • ovariohisterectomia

Log in to answer this question.

More posts like this