Question: Setting FDR at 0.05
0
gravatar for ovariohisterectomia
5 months ago by
ovariohisterectomia20 wrote:

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 • 278 views
ADD COMMENTlink modified 5 months ago by RamRS27k • written 5 months ago by ovariohisterectomia20
1

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 REPLYlink modified 5 months ago • written 5 months ago by ATpoint35k

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

ADD REPLYlink written 5 months ago by ovariohisterectomia20

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 REPLYlink written 5 months ago by ATpoint35k
1
gravatar for zubenel
5 months ago by
zubenel80
zubenel80 wrote:

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 COMMENTlink modified 5 months ago • written 5 months ago by zubenel80

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 REPLYlink modified 5 months ago • written 5 months ago by ovariohisterectomia20
1

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

ADD REPLYlink written 5 months ago by ovariohisterectomia20

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 REPLYlink modified 5 months ago • written 5 months ago by zubenel80

thank you so much!!!

ADD REPLYlink written 5 months ago by ovariohisterectomia20
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: 2002 users visited in the last hour