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.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]]
>
>
> _______________________________________________
> 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]]

