Question

Biological coefficient of variation(BCV) calculation in edgeR

0

Entering edit mode

n.naharfancy • 0

@nnaharfancy-20022

Last seen 3.5 years ago

Imperial College London

Hi

I doing DE analysis using edgeR. I have 12 paired samples from human blood cells. 6 treated with drug and 6 untreated. After doing DE analysis I get mostly down-regulated genes (53 down vs 1 up) which seems a bit abnormal due to the fact that I am getting mostly down-regulated genes. I checked the BCV it is 0.1 which is also seems a bit lower for human samples. Can it be the reason of the results that I get. How else can I check. The PCA analysis also shows the samples are not clustered by treatments. Any suggestions? my design matrix is (~Patients+Treatment). Thanks in advance.

Nurun

edger bcv edgeR • 1.3k views

ADD COMMENT • link updated 5.2 years ago by Aaron Lun ★ 28k • written 5.2 years ago by n.naharfancy • 0

Martin Morgan · Answer 1 · 2019-02-27

1

Entering edit mode

Aaron Lun ★ 28k

@alun

Last seen 4 hours ago

The city by the bay

You don't show any analysis code, so it's hard to know whether you made a mistake anywhere.

Generally speaking, 53 down vs 1 up is not a strong imbalance. If you had 5000 down and 1 up, that would be more concerning. But your current results do not seem inherently problematic.

A BCV of 0.1 for human samples is excellent; usually replicates from human samples are quite variable. So there's no (obvious) problem with power here. Maybe the data's trying to tell you that there just aren't many DE genes.

ADD COMMENT • link 5.2 years ago Aaron Lun ★ 28k

0

Entering edit mode

Hi Aaron,

Thank you so much for your reply. Please check the code below.

x<-read.delim("Complete_Count.txt")
XBD<-x[,c(1,4,5,6,9,13,12,17,14,19,18,23,22)]
colnames(XBD)<-c("GeneID", "Pre_1", "Post1", "Pre_2", "Post2", "Pre_3","Post3", "Pre_4","Post4", "Pre_5","Post5", "Pre_6", "Post6")
protein_coding<-read.delim("protein_coding_genelist.txt")
XBD_proteinCoding<-merge(XBD, protein_coding, by="GeneID", all.x=FALSE, sort=FALSE)
TreatmentXBD <- factor(substring(colnames(XBD_proteinCoding[,c(2:7, 12,13)]),1,3))
TreatmentXBD <- relevel(TreatmentXBD, ref="Pre")
PatientXBD <- factor(substring(colnames(XBD_proteinCoding[,c(2:7, 12, 13)]),5,5))
yXBD<-DGEList(counts =XBD_proteinCoding[, c(2:7, 12, 13)], genes = XBD_proteinCoding[,c(1, 15, 16)], group=TreatmentXBD)
keepXBD <- rowSums(cpm(yXBD)>1) >= 2
table(keepXBD)
yXBD<-yXBD[keepXBD, , keep.lib.sizes=FALSE]
yXBD<-calcNormFactors(yXBD)
plotMDS(yXBD, col=rep(1:2, each=1))
designXBD<-model.matrix(~PatientXBD+TreatmentXBD)
yXBD<-estimateDisp(yXBD,designXBD)
fit_glm_XBD <- glmFit(yXBD,designXBD)
lrtXBD <- glmLRT(fit_glm_XBD)

ADD REPLY • link updated 5.2 years ago by Martin Morgan 25k • written 5.2 years ago by n.naharfancy • 0