In this video, I explain the basics of analyzing Affymetrix microarray datasets using R markdown and the limma package.
The example is for a single-channel Affy microarray chip to analyze gene expression differences between two conditions.
Here is the link to the annotation data: www.thermofisher.com/us/en/ho...
Here is the script:
title: "Arrays data"
library(tidyverse)
library(limma)
library(oligo)
pd = read.AnnotatedDataFrame("pdata.txt")
affyData = read.celfiles(rownames(pData(pd)))
eset = rma(affyData, target="core")
expres_eset = exprs(eset)
Type = as.factor(pd$type)
design = model.matrix(~0+Type)
colnames(design) = levels(Type)
fit = lmFit(expres_eset, design)
contrast.matrix = makeContrasts(contrasts = "exp-ctl", levels=design)
fit2 = eBayes(contrasts.fit(fit, contrast.matrix))
allstats = topTable(fit2, number=nrow(expres_eset))
annotation = read.table("MoGene-2_0-st-v1.na36.mm10.transcript.csv",
sep=",", header=TRUE,
row.names=1,as.is=TRUE)
simple = annotation %% filter(Status == "protein_coding")
filtered = expres_eset[as.character(simple$probeset_id),]
fit = lmFit(filtered, design)
fit2 = eBayes(contrasts.fit(fit, contrast.matrix))
allstats = topTable(fit2, number=nrow(filtered))
allstats = allstats %% mutate(probeset_id = as.integer(rownames(allstats)))
ids = simple %% select(probeset_id, GeneName, GeneDesc, Status)
allstats = allstats %% left_join(ids) %%
relocate(probeset_id, GeneName, GeneDesc, Status)
filtered = as.data.frame(filtered) %%
mutate(probeset_id = rownames(filtered)) %%
relocate(probeset_id)
22 май 2024