Question: Setting FDR at 0.05
0
gravatar for ovariohisterectomia
8 weeks ago by
ovariohisterectomia10 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 • 188 views
ADD COMMENTlink modified 8 weeks ago by RamRS25k • written 8 weeks ago by ovariohisterectomia10
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 8 weeks ago • written 8 weeks ago by ATpoint29k

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

ADD REPLYlink written 8 weeks ago by ovariohisterectomia10

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 8 weeks ago by ATpoint29k
1
gravatar for zubenel
8 weeks ago by
zubenel30
zubenel30 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 8 weeks ago • written 8 weeks ago by zubenel30

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 8 weeks ago • written 8 weeks ago by ovariohisterectomia10
1

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

ADD REPLYlink written 8 weeks ago by ovariohisterectomia10

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 8 weeks ago • written 8 weeks ago by zubenel30

thank you so much!!!

ADD REPLYlink written 8 weeks ago by ovariohisterectomia10
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: 2360 users visited in the last hour