Question: edgeR Code syntax for without replicate
0
gravatar for Sushant Pawar
4.0 years ago by
Nashik
Sushant Pawar0 wrote:

Hi all ,

We want only the Gene Expression for the experiment without replicate, So i just want to confirm that is this code is wright or not . I am providing code below.

And what should be the dispersion in "exactTest" . Please help me in that.

 

==============================================

library(edgeR)
Read_counts=read.delim("Table.txt",row.names = "Gene")
counts=subset(Read_counts,Read_counts$Control >= 10 | Read_counts$Case >=10)
y=DGEList(counts=counts,group = 1:2)
et=exactTest(y,dispersion = 0.04)
et
 
fdr=p.adjust(et$table$PValue)
cpmres <- cpm(y)[rownames(et),]
summary(de <- decideTestsDGE(et))
write.csv(et$table, file="DEresults.csv")
write.csv(fdr, file="fdrcorrection.csv")
raw=data.frame(et$table)
DEResult_Counts=merge(counts,raw,by="row.names",row.names=FALSE)
write.table(as.data.frame(DEResult_Counts),"DEResult_with_Count.tsv",row.names = FALSE)
detags <- rownames(y)[as.logical(de)]
plotSmear(et, de.tags=detags)
pdf("Volcano_Plot.pdf")
plot(raw$logFC,-log(raw$PValue,10),
     xlim = c(-6,6),ylim = c(0,9),
     col=ifelse(abs(raw$logFC)<=2 | -log(raw$PValue,10)<1 ,"black","red"),
     pch=".",cex=2,xlab="logFoldChange",ylab="-Log10(Pvalue)",main="Volcano Plot")   
abline(v=c(-2,2),h=1,lty=2)

legend(3,3,"Significant Differentially \n Expressed Genes",lty=1, col=c('red'), bty='n',
       cex=0.5)
text(-5,8,"Down Regulated Genes",cex=0.5)
text(5,8,"Up Regulated Genes",cex=0.5)
dev.off()

==========================================================================

edger • 1.0k views
ADD COMMENTlink modified 4.0 years ago by Aaron Lun25k • written 4.0 years ago by Sushant Pawar0
Answer: edgeR Code syntax for without replicate
0
gravatar for Aaron Lun
4.0 years ago by
Aaron Lun25k
Cambridge, United Kingdom
Aaron Lun25k wrote:

You should probably run:

y <- calcNormFactors(y)

... before running exactTest. We can't really help you with the value of the dispersion, because you don't have any replicates from which to estimate it for your data set. For well-controlled experiments in the lab (e.g., with mice or cell lines), values below 0.05 are usually observed. For human patient data, this is often higher (0.1 to 0.5) and for single-cell data, it is higher still (>1). It really depends on what system you're working with, and how good your sample preparation is. That's the whole point of doing experiments with replicates.

Some other comments; your p.adjust call will use the default Holm correction to control the FWER, not the BH correction to control the FDR. You could fix this by setting method="BH", but it would be simpler to just use:

raw <- topTags(et, n=Inf, sort.by="none")$table

Incidentally, this means you can just cbind it with the counts to get DEResult_Counts, as the row order is the same.

ADD COMMENTlink modified 4.0 years ago • written 4.0 years ago by Aaron Lun25k

Dear Aaron,

Thank you for reply.

I have add the syntax that you have suggested ,

And "Output<- topTags(et, n=Inf, sort.by="none")$table" will be my final  table for further analysis

please take a look at syntax once again.

===============================================================

library(edgeR)
Read_counts=read.delim("Table.txt",row.names = "Gene")
counts=subset(Read_counts,Read_counts$Control >= 10 | Read_counts$Case >=10)
y=DGEList(counts=counts,group = 1:2)
y <- calcNormFactors(y)

et=exactTest(y,dispersion = 0.04)
et
Output<- topTags(et, n=Inf, sort.by="none")$table

fdr=p.adjust(et$table$PValue)
cpmres <- cpm(y)[rownames(et),]
summary(de <- decideTestsDGE(et))
write.csv(et$table, file="DEresults.csv")
write.csv(fdr, file="fdrcorrection.csv")
raw=data.frame(et$table)
DEResult_Counts=merge(counts,raw,by="row.names",row.names=FALSE)
write.table(as.data.frame(DEResult_Counts),"DEResult_with_Count.tsv",row.names = FALSE)
detags <- rownames(y)[as.logical(de)]
plotSmear(et, de.tags=detags)


pdf("Volcano_Plot.pdf")
plot(raw$logFC,-log(raw$PValue,10),
     xlim = c(-6,6),ylim = c(0,9),
     col=ifelse(abs(raw$logFC)<=2 | -log(raw$PValue,10)<1 ,"black","red"),
     pch=".",cex=2,xlab="logFoldChange",ylab="-Log10(Pvalue)",main="Volcano Plot")   
abline(v=c(-2,2),h=1,lty=2)

legend(3,3,"Significant Differentially \n Expressed Genes",lty=1, col=c('red'), bty='n',
       cex=0.5)
text(-5,8,"Down Regulated Genes",cex=0.5)
text(5,8,"Up Regulated Genes",cex=0.5)
dev.off()


====================================================

ADD REPLYlink written 4.0 years ago by Sushant Pawar0

You haven't fixed the problem with fdr=p.adjust(et$table$PValue).

ADD REPLYlink written 4.0 years ago by Aaron Lun25k

Hello Aaron ,

Thank you for your support.

i have done with the "fdr" command please once again go through it.

And suggest if any thing require to complete the programme.

===============================================

library(edgeR)
Read_counts=read.delim("Table.txt",row.names = "Gene")
counts=subset(Read_counts,Read_counts$Control >= 10 | Read_counts$Case >=10)
y=DGEList(counts=counts,group = 1:2)
y <- calcNormFactors(y)

et=exactTest(y,dispersion = 0.04)
et
fdr=p.adjust(et$table$PValue,method="BH")
Output <- topTags(et, n=Inf, sort.by="none")$table

cpmres <- cpm(y)[rownames(et),]
summary(de <- decideTestsDGE(et))
write.csv(et$table, file="DEresults.csv")
write.csv(fdr, file="fdrcorrection.csv")
raw=data.frame(Output)
DEResult_Counts=merge(counts,raw,by="row.names",row.names=FALSE)
write.table(as.data.frame(DEResult_Counts),"DEResult_with_Count_0.04.tsv",row.names = FALSE)
detags <- rownames(y)[as.logical(de)]
plotSmear(et, de.tags=detags)

FileforiPathway=data.frame(DEResult_Counts[c(-2,-3,-5,-6)])
write.table(as.data.frame(FileforiPathway),"FileforiPathway0.04.tsv",row.names = FALSE,sep="\t")


pdf("Volcano_Plot.pdf")
plot(raw$logFC,-log(raw$PValue,10),
     xlim = c(-6,6),ylim = c(0,9),
     col=ifelse(abs(raw$logFC)<=2 | -log(raw$PValue,10)<1 ,"black","red"),
     pch=".",cex=2,xlab="logFoldChange",ylab="-Log10(Pvalue)",main="Volcano Plot")   
abline(v=c(-2,2),h=1,lty=2)

legend(3,3,"Significant Differentially \n Expressed Genes",lty=1, col=c('red'), bty='n',
       cex=0.5)
text(-5,8,"Down Regulated Genes",cex=0.5)
text(5,8,"Up Regulated Genes",cex=0.5)
dev.off()

=============================================

ADD REPLYlink written 4.0 years ago by Sushant Pawar0
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 16.09
Traffic: 186 users visited in the last hour