Search
Question: report genes in intersection of Venn diagram with topTable?
0
gravatar for MBWatson
3.1 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")
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")
ADD COMMENTlink modified 3.1 years ago by Gordon Smyth34k • written 3.1 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.

ADD REPLYlink modified 3.1 years ago • written 3.1 years ago by Aaron Lun20k
2
gravatar for Gordon Smyth
3.1 years ago by
Gordon Smyth34k
Walter and Eliza Hall Institute of Medical Research, Melbourne, Australia
Gordon Smyth34k 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")
ADD COMMENTlink written 3.1 years ago by Gordon Smyth34k
0
gravatar for James W. MacDonald
3.1 years ago by
United States
James W. MacDonald46k 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).

ADD COMMENTlink written 3.1 years ago by James W. MacDonald46k

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 REPLYlink written 3.1 years ago by Aaron Lun20k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.2.0
Traffic: 258 users visited in the last hour