Question: edgeR Code syntax for without replicate
0
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.

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

library(edgeR)
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
modified 4.0 years ago by Aaron Lun25k • written 4.0 years ago by Sushant Pawar0
Answer: edgeR Code syntax for without replicate
0
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)
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()

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