Search
Question: report genes in intersection of Venn diagram with topTable?
0
3.3 years ago by
MBWatson0
United States
MBWatson0 wrote:

Good afternoon,

I am analyzing RNA-Seq data generated from an experiment to test the effects of salinity (15ppt, 35ppt, 60ppt) on different populations (BR and SD). My code for reporting the genes differentially expressed in each contrast is working well but I am wondering if topTable can report the list of genes that occur in the intersection of the two circles in the Venn diagram. My code is below. Thank you in advance for any advice.

source("http://bioconductor.org/biocLite.R")
biocLite("edgeR")
library("edgeR")
biocLite("limma")
biocLite("statmod")
z <-x[ ,c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18)]
keep <- rowSums(cpm(z)>1) >=3
z<- z[keep,]
dge <- DGEList(counts=z)
dgeN <-calcNormFactors(dge)
attach(targets)
pop.sal <-paste(targets$pop, targets$sal, sep=".")
pop.sal <-factor(pop.sal, level=c("br.15", "br.35", "br.60", "sd.15", "sd.35", "sd.60"))
design <-model.matrix(~0+pop.sal)
colnames(design) <- levels(pop.sal)
v <- voom(dgeN,design,plot=TRUE)
fit <- lmFit(v, design)
cont.matrix <- makeContrasts(br.low=br.15-br.35, sd.low=sd.15-sd.35, levels=design)
fit2 <- contrasts.fit(fit, cont.matrix1)
fit2 <- eBayes(fit2)
results <- decideTests(fit2)
vennDiagram(results)
BRlow <- topTable(fit2, 1, adjust.method="BH", number=30000, p.value=0.05, sort.by="p")
write.csv(BRlow, file="BR_SA.csv")
SDlow <- topTable(fit2, 2, adjust.method="BH", number=30000, p.value=0.05, sort.by="p")
write.csv(SDlow, file="SDlow_SA.csv")
modified 3.2 years ago by Gordon Smyth35k • written 3.3 years ago by MBWatson0

Set n=Inf if you want to get all genes from topTable. It's less typing and it's more obvious to the reader.

2
3.2 years ago by
Gordon Smyth35k
Walter and Eliza Hall Institute of Medical Research, Melbourne, Australia
Gordon Smyth35k wrote:

Yes, it's easy. As James and Aaron have said, this is what decideTests() is for:

intersection <- results[,"br.low"] & results[,"sd.low"]
topTable(fit2[intersection,], n=Inf, sort.by="p")
0
3.3 years ago by
United States
James W. MacDonald47k wrote:

The topTable function isn't designed to know anything about intersections. That is what decideTests is for. I don't think there is anything in limma to report the list of genes in an intersection, so I have code in affycoretools that I use to do that sort of thing.

I have everything primarily wrapped up in two high level functions, makeVenn() and plotVenn(). The former uses ReportingTools to make HTML pages (based on topTable output) for each of the 'cells' of the venn diagram, and plotVenn() makes an HTML page with a Venn diagram where the numbers in each cell are clickable links to the HTML pages containing those genes.

As a side effect, makeVenn() also outputs csv files for each cell of the Venn diagram as well, although given the apparent awesomeness of the openxlsx package I may end up converting to outputting Excel workbooks, as that is what most people do with the csv files anyway (and directly writing to xlsx files eliminates Excel's nasty habit of converting date looking things to actual date format).

To follow up on decideTests, you should be able to do:

intersected <- rowSums(results[, c("br.low", "sd.low")] != 0L) == 2L

to identify those genes that are DE in both of your contrasts.