In this video, I show you how to plot microarray data using R. I cover the
basics on how to get the data from .CEL files, and then I show you how to
generate volcano plots and expression plots.
0:00 Introduction
0:13 Data Preparation
3:05 Volcano plots
4:05 Expression plots
5:35 Outro
library(tidyverse)
library(limma)
library(oligo)
library(ggplot2)
pd = read.AnnotatedDataFrame("pdata.txt")
affyData = read.celfiles(rownames(pData(pd)))
eset = rma(affyData, target="core")
expres_eset = exprs(eset)
annotation = read.table("MoGene-2_0-st-v1.na36.mm10.transcript.csv",
sep=",", header=TRUE,
row.names=1,as.is=TRUE)
annotation = annotation %% rowwise() %%
mutate(GeneName = trimws(unlist(strsplit(gene_assignment,"//"))[2])) %%
mutate(GeneDesc = unlist(strsplit(gene_assignment,"//"))[3]) %%
mutate(Status = unlist(strsplit(unlist(strsplit(mrna_assignment,"gene_biotype:"))[2]," "))[1]) %%
select(probeset_id, GeneName, GeneDesc, Status)
Type = as.factor(pd$type)
design = model.matrix(~0+Type)
colnames(design) = levels(Type)
fit = lmFit(eset, design)
contrast.matrix = makeContrasts(contrasts = "exp-ctl", levels=design)
fit2 = eBayes(contrasts.fit(fit, contrast.matrix))
allstats = topTable(fit2, number=nrow(eset))
allstats = allstats %% mutate(probeset_id = as.integer(rownames(allstats))) %%
left_join(annotation) %%
relocate(probeset_id, GeneName, GeneDesc, Status) %%
mutate(Significant = ifelse(adj.P.Val x 0.05 & abs(logFC) 1, "Yes", "No"))
ggplot(allstats, aes(x=logFC, y=-log10(adj.P.Val), color=Significant)) +
geom_point() +
theme_minimal() +
labs(x="logFC", y="-log10(adj.P.Val)")
pdd = read.table("pdata.txt", header=TRUE, sep="\t")
plot_table = as.data.frame(expres_eset) %%
mutate(probeset_id = as.integer(rownames(expres_eset))) %%
relocate(probeset_id) %%
pivot_longer(-probeset_id, names_to="filename", values_to="Expression") %%
left_join(pdd) %%
group_by(probeset_id, type) %%
summarise(mean = mean(Expression)) %%
left_join(allstats %%
select(probeset_id, GeneName, logFC, adj.P.Val, Significant)) %%
pivot_wider(names_from=type, values_from=mean) %%
filter(!is.na(GeneName))
ggplot(plot_table, aes(x=ctl, y=exp, color=Significant)) +
geom_point() +
theme_minimal() +
labs(x="Control", y="Experimental")
11 июл 2024
