I am performing differential gene expression analysis to microarray data for type 2 diabetes donors and nondiabetic donors. When I run the code I get some different results in each time (about 50 or so different genes in each run). Is this normal? I am submitting a project and it could be very problematic if my results change in each run, since I will perform downstream analysis for the DE genes and use web tools such as DAVID for further analysis.
Is there a way to keep the results from changing in each run?
g <- getGEO("GSE20966",GSEMatrix = TRUE)
g<- g[[1]]
exprs(g)
exprs(g) <- log2(exprs(g))
#Diffrential gene exp
library(limma)
design <- model.matrix(~0+sampleInfo$characteristics_ch1.1)
colnames(design) <- c("NonDiabetic","T2D")
summary(exprs(g))
exp <- na.omit(exprs(g))
cutoff <- median(exp)
## TRUE or FALSE for whether each gene is "expressed" in each sample
is_expressed <- exprs(g) > cutoff
## Identify genes expressed in more than 2 samples
keep <- rowSums(is_expressed) > 2
## check how many genes are removed / retained.
table(keep)
keep<- na.omit(keep)
## subset to just those expressed genes
g<- g[keep,]
fit <- lmFit(exp, design)
head(fit$coefficients)
contrasts <- makeContrasts(T2D - NonDiabetic, levels=design)
fit2 <- contrasts.fit(fit, contrasts)
fit2 <- eBayes(fit2)
topTable(fit2)
topTable(fit2, coef=1)
decideTests(fit2)
table(decideTests(fit2))
aw <- arrayWeights(exp,design)
aw
fit <- lmFit(exp, design,
weights = aw)
contrasts <- makeContrasts(T2D - NonDiabetic, levels=design)
fit2 <- contrasts.fit(fit, contrasts)
fit2 <- eBayes(fit2)
g <- na.omit(g)
anno <- fData(g)
anno
anno <- select(anno,`Gene Symbol`,ENTREZ_GENE_ID)
fit2$genes <- anno
topTable(fit2)
full_results <- topTable(fit2, number=Inf)
full_results <- tibble::rownames_to_column(full_results,"ID")
library(ggplot2)
dev.off()
ggplot(full_results,aes(x = logFC, y=B)) + geom_point()
p_cutoff <- 0.05
fc_cutoff <- 1
full_results %>%
mutate(Significant = adj.P.Val < p_cutoff, abs(logFC) > fc_cutoff )
%>%
mutate(Rank = 1:n(), Label = ifelse(Rank < topN, anno$`Gene
Symbol`,"")) %>%
ggplot(aes(x = logFC, y = B, col=Significant,label=Label)) +
geom_point() + geom_text_repel(col="black")
p_cutoff <- 0.05
fc_cutoff <- 1
filter(full_results, adj.P.Val < 0.05, abs(logFC) > 1)
library(pheatmap)
topN <- 20
top<-150
gene_names <- mutate(full_results, Rank = 1:n()) %>%
filter(Rank < top) %>%
pull(Gene.Symbol)