Entering edit mode
Moritz Kebschull
▴
100
@moritz-kebschull-4339
Last seen 9.6 years ago
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]]
