0
$\begingroup$

I'm new to R, so be patient with me. 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 although (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 ? Thanks in advance,

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) 

$\endgroup$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.