report genes in intersection of Venn diagram with topTable?
2
0
Entering edit mode
MBWatson • 0
@mbwatson-7475
Last seen 8.1 years ago
United States

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")
x <- read.delim("~mapped2SA.matrix", header=TRUE, row.names="Symbols") 
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) 
targets=read.delim("~targets.txt", header=TRUE)
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")
rnaseq venn diagram limma toptable edgeR • 3.9k views
ADD COMMENT
0
Entering edit mode

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

ADD REPLY
2
Entering edit mode
@gordon-smyth
Last seen 8 hours ago
WEHI, Melbourne, Australia

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")
ADD COMMENT
0
Entering edit mode
@james-w-macdonald-5106
Last seen 13 hours ago
United States

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).

ADD COMMENT
0
Entering edit mode

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.

ADD REPLY

Login before adding your answer.

Traffic: 883 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