GSE118337
29 April 2026
Abstract
This report analyses the GSE118337 microarray dataset using the limma package.
Contents
1 Read Target Sample Information
## Accession ID Condition FileName Replicate Cell
## GSM3325550 GSM3325550 465 HK2_ctrl GSM3325550_465.txt.gz 1 HK2
## GSM3325551 GSM3325551 466 HK2_TGFB GSM3325551_466.txt.gz 1 HK2
## GSM3325552 GSM3325552 467 HK2_Empa GSM3325552_467.txt.gz 1 HK2
## GSM3325553 GSM3325553 468 HK2_Cana GSM3325553_468.txt.gz 1 HK2
## GSM3325554 GSM3325554 469 RPTEC_ctrl GSM3325554_469.txt.gz 1 RPTEC
## GSM3325555 GSM3325555 470 RPTEC_TGFB GSM3325555_470.txt.gz 1 RPTEC
## GSM3325556 GSM3325556 471 RPTEC_Emp GSM3325556_471.txt.gz 1 RPTEC
## GSM3325557 GSM3325557 472 RPTEC_Cana GSM3325557_472.txt.gz 1 RPTEC
## GSM3325558 GSM3325558 473 HK2_ctrl2 GSM3325558_473.txt.gz 2 HK2
## GSM3325559 GSM3325559 474 HK2_TGFB2 GSM3325559_474.txt.gz 2 HK2
## GSM3325560 GSM3325560 475 HK2_Empa2 GSM3325560_475.txt.gz 2 HK2
## GSM3325561 GSM3325561 476 HK2_Cana2 GSM3325561_476.txt.gz 2 HK2
## GSM3325562 GSM3325562 477 RPTEC_ctrl2 GSM3325562_477.txt.gz 2 RPTEC
## GSM3325563 GSM3325563 478 RPTEC_TGFB2 GSM3325563_478.txt.gz 2 RPTEC
## GSM3325564 GSM3325564 479 RPTEC_Empa2 GSM3325564_479.txt.gz 2 RPTEC
## GSM3325565 GSM3325565 480 RPTEC_Cana2 GSM3325565_480.txt.gz 2 RPTEC
## Treatment
## GSM3325550 ctrl
## GSM3325551 TGFB
## GSM3325552 Empa
## GSM3325553 Cana
## GSM3325554 ctrl
## GSM3325555 TGFB
## GSM3325556 Empa
## GSM3325557 Cana
## GSM3325558 ctrl
## GSM3325559 TGFB
## GSM3325560 Empa
## GSM3325561 Cana
## GSM3325562 ctrl
## GSM3325563 TGFB
## GSM3325564 Empa
## GSM3325565 Cana2 Read Probe Intensities
convert the data to an EListRaw object, which is a data object for single channel data, specify green.only = TRUE for Agilent ,retain information about background via gIsWellAboveBG.
project <- read.maimages(
targets,
source = 'agilent',
green.only = TRUE,
other.columns = 'gIsWellAboveBG',
path = "data/GSE118337/")## Read data/GSE118337//GSM3325550_465.txt.gz
## Read data/GSE118337//GSM3325551_466.txt.gz
## Read data/GSE118337//GSM3325552_467.txt.gz
## Read data/GSE118337//GSM3325553_468.txt.gz
## Read data/GSE118337//GSM3325554_469.txt.gz
## Read data/GSE118337//GSM3325555_470.txt.gz
## Read data/GSE118337//GSM3325556_471.txt.gz
## Read data/GSE118337//GSM3325557_472.txt.gz
## Read data/GSE118337//GSM3325558_473.txt.gz
## Read data/GSE118337//GSM3325559_474.txt.gz
## Read data/GSE118337//GSM3325560_475.txt.gz
## Read data/GSE118337//GSM3325561_476.txt.gz
## Read data/GSE118337//GSM3325562_477.txt.gz
## Read data/GSE118337//GSM3325563_478.txt.gz
## Read data/GSE118337//GSM3325564_479.txt.gz
## Read data/GSE118337//GSM3325565_480.txt.gzcolnames(project) <- targets$Condition
dim(project)## [1] 44495 16length(unique(project$genes$ProbeName)) ## [1] 34183in the article they refer these same number of probes as “genes”
3 Background correct and normalize
# Correction
project.bgcorrect = backgroundCorrect(project, method="normexp")## Array 1 corrected
## Array 2 corrected
## Array 3 corrected
## Array 4 corrected
## Array 5 corrected
## Array 6 corrected
## Array 7 corrected
## Array 8 corrected
## Array 9 corrected
## Array 10 corrected
## Array 11 corrected
## Array 12 corrected
## Array 13 corrected
## Array 14 corrected
## Array 15 corrected
## Array 16 corrected# Then normalize and log-transformed the data.
project.bgcorrect.norm = normalizeBetweenArrays(project.bgcorrect,
method="quantile")
# Before normalization.
boxplot(log(as.matrix(project.bgcorrect)),las=2,ylab="Log2(Intensity)")
# After it.
boxplot(as.matrix(project.bgcorrect.norm), las=2, ylab = "Log2(Intensity)")
4 Gene Filtering
# filter out control probes
Control <- project.bgcorrect.norm$genes$ControlType==1L
# Keep probes that are above background on at least two arrays
IsExpr <- rowSums(project.bgcorrect.norm$other$gIsWellAboveBG > 0) >= 2
project.bgcorrect.norm.filt <- project.bgcorrect.norm[!Control & IsExpr, ]
dim(project.bgcorrect.norm)## [1] 44495 16dim(project.bgcorrect.norm.filt)## [1] 28685 16# remove annotation columns we no longer need
project.bgcorrect.norm.filt$genes <- project.bgcorrect.norm.filt$genes[
,c("ProbeName", "GeneName", "SystematicName")]
head(project.bgcorrect.norm.filt$genes)## ProbeName GeneName SystematicName
## 13 A_23_P42935 BRAF NM_004333
## 14 A_23_P117082 HEBP1 NM_015987
## 15 A_23_P2683 RPAP3 NM_024604
## 16 A_24_P358131 ENST00000404956 ENST00000404956
## 19 A_32_P14850 NPIPB15 ENST00000429990
## 20 A_23_P158596 AGTRAP NM_001040196#Finally calculate the average intensity values from the probes of each gene.
# project.bgcorrect.norm.mean = avereps(project.bgcorrect.norm.filt,
# ID=project.bgcorrect.norm.filt$genes$ProbeName)
#
# dim(project.bgcorrect.norm.mean)5 Data Exploration
plotMDS(project.bgcorrect.norm.filt)
# H_TGFB1 seems like a big outlier so we will remove it
project.bgcorrect.norm.filt <- project.bgcorrect.norm.filt[,-2]
plotMDS(project.bgcorrect.norm.filt)
# RPTEC samples all cluster together, while HK2 separate6 Differential Expression
6.1 Using Group
Replicate <- as.factor(project.bgcorrect.norm.filt$targets$Replicate)
Group <- paste(project.bgcorrect.norm.filt$targets$Cell, sep=".",
project.bgcorrect.norm.filt$targets$Treatment)
Group <- as.factor(Group)
design <- model.matrix(~ 0 + Group + Replicate)
fit <- lmFit(project.bgcorrect.norm.filt, design)
# cotrasts of drug should be against TGFB but i took out at HK2
contrast_matrix <- makeContrasts(
HK2_CANA_vs_CONTROL = GroupHK2.Cana - GroupHK2.ctrl,
HK2_EMPA_vs_CONTROL = GroupHK2.Empa - GroupHK2.ctrl,
HK2_TGFB_vs_CONTROL = GroupHK2.TGFB - GroupHK2.ctrl,
HK2_CANA_vs_TGFB = GroupHK2.Cana - GroupHK2.TGFB,
HK2_EMPA_vs_TGFB = GroupHK2.Empa - GroupHK2.TGFB,
HK2_CANA_vs_EMPA = GroupHK2.Cana - GroupHK2.Empa,
RPTEC_CANA_vs_CONTROL = GroupRPTEC.Cana - GroupRPTEC.ctrl,
RPTEC_EMPA_vs_CONTROL = GroupRPTEC.Empa - GroupRPTEC.ctrl,
RPTEC_TGFB_vs_CONTROL = GroupRPTEC.TGFB - GroupRPTEC.ctrl,
RPTEC_CANA_vs_TGFB = GroupRPTEC.Cana - GroupRPTEC.TGFB,
RPTEC_EMPA_vs_TGFB = GroupRPTEC.Empa - GroupRPTEC.TGFB,
RPTEC_CANA_vs_EMPA = GroupRPTEC.Cana - GroupRPTEC.Empa,
levels = design
)
fit2 <- contrasts.fit(fit, contrast_matrix)
fit2 <- eBayes(fit2, trend = TRUE, robust = TRUE)
plotSA(fit2)
summary(decideTests(fit2))## HK2_CANA_vs_CONTROL HK2_EMPA_vs_CONTROL HK2_TGFB_vs_CONTROL
## Down 7 5 318
## NotSig 28673 28679 27849
## Up 5 1 518
## HK2_CANA_vs_TGFB HK2_EMPA_vs_TGFB HK2_CANA_vs_EMPA RPTEC_CANA_vs_CONTROL
## Down 710 706 0 3
## NotSig 27612 27566 28685 28682
## Up 363 413 0 0
## RPTEC_EMPA_vs_CONTROL RPTEC_TGFB_vs_CONTROL RPTEC_CANA_vs_TGFB
## Down 12 0 279
## NotSig 28672 28684 28364
## Up 1 1 42
## RPTEC_EMPA_vs_TGFB RPTEC_CANA_vs_EMPA
## Down 301 0
## NotSig 28349 28685
## Up 35 06.2 Separate for each Cell
6.2.1 HK2
y.HK2 <- project.bgcorrect.norm.filt[,project.bgcorrect.norm.filt$targets$Cell=="HK2"]
plotMDS(y.HK2, main="HK2", label=y.HK2$targets$Condition)
Replicate <- as.factor(y.HK2$targets$Replicate)
Treatment <- factor(y.HK2$targets$Treatment, levels = c("ctrl", "TGFB", "Empa", "Cana"))
design <- model.matrix(~ 0 + Treatment + Replicate)
fit.HK2 <- lmFit(y.HK2, design)
contrast_matrix <- makeContrasts(
HK2_CANA_vs_CONTROL = TreatmentCana - Treatmentctrl,
HK2_EMPA_vs_CONTROL = TreatmentEmpa - Treatmentctrl,
HK2_TGFB_vs_CONTROL = TreatmentTGFB - Treatmentctrl,
HK2_CANA_vs_TGFB = TreatmentCana - TreatmentTGFB,
HK2_EMPA_vs_TGFB = TreatmentEmpa - TreatmentTGFB,
HK2_CANA_vs_EMPA = TreatmentCana - TreatmentEmpa,
levels = design)
fit.HK2 <- contrasts.fit(fit.HK2, contrast_matrix)
fit.HK2 <- eBayes(fit.HK2, trend = TRUE, robust = TRUE)
plotSA(fit.HK2)
summary(decideTests(fit.HK2))## HK2_CANA_vs_CONTROL HK2_EMPA_vs_CONTROL HK2_TGFB_vs_CONTROL
## Down 626 597 955
## NotSig 27537 27565 26746
## Up 522 523 984
## HK2_CANA_vs_TGFB HK2_EMPA_vs_TGFB HK2_CANA_vs_EMPA
## Down 1416 1391 0
## NotSig 25991 25969 28685
## Up 1278 1325 06.2.2 RPTEC
y.RPTEC <- project.bgcorrect.norm.filt[,project.bgcorrect.norm.filt$targets$Cell=="RPTEC"]
plotMDS(y.RPTEC, main="RPTEC", label=y.RPTEC$targets$Condition)
Replicate <- as.factor(y.RPTEC$targets$Replicate)
Treatment <- factor(y.RPTEC$targets$Treatment, levels = c("ctrl", "TGFB", "Empa", "Cana"))
design <- model.matrix(~ 0 + Treatment + Replicate)
fit.RPTEC <- lmFit(y.RPTEC, design)
contrast_matrix <- makeContrasts(
RPTEC_CANA_vs_CONTROL = TreatmentCana - Treatmentctrl,
RPTEC_EMPA_vs_CONTROL = TreatmentEmpa - Treatmentctrl,
RPTEC_TGFB_vs_CONTROL = TreatmentTGFB - Treatmentctrl,
RPTEC_CANA_vs_TGFB = TreatmentCana - TreatmentTGFB,
RPTEC_EMPA_vs_TGFB = TreatmentEmpa - TreatmentTGFB,
RPTEC_CANA_vs_EMPA = TreatmentCana - TreatmentEmpa,
levels = design)
fit.RPTEC.2 <- contrasts.fit(fit.RPTEC, contrast_matrix)
fit.RPTEC.2 <- eBayes(fit.RPTEC.2, trend = TRUE, robust = TRUE)
plotSA(fit.RPTEC.2)
summary(decideTests(fit.RPTEC.2))## RPTEC_CANA_vs_CONTROL RPTEC_EMPA_vs_CONTROL RPTEC_TGFB_vs_CONTROL
## Down 561 662 58
## NotSig 27719 27585 28290
## Up 405 438 337
## RPTEC_CANA_vs_TGFB RPTEC_EMPA_vs_TGFB RPTEC_CANA_vs_EMPA
## Down 1823 1598 0
## NotSig 25067 25941 28685
## Up 1795 1146 0R.cana.vs.TGFB.res <- topTable(fit.RPTEC.2, coef = "RPTEC_CANA_vs_TGFB", number = Inf)
sum(duplicated(R.cana.vs.TGFB.res$ProbeName))## [1] 6888R.cana.vs.TGFB.res %>% filter(GeneName == "NegativeControl")## ProbeName GeneName SystematicName logFC AveExpr t
## 35476 (-)3xSLv1 NegativeControl NegativeControl 2.1971974 4.097197 4.8019783
## 37333 (-)3xSLv1 NegativeControl NegativeControl 1.3851634 4.239748 4.1008296
## 12706 (-)3xSLv1 NegativeControl NegativeControl -0.7668852 4.295831 -1.1219886
## 36054 (-)3xSLv1 NegativeControl NegativeControl 0.8937647 4.341827 0.9190497
## 11076 (-)3xSLv1 NegativeControl NegativeControl 0.1976275 4.154849 0.2211915
## P.Value adj.P.Val B
## 35476 0.001328503 0.02031349 -0.9283979
## 37333 0.002713618 0.03054712 -1.7535971
## 12706 0.305172665 0.52900711 -6.2642226
## 36054 0.397417188 0.61501468 -6.4637796
## 11076 0.832928044 0.91808465 -6.8957330