Search
Question: How to extract lfc & p-values for probes identified with a Venn diagram
0
gravatar for Moritz Kebschull
6.6 years ago by
Moritz Kebschull100 wrote:
Dear list, I am trying to generate a list of probes that are differentially regulated after stimulation in cells with a gene knockdown vs. the wildtype controls. I generated two constrasts and looked at the overlap with a Venn diagram. It looks like this: diff_con(2288(3069)1010)diff_kd I am interested in what´s going on in response to the gene knockdown - likely the 1010 probes regulated exclusively by the knockdown. Or would you recommend to take a higher p-value cut-off to reduce the number of significant genes. I can extract the IDs of the up- or down-regulated probes, and possibly annotate them, but how do I attach the individual fold changes and stats for both contrasts? I would also love look into the GO categories in these groups - any idea how to do that? Many thanks in advance for your help! Moritz (Fellow, University of Bonn, Germany) Samples (all n=3, Illumina Human HT12v3) control cells + vehicle control cells + CPT (stimulation) knockdown cells + vehicle knockdown cells + CPT What I did: library(beadarray) dataFile = "Spag4_data.txt" BSData = readBeadSummaryData(dataFile = dataFile, dec=",", qcFile = NULL, skip = 0, qc.skip = 0) BSData.quantile = normaliseIllumina(BSData, method = "quantile", transform = "log2") BSData.genes = BSData.quantile[which(fData(BSData)$Status == "Gene"), ] expressed = apply(Detection(BSData.genes) < 0.05, 1, any) BSData.filt = BSData.genes[expressed,] library(limma) samples = c(rep("con_veh",3), rep("con_cpt",3), rep("kd_veh",3), rep("kd_cpt",3) ) samples = as.factor(samples) design=model.matrix(~0 + samples) colnames(design) = levels(samples) fit = lmFit(exprs(BSData.filt), design) cont.matrix=makeContrasts(diff_con=con_veh - con_cpt, diff_kd=kd_veh - kd_cpt, levels=design) fit=contrasts.fit(fit, cont.matrix) fit$genes=fData(BSData.filt) ebFit=eBayes(fit) library(illuminaHumanv3.db) illuminaHumanv3() ids = as.character(fData(BSData)[, 1]) chr = mget(ids, illuminaHumanv3CHR, ifnotfound = NA) chrloc = mget(ids, illuminaHumanv3CHRLOC, ifnotfound = NA) refseq = mget(ids, illuminaHumanv3REFSEQ, ifnotfound = NA) genename = mget(ids, illuminaHumanv3GENENAME, ifnotfound = NA) symbol = mget(ids, illuminaHumanv3SYMBOL, ifnotfound = NA) aidsnno = cbind(Ill_ID = as.character(ids), Chr = as.character(chr), Loc = as.character(chrloc), RefSeq = as.character(refseq), Name = as.character(genename), Symbol = as.character(symbol)) ebFit$genes = anno con=topTable(ebFit, coef=1, number=10000, p.value=0.05, lfc=1) kd=topTable(ebFit, coef=2, number=10000, p.value=0.05, lfc=1) write.table(con, file="diff_exp_con.txt", sep="\t", dec=",") write.table(kd, file="diff_exp_kd.txt", sep="\t", dec=",") results=decideTests(ebFit,method="separate",adjust.method="BH",p.value =0.05) a <- vennCounts(results) print(a) vennDiagram(a) kd_only_up <- which(results[,1] == 0 & results[,2] == 1) kd_only_down <- which(results[,1] == 0 & results[,2] == -1) ids_kd_only_up = as.character(names(kd_only_up)) ids_kd_only_down = as.character(names(kd_only_down)) [[alternative HTML version deleted]]
go
ADD COMMENTlink modified 6.6 years ago by René Dreos80 • written 6.6 years ago by Moritz Kebschull100
0
gravatar for René Dreos
6.6 years ago by
René Dreos80
René Dreos80 wrote:
Hi Moritz to extract the fold changes of your target genes you should do: > ebFit$coefficients[kd_only_up,] and to extract all the stats have a look at the command topTable (and then use the result as the above example). For the GO enrichment I normally use the topGO library best r On 15 November 2010 13:22, Moritz Kebschull <endothel@gmail.com> wrote: > Dear list, > > I am trying to generate a list of probes that are differentially regulated > after stimulation in cells with a gene knockdown vs. the wildtype controls. > > I generated two constrasts and looked at the overlap with a Venn diagram. > > It looks like this: diff_con(2288(3069)1010)diff_kd > I am interested in what´s going on in response to the gene knockdown - > likely the 1010 probes regulated exclusively by the knockdown. Or would you > recommend to take a higher p-value cut-off to reduce the number of > significant genes. > > I can extract the IDs of the up- or down-regulated probes, and possibly > annotate them, but how do I attach the individual fold changes and stats > for > both contrasts? > > I would also love look into the GO categories in these groups - any idea > how > to do that? > > Many thanks in advance for your help! > > Moritz (Fellow, University of Bonn, Germany) > > > Samples (all n=3, Illumina Human HT12v3) > control cells + vehicle > control cells + CPT (stimulation) > knockdown cells + vehicle > knockdown cells + CPT > > > What I did: > > library(beadarray) > dataFile = "Spag4_data.txt" > BSData = readBeadSummaryData(dataFile = dataFile, dec=",", qcFile = NULL, > skip = 0, qc.skip = 0) > BSData.quantile = normaliseIllumina(BSData, method = "quantile", transform > = > "log2") > BSData.genes = BSData.quantile[which(fData(BSData)$Status == "Gene"), ] > expressed = apply(Detection(BSData.genes) < 0.05, 1, any) > BSData.filt = BSData.genes[expressed,] > > library(limma) > samples = c(rep("con_veh",3), rep("con_cpt",3), rep("kd_veh",3), > rep("kd_cpt",3) ) > samples = as.factor(samples) > design=model.matrix(~0 + samples) > colnames(design) = levels(samples) > fit = lmFit(exprs(BSData.filt), design) > cont.matrix=makeContrasts(diff_con=con_veh - con_cpt, diff_kd=kd_veh - > kd_cpt, levels=design) > fit=contrasts.fit(fit, cont.matrix) > fit$genes=fData(BSData.filt) > ebFit=eBayes(fit) > > library(illuminaHumanv3.db) > illuminaHumanv3() > ids = as.character(fData(BSData)[, 1]) > chr = mget(ids, illuminaHumanv3CHR, ifnotfound = NA) > chrloc = mget(ids, illuminaHumanv3CHRLOC, ifnotfound = NA) > refseq = mget(ids, illuminaHumanv3REFSEQ, ifnotfound = NA) > genename = mget(ids, illuminaHumanv3GENENAME, ifnotfound = NA) > symbol = mget(ids, illuminaHumanv3SYMBOL, ifnotfound = NA) > aidsnno = cbind(Ill_ID = as.character(ids), Chr = as.character(chr), Loc = > as.character(chrloc), RefSeq = as.character(refseq), Name = > as.character(genename), Symbol = as.character(symbol)) > ebFit$genes = anno > > > con=topTable(ebFit, coef=1, number=10000, p.value=0.05, lfc=1) > kd=topTable(ebFit, coef=2, number=10000, p.value=0.05, lfc=1) > write.table(con, file="diff_exp_con.txt", sep="\t", dec=",") > write.table(kd, file="diff_exp_kd.txt", sep="\t", dec=",") > > > results=decideTests(ebFit,method="separate",adjust.method="BH",p.val ue=0.05) > a <- vennCounts(results) > print(a) > vennDiagram(a) > kd_only_up <- which(results[,1] == 0 & results[,2] == 1) > kd_only_down <- which(results[,1] == 0 & results[,2] == -1) > > ids_kd_only_up = as.character(names(kd_only_up)) > ids_kd_only_down = as.character(names(kd_only_down)) > > [[alternative HTML version deleted]] > > > _______________________________________________ > Bioconductor mailing list > Bioconductor@stat.math.ethz.ch > https://stat.ethz.ch/mailman/listinfo/bioconductor > Search the archives: > http://news.gmane.org/gmane.science.biology.informatics.conductor > [[alternative HTML version deleted]]
ADD COMMENTlink written 6.6 years ago by René Dreos80
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: 357 users visited in the last hour