Pairwise comparison for DEG analysis
0
0
Entering edit mode
Melissa • 0
@33edf2a1
Last seen 10 months ago
Malaysia

I want to perform pairwise comparisons between the different groups. The output (heatmap) did not show any visible demarcation between the groups. Is there something wrong with my code? The input is raw RNA-seq count.

KIRP <- factor(clin.info$subtype, levels=c("1a", "1b", "1c", "2a", "2b", "2c","2d")) 
batch <- clin.info$admin.batch_number
design <- model.matrix(~0+KIRP+batch) 
rownames(design) <- rownames(clin.info) 

mrna.dge <- DGEList(counts=mat, group=KIRP, remove.zeros=T) 
keep <- filterByExpr(mrna.dge, design, min.count=50) 
mrna.dge <- mrna.dge[keep,,keep.lib.sizes=FALSE] 
mrna.dge <- calcNormFactors(mrna.dge)

v <- voom(mrna.dge, design, plot=TRUE) 
fit <- lmFit(v, design)
contrasts.matrix <- makeContrasts("KIRP1a-KIRP1b", "KIRP1a-KIRP1c","KIRP1a-KIRP2a", "KIRP1a-KIRP2b", "KIRP1a-KIRP2c", "KIRP1a-KIRP2d",
                                  "KIRP1b-KIRP1c","KIRP1b-KIRP2a", "KIRP1b-KIRP2b", "KIRP1b-KIRP2c", "KIRP1b-KIRP2d",
                                  "KIRP1c-KIRP2a","KIRP1c-KIRP2b", "KIRP1c-KIRP2c", "KIRP1c-KIRP2d",
                                  "KIRP2a-KIRP2b","KIRP2a-KIRP2c", "KIRP2a-KIRP2d",
                                  "KIRP2b-KIRP2c", "KIRP2b-KIRP2d",
                                  "KIRP2c-KIRP2d", 
                                  levels=design) 
fit2 <- contrasts.fit(fit, contrasts.matrix) 
efit <- eBayes(fit2)

Top DEGs

top.genes.all <- topTable(efit, number=Inf, adjust.method="BH", sort.by="F", p.value=0.05)

Input:

> dput(mat[1:10,1:10])
structure(c(999, 1009.61, 4077.44, 9350, 3401, 7964, 29, 46, 
38222, 6805, 536, 1569.99, 2000.28, 5853, 4190, 7052, 90, 33, 
11212, 11684, 2220, 2451.35, 3603.86, 8557, 7820, 9805, 22, 11, 
9827, 8669, 659, 2012.63, 1181.37, 4454, 2527, 7454, 23, 33, 
915, 13502, 790, 1869.43, 1522.12, 5264, 2777, 4450, 5, 26, 32164, 
5425, 1344, 5643.93, 1274.56, 14406, 6390, 13868, 4, 68, 400, 
20151, 1145, 2608.79, 4077.32, 11937, 2831, 14999, 17, 31, 859, 
9118, 1880, 1323.13, 7910.62, 16025, 11039, 16554, 11, 208, 8274, 
4014, 90, 535.66, 1835.07, 1239, 208, 2941, 19, 226, 4453, 3789, 
2140, 2365.18, 2375.95, 7034, 3195, 8222, 69, 319, 21306, 4883
), dim = c(10L, 10L), dimnames = list(c("ACE", "ACSM3", "ACTA2", 
"ACTR2", "ADAM10", "ADAR", "ADRA2B", "ADRA2C", "AEBP1", "AIFM1"
), c("TCGA.2K.A9WE.01", "TCGA.2Z.A9J1.01", "TCGA.2Z.A9J3.01", 
"TCGA.2Z.A9J6.01", "TCGA.2Z.A9J7.01", "TCGA.2Z.A9J8.01", "TCGA.2Z.A9JI.01", 
"TCGA.2Z.A9JJ.01", "TCGA.2Z.A9JQ.01", "TCGA.4A.A93W.01")))

> dput(clin.info[,c("survival", "subtype")][1:10,])
structure(list(survival = c("lts", "lts", "lts", "lts", "lts", 
"lts", "lts", "lts", "lts", "lts"), subtype = c("2a", "1a", "1a", 
"1a", "1b", "2b", "1a", "2a", "1a", "1c")), row.names = c("TCGA.Y8.A8S1.01", 
"TCGA.Y8.A8RZ.01", "TCGA.Y8.A8RY.01", "TCGA.Y8.A897.01", "TCGA.Y8.A896.01", 
"TCGA.Y8.A895.01", "TCGA.Y8.A894.01", "TCGA.WN.A9G9.01", "TCGA.V9.A7HT.01", 
"TCGA.UZ.A9Q0.01"), class = "data.frame")
limma • 594 views
ADD COMMENT
0
Entering edit mode

Isn't this the same question you posted before: Pairwise differential expression using limma ? Why do you think there might be something wrong with your code?

Please don't delete questions after someone has answered them. That really discourages people from answering any new questions.

ADD REPLY
0
Entering edit mode

The top.genes.all return significant DEGs but the pattern doesn't show up clearly in my heatmap.

ADD REPLY
0
Entering edit mode

Sorry, we can't trouble-shoot problems with your heatmap. The DE analysis looks ok.

ADD REPLY
0
Entering edit mode
ADD REPLY

Login before adding your answer.

Traffic: 381 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6