Setting FDR at 0.05
2
0
Entering edit mode
4.4 years ago

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 • 1.9k views
ADD COMMENT
1
Entering edit mode

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.

ADD REPLY
0
Entering edit mode

THank you! ive modified the post! its clearer to see it now. Thank you

ADD REPLY
0
Entering edit mode

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.

ADD REPLY
1
Entering edit mode
4.4 years ago
zubenel ▴ 120

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)
ADD COMMENT
0
Entering edit mode

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
ADD REPLY
1
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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.

ADD REPLY
0
Entering edit mode

thank you so much!!!

ADD REPLY

Login before adding your answer.

Traffic: 2806 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6