1803IMA
Brandt Pence
October 9, 2018
Instructions
In order to use the .Rmd file, the following datasets are required: gst.csv, gst_counts.csv, infl_data.csv, infl_counts.csv, lps_data.csv, lps_counts.csv, qpcr_data_infl.csv, qpcr_data_lps.csv, and demographics.csv. Add the .Rmd file to the same folder as the datasets to run. Otherwise, use setwd(). Required libraries are dplyr, plotrix, DescTools, and car.
Custom Functions
The Seahorse analysis for this project relies on a group of custom functions for data normalization, replicate grouping, glycolysis stress test calculations, and graphing. These are included below, but are not shown due to space constraints. Download the R Markdown (.Rmd) file for the code.
Options and libraries
options(scipen=999)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(car)##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recodelibrary(plotrix)
library(DescTools)##
## Attaching package: 'DescTools'## The following object is masked from 'package:car':
##
## RecodeDemographics
The following data are demographic and anthropometric.
screen <- read.csv("demographics.csv", header=T)Cohort 1
Subset data for cohort 1 subjects only
cohort1 <- subset(screen, screen[,"cohort"]==1)Chi-square analysis for categorical variables:
Sex:
tblsex <- table(cohort1$sex, cohort1$group)
chisq.test(tblsex, correct=F)## Warning in chisq.test(tblsex, correct = F): Chi-squared approximation may
## be incorrect##
## Pearson's Chi-squared test
##
## data: tblsex
## X-squared = 0.22222, df = 1, p-value = 0.6374Race:
tblrace <- table(cohort1$group, cohort1$race)
chisq.test(tblrace, correct=F)## Warning in chisq.test(tblrace, correct = F): Chi-squared approximation may
## be incorrect##
## Pearson's Chi-squared test
##
## data: tblrace
## X-squared = NaN, df = 3, p-value = NASubsets for aged and young subjects
c1aged <- subset(cohort1, group=="A")
c1young <- subset(cohort1, group=="Y")Exploratory boxplots for continuous variables
par(mfrow=c(2,2))
for (i in 4:7) {
y <- colnames(cohort1)[i]
variab <- cohort1[,i]
plot(variab~group, data=cohort1, ylab=y, xlab=NULL)
stripchart(variab~group, data=cohort1, vertical=T, method="jitter", add=T, pch=20, col="red")
}
Tests for equality of variances. All variances are equal by Levene’s test.
eqvartest <- apply(cohort1[4:7], 2, function(x) {leveneTest(x,cohort1$group)$'Pr(>F)'})[1,]
eqvartest<0.05## age height weight bmi
## FALSE FALSE FALSE FALSET-tests with output of p-values
c1ttests <- apply(cohort1[c(4:7)], 2, function(x) {t.test(x~group, data=cohort1, var.equal=T)$p.value})
c1ttests2 <- apply(cohort1[c(4:7)], 2, function(x) {t.test(x~group, data=cohort1, var.equal=T)})
round(c1ttests, 4)## age height weight bmi
## 0.0000 0.1663 0.8137 0.6881Calculations for means and SEs
c1meansyoung <- apply(c1young[c(4:7)], 2, function(x) {mean(x, na.rm=T)})
c1sesyoung <- apply(c1young[c(4:7)], 2, function(x) {std.error(x, na.rm=T)})
c1meansaged <- apply(c1aged[c(4:7)], 2, function(x) {mean(x, na.rm=T)})
c1sesaged <- apply(c1aged[c(4:7)], 2, function(x) {std.error(x, na.rm=T)})Means for continuous variables
c1means <- rbind(c1meansyoung, c1meansaged)
row.names(c1means) <- c("young","aged")
c1means## age height weight bmi
## young 25.77778 169.2222 77.31111 26.43333
## aged 65.00000 176.6667 79.55556 25.51111SEMs for continuous variables
c1ses <- rbind(c1sesyoung, c1sesaged)
row.names(c1ses) <- c("young", "aged")
c1ses## age height weight bmi
## young 1.913242 4.162591 8.648594 1.949430
## aged 1.166667 3.004626 3.603398 1.135714Cohort 2
Subset data for cohort 2 subjects only
cohort2 <- subset(screen, screen[,"cohort"]==2)Sex:
c2tblsex <- table(cohort2$sex, cohort2$group)
chisq.test(c2tblsex, correct=F)## Warning in chisq.test(c2tblsex, correct = F): Chi-squared approximation may
## be incorrect##
## Pearson's Chi-squared test
##
## data: c2tblsex
## X-squared = 0.13468, df = 1, p-value = 0.7136Race:
c2tblrace <- table(cohort2$group, cohort2$race)
c2tblrace <- c2tblrace[,c(1:2,4)]
chisq.test(c2tblrace, correct=F)## Warning in chisq.test(c2tblrace, correct = F): Chi-squared approximation
## may be incorrect##
## Pearson's Chi-squared test
##
## data: c2tblrace
## X-squared = 5.6749, df = 2, p-value = 0.05857Subsets for aged and young subjects
c2aged <- subset(cohort2, group=="A")
c2young <- subset(cohort2, group=="Y")Exploratory boxplots for continuous variables
par(mfrow=c(2,2))
for (i in 4:7) {
y <- colnames(cohort2)[i]
variab <- cohort2[,i]
plot(variab~group, data=cohort2, ylab=y, xlab=NULL)
stripchart(variab~group, data=cohort2, vertical=T, method="jitter", add=T, pch=20, col="red")
}
Tests for equality of variances. All variances are equal by Levene’s test.
c2eqvartest <- apply(cohort2[4:7], 2, function(x) {leveneTest(x,cohort2$group)$'Pr(>F)'})[1,]
c2eqvartest<0.05## age height weight bmi
## FALSE FALSE FALSE FALSET-tests with output of p-values
c2ttests <- apply(cohort2[c(4:7)], 2, function(x) {t.test(x~group, data=cohort2, var.equal=T)$p.value})
c2ttests2 <- apply(cohort2[c(4:7)], 2, function(x) {t.test(x~group, data=cohort2, var.equal=T)})
round(c2ttests, 4)## age height weight bmi
## 0.0000 0.3749 0.7889 0.9800Calculations for means and SEs
c2meansyoung <- apply(c2young[c(4:7)], 2, function(x) {mean(x, na.rm=T)})
c2sesyoung <- apply(c2young[c(4:7)], 2, function(x) {std.error(x, na.rm=T)})
c2meansaged <- apply(c2aged[c(4:7)], 2, function(x) {mean(x, na.rm=T)})
c2sesaged <- apply(c2aged[c(4:7)], 2, function(x) {std.error(x, na.rm=T)})Means for continuous variables
c2means <- rbind(c2meansyoung, c2meansaged)
row.names(c2means) <- c("young","aged")
c2means## age height weight bmi
## young 27.54545 169.7273 74.78182 25.66364
## aged 66.55556 173.7778 77.25556 25.60000SEMs for continuous variables
c2ses <- rbind(c2sesyoung, c2sesaged)
row.names(c2ses) <- c("young", "aged")
c2ses## age height weight bmi
## young 1.337322 3.182597 7.533004 2.029387
## aged 1.415304 3.008219 3.977894 1.202775Metabolism Data
The following data are for metabolic assays by Seahorse extracellular flux analysis
Glycolysis Stress Test
Read in data
gst_data <- read.csv("gst.csv", header=T)
gst_ct <- "gst_counts.csv"Normalize data to cell number
gst_norm <- normXFp(gst_data, gst_ct)Group replicates
gst <- groupXFp(gst_norm, "Group")Calculate GST parameters
age <- gst_norm %>% distinct(Group, Age)
gst_param <- gstXFp(gst)
gst_param <- merge(gst_param, age, by="Group")Test for equality of variances. All variances are equal by Levene’s test.
eqvartest <- apply(gst_param[2:8], 2, function(x) {leveneTest(x,gst_param$Age)$'Pr(>F)'})[1,]
eqvartest<0.05 ## UnstimECAR GlucECAR OligoECAR DgECAR Glyc Capacity
## FALSE FALSE FALSE TRUE FALSE FALSE
## Reserve
## FALSEGlycolysis T-test:
t.test(Glyc~Age, data=gst_param)##
## Welch Two Sample t-test
##
## data: Glyc by Age
## t = -0.092878, df = 14.28, p-value = 0.9273
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -23.35284 21.41081
## sample estimates:
## mean in group A mean in group Y
## 60.18978 61.16079Glycolytic Capacity T-test
t.test(Capacity~Age, data=gst_param)##
## Welch Two Sample t-test
##
## data: Capacity by Age
## t = -0.68051, df = 13.582, p-value = 0.5076
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -30.00869 15.58448
## sample estimates:
## mean in group A mean in group Y
## 70.67801 77.89011Glycolytic Reserve T-test
t.test(Reserve~Age, data=gst_param)##
## Welch Two Sample t-test
##
## data: Reserve by Age
## t = -1.2931, df = 9.099, p-value = 0.2278
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -17.140972 4.658781
## sample estimates:
## mean in group A mean in group Y
## 10.48823 16.72932Calculate Means and SEMs
gst_desc <- gst_param %>%
group_by(Age) %>%
summarise(Glyc_mean=mean(Glyc),Glyc_sem=std.error(Glyc),
Capacity_mean=mean(Capacity), Capacity_sem=std.error(Capacity),
Reserve_mean=mean(Reserve), Reserve_sem=std.error(Reserve))
gst_desc## # A tibble: 2 × 7
## Age Glyc_mean Glyc_sem Capacity_mean Capacity_sem Reserve_mean
## <fctr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 A 60.18978 6.786212 70.67801 6.465194 10.48823
## 2 Y 61.16079 7.952837 77.89011 8.397748 16.72932
## # ... with 1 more variables: Reserve_sem <dbl>2 Hour LPS Assay
Read in data
lps_data <- read.csv("lps_data.csv", header=T)
lps_ct <- "lps_counts.csv"Normalize to cell number
lps_norm <- normXFp(lps_data, lps_ct)Group by triplicates
lps <- groupXFp(lps_norm, "Group")Age data frame for further calculations
age2 <- lps_norm %>% distinct(Group, Age)
age2$Cond <- as.factor(substr(age2$Group, 5, 15))
age2$Group <- as.factor(substr(age2$Group, 1, 3))
age2 <- filter(age2, Cond=="LPS Gluc")
age2 <- age2 %>% select(Group, Age)Break condition (LPS vs. MED) and group (AGED vs. YOUNG) into separate variables
lps$Cond <- as.factor(substr(lps$Group, 5, 15))
lps$Group <- as.factor(substr(lps$Group, 1, 3))Subset only relevant data (condition with glucose present)
gluc <- filter(lps, Cond=="LPS Gluc")
gluc <- merge(gluc, age2, by="Group")Group data by age
gluc2 <- groupXFp(gluc, "Age")Graph data by age group
graphAllXFp(gluc2, "GST")
Process data for statistics
gluc_tab <- tbl_df(gluc)
gluc_tab2 <- gluc_tab %>%
group_by(Group, Age, Cond) %>%
summarise(maxECAR=max(ECAR),
minECAR=min(ECAR),
diffECAR=maxECAR-minECAR,
AUC=AUC(Time, ECAR, method="trapezoid"))
gluc3 <- as.data.frame(select(gluc_tab2, Group, Age, Cond, minECAR, maxECAR, diffECAR, AUC))Test for equality of variances. All variances are equal by Levene’s test
eqvartest <- apply(gluc3[4:7], 2, function(x) {leveneTest(x,gluc3$Age)$'Pr(>F)'})[1,]
eqvartest<0.05 ## minECAR maxECAR diffECAR AUC
## FALSE FALSE FALSE FALSEDifference in ECAR (max - min) T-test
t.test(diffECAR ~ Age, data=gluc3) ##
## Welch Two Sample t-test
##
## data: diffECAR by Age
## t = -0.95443, df = 12.514, p-value = 0.3579
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -45.05295 17.51848
## sample estimates:
## mean in group A mean in group Y
## 64.11809 77.88532AUC T-test
t.test(AUC ~ Age, data=gluc3)##
## Welch Two Sample t-test
##
## data: AUC by Age
## t = 0.23832, df = 12.801, p-value = 0.8154
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -4998.909 6236.335
## sample estimates:
## mean in group A mean in group Y
## 15279.61 14660.9024 Hour LPS Assay
Read data
infl_data <- read.csv("infl_data.csv")
infl_ct <- "infl_counts.csv"Normalize to cell number
infl_norm <- normXFp(infl_data,infl_ct)Group by triplicate
infl <- groupXFp(infl_data, "Group")New variables to split condiion from subject ID and to add age
infl$Cond <- as.factor(substr(infl$Group, 1, 3))
infl$Group <- as.factor(substr(infl$Group, 5, 8))
infl$Age <- as.factor(substr(infl$Group, 1, 1))Group data by age
med <- subset(infl, Cond=="MED")
lps <- subset(infl, Cond=="LPS")
infl_med <- groupXFp(med, "Age")
infl_med$Cond <- rep("MED", nrow(infl_med))
infl_lps <- groupXFp(lps, "Age")
infl_lps$Cond <- rep("LPS", nrow(infl_lps))
infl2 <- rbind(infl_med, infl_lps)Calculate Means and SEMs
infl2 <- tbl_df(infl2)
infl2 <- infl2 %>%
group_by(Group, Cond) %>%
summarise(maxECAR=max(ECAR),
minECAR=min(ECAR),
seECAR=std.error(ECAR),
meanECAR=mean(ECAR))
infl2## Source: local data frame [4 x 6]
## Groups: Group [?]
##
## Group Cond maxECAR minECAR seECAR meanECAR
## <fctr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 A LPS 35.89942 22.96032 2.1153775 25.34420
## 2 A MED 25.46341 17.32251 1.2527595 19.33470
## 3 Y LPS 30.09795 22.03252 1.3092975 23.55762
## 4 Y MED 21.29150 17.66277 0.4811749 19.50801Statistics
infl <- tbl_df(infl)
infl3 <- infl %>%
group_by(Group, Cond, Age) %>%
summarise(maxECAR=max(ECAR),
minECAR=min(ECAR),
meanECAR=mean(ECAR))
infl_ANOVA <- aov(meanECAR~Cond*Age, data=infl3)
summary(infl_ANOVA)## Df Sum Sq Mean Sq F value Pr(>F)
## Cond 1 204.6 204.60 2.933 0.0987 .
## Age 1 4.1 4.06 0.058 0.8112
## Cond:Age 1 6.9 6.88 0.099 0.7559
## Residuals 26 1813.5 69.75
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Gene Expression Data
The following data are for cytokine gene expression from Seahorse cell lysates
2 Hour LPS Assay
Read Ct data
lps_gene <- read.csv("qpcr_data_lps.csv", header=T)
lps_gene_1 <- subset(lps_gene, lps_gene[,"Run"]==1)
lps_gene_2 <- subset(lps_gene, lps_gene[,"Run"]==2)Calculate dCts
dCt1 <- transmute(lps_gene_1, ID=ID, ID2=ID2, Group=Group, Run=Run,
IL10=IL10-B2M, IL6=IL6-B2M, IL12A=IL12A-B2M, TNF=TNF-B2M, IL1B=IL1B-B2M)
dCt2 <- transmute(lps_gene_2, ID=ID, ID2=ID2, Group=Group, Run=Run,
IL10=IL10-B2M, IL6=IL6-B2M, IL12A=IL12A-B2M, TNF=TNF-B2M, IL1B=IL1B-B2M)
dCt1$IL12A[is.na(dCt1$IL12A)] <- dCt2$IL12A[match(dCt1$ID, dCt2$ID)][which(is.na(dCt1$IL12A))]
dCt1$TNF[is.na(dCt1$TNF)] <- dCt2$TNF[match(dCt1$ID, dCt2$ID)][which(is.na(dCt1$TNF))]
dCt1$IL1B[is.na(dCt1$IL1B)] <- dCt2$IL1B[match(dCt1$ID, dCt2$ID)][which(is.na(dCt1$IL1B))]Calculate ddCts by 2^(-ddCt) method
young <- subset(dCt1, dCt1[,"Group"]=="Y")
hk2 <- apply(young[,5:9], 2, mean)
ddCt2 <- transmute(dCt1, ID=ID, ID2=ID2, Group=Group, Run=Run,
IL10=2^(-(IL10-hk2[1])), IL6=2^(-(IL6-hk2[2])),
IL12A=2^(-(IL12A-hk2[3])), TNF=2^(-(TNF-hk2[4])),
IL1B=2^(-(IL1B-hk2[5])))Test for equality of variances
eqvartest <- apply(ddCt2[5:9], 2, function(x) {leveneTest(x,ddCt2$Group)$'Pr(>F)'})[1,]
eqvartest<0.05## IL10 IL6 IL12A TNF IL1B
## FALSE FALSE FALSE FALSE FALSEIL10 T-test
t.test(IL10~Group, data=ddCt2)##
## Welch Two Sample t-test
##
## data: IL10 by Group
## t = 0.50089, df = 14.802, p-value = 0.6238
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -1.207331 1.947949
## sample estimates:
## mean in group A mean in group Y
## 2.042032 1.671722IL6 T-test
t.test(IL6~Group, data=ddCt2)##
## Welch Two Sample t-test
##
## data: IL6 by Group
## t = -0.93168, df = 15.987, p-value = 0.3654
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -1.2588506 0.4902049
## sample estimates:
## mean in group A mean in group Y
## 0.8236564 1.2079793IL12A T-test
t.test(IL12A~Group, data=ddCt2)##
## Welch Two Sample t-test
##
## data: IL12A by Group
## t = -0.47757, df = 9.9032, p-value = 0.6433
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -4.160283 2.693260
## sample estimates:
## mean in group A mean in group Y
## 1.836114 2.569626TNF T-test
t.test(TNF~Group, data=ddCt2)##
## Welch Two Sample t-test
##
## data: TNF by Group
## t = -1.4017, df = 15.727, p-value = 0.1804
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -1.2746000 0.2607995
## sample estimates:
## mean in group A mean in group Y
## 0.6906747 1.1975749IL1B T-test
t.test(IL1B~Group, data=ddCt2)##
## Welch Two Sample t-test
##
## data: IL1B by Group
## t = -0.66453, df = 10.974, p-value = 0.5201
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.9390071 0.5035838
## sample estimates:
## mean in group A mean in group Y
## 0.8675646 1.0852763Calculate summary statistics
ddCt2_tab <- tbl_df(ddCt2)
ddCt_tab2 <- ddCt2_tab %>%
na.omit() %>%
group_by(Group) %>%
summarise(IL1Bse=std.error(IL1B),
IL1Bmean=mean(IL1B),
IL6se=std.error(IL6),
IL6mean=mean(IL6),
IL10se=std.error(IL10),
IL10mean=mean(IL10),
TNFse=std.error(TNF),
TNFmean=mean(TNF))
means2 <- ddCt_tab2 %>%
select(Group, IL1Bmean, IL6mean, IL10mean, TNFmean)
norm_means_A <- as.data.frame(means2[1,2:5])/as.data.frame(means2[2,2:5])Normalized means
norm_means_A## IL1Bmean IL6mean IL10mean TNFmean
## 1 0.7993952 0.6407549 1.098791 0.5236554SEMs
sems <- ddCt_tab2 %>%
select(Group, IL1Bse, IL6se, IL10se, TNFse)
sems## # A tibble: 2 × 5
## Group IL1Bse IL6se IL10se TNFse
## <fctr> <dbl> <dbl> <dbl> <dbl>
## 1 A 0.2872693 0.3210315 0.6302549 0.2604002
## 2 Y 0.1575156 0.2958658 0.4421980 0.272021824 Hour LPS Assay
Read data
infl_gene <- read.csv("qpcr_data_infl.csv", header=T)Calculate dCt
dCt <- transmute(infl_gene, ID=ID, ID2=ID2, Group=Group, Cond=Cond, Set=Set,
IL1B=IL1B-B2M, IL6=IL6-B2M, IL10=IL10-B2M, IL12A=IL12A-B2M, TNF=TNF-B2M)Calculate ddCt by 2^(-ddCt) method
young_med <- subset(dCt, dCt[,"Group"]=="Y" & dCt[,"Cond"]=="Med")
hk <- apply(young_med[,6:10], 2, function(x){mean(x, na.rm=T)})
ddCt <- transmute(dCt, ID=ID, ID2=ID2, Group=Group, Cond=Cond, Set=Set,
IL1B=2^(-(IL1B-hk[1])), IL6=2^(-(IL6-hk[2])),
IL10=2^(-(IL10-hk[3])), IL12A=2^(-(IL12A-hk[4])),
TNF=2^(-(TNF-hk[5])))Test for equality of variances
eqvartest <- apply(ddCt[6:10], 2, function(x) {leveneTest(x,ddCt$Group)$'Pr(>F)'})[1,]
eqvartest<0.05 ## IL1B IL6 IL10 IL12A TNF
## FALSE FALSE FALSE FALSE FALSEIL1B ANOVA
il1b <- aov(IL1B~Group*Cond, data=ddCt)
summary(il1b)## Df Sum Sq Mean Sq F value Pr(>F)
## Group 1 25.27 25.27 4.225 0.05190 .
## Cond 1 76.49 76.49 12.786 0.00169 **
## Group:Cond 1 5.40 5.40 0.902 0.35249
## Residuals 22 131.60 5.98
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1IL6 ANOVA
il6 <- aov(IL6~Group*Cond, data=ddCt)
summary(il6)## Df Sum Sq Mean Sq F value Pr(>F)
## Group 1 130.5 130.48 5.912 0.02364 *
## Cond 1 199.6 199.60 9.044 0.00648 **
## Group:Cond 1 77.9 77.91 3.530 0.07358 .
## Residuals 22 485.6 22.07
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1IL10 ANOVA
il10 <- aov(IL10~Group*Cond, data=ddCt)
summary(il10)## Df Sum Sq Mean Sq F value Pr(>F)
## Group 1 31.33 31.33 3.620 0.0703 .
## Cond 1 113.73 113.73 13.140 0.0015 **
## Group:Cond 1 11.29 11.29 1.305 0.2657
## Residuals 22 190.41 8.66
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1IL12A ANOVA
il12a <- aov(IL12A~Group*Cond, data=ddCt)
summary(il12a)## Df Sum Sq Mean Sq F value Pr(>F)
## Group 1 0.103 0.1029 0.135 0.717
## Cond 1 0.001 0.0014 0.002 0.967
## Group:Cond 1 0.271 0.2715 0.355 0.558
## Residuals 20 15.275 0.7637
## 2 observations deleted due to missingnessTNF ANOVA
tnf <- aov(TNF~Group*Cond, data=ddCt)
summary(tnf)## Df Sum Sq Mean Sq F value Pr(>F)
## Group 1 0.171 0.1713 0.812 0.3773
## Cond 1 0.782 0.7818 3.706 0.0673 .
## Group:Cond 1 0.155 0.1552 0.736 0.4003
## Residuals 22 4.642 0.2110
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Means and SEMs
Y_med <- subset(ddCt, ddCt[,"Set"]==3)
Y_lps <- subset(ddCt, ddCt[,"Set"]==4)
A_med <- subset(ddCt, ddCt[,"Set"]==1)
A_lps <- subset(ddCt, ddCt[,"Set"]==2)
infl_means_young_med <- apply(Y_med[c(6:10)], 2, function(x) {mean(x, na.rm=T)})
infl_ses_young_med <- apply(Y_med[c(6:10)], 2, function(x) {std.error(x, na.rm=T)})
infl_means_aged_med <- apply(A_med[c(6:10)], 2, function(x) {mean(x, na.rm=T)})
infl_ses_aged_med <- apply(A_med[c(6:10)], 2, function(x) {std.error(x, na.rm=T)})
infl_means_young_lps <- apply(Y_lps[c(6:10)], 2, function(x) {mean(x, na.rm=T)})
infl_ses_young_lps <- apply(Y_lps[c(6:10)], 2, function(x) {std.error(x, na.rm=T)})
infl_means_aged_lps <- apply(A_lps[c(6:10)], 2, function(x) {mean(x, na.rm=T)})
infl_ses_aged_lps <- apply(A_lps[c(6:10)], 2, function(x) {std.error(x, na.rm=T)})
means <- rbind(infl_means_young_med, infl_means_aged_med, infl_means_young_lps,
infl_means_aged_lps)
norm_means <- sweep(means, 2, infl_means_young_med, `/`)Normalized means
row.names(norm_means) <- c("Young Med", "Aged Med", "Young LPS", "Aged LPS")
norm_means## IL1B IL6 IL10 IL12A TNF
## Young Med 1.0000000 1.000000 1.0000000 1.0000000 1.0000000
## Aged Med 0.2850921 0.280970 0.3921155 0.7193305 0.7278581
## Young LPS 3.9670746 7.534031 4.8721800 0.8324078 0.5601089
## Aged LPS 2.0236257 1.925715 2.4383328 0.8824156 0.5533877SEMs
sems2 <- rbind(infl_ses_young_med, infl_ses_aged_med, infl_ses_young_lps, infl_ses_aged_lps)
row.names(sems2) <- c("Young Med", "Aged Med", "Young LPS", "Aged LPS")
sems2## IL1B IL6 IL10 IL12A TNF
## Young Med 0.4597683 0.4069962 0.6375942 0.4233705 0.32252249
## Aged Med 0.2086043 0.2043841 0.1506255 0.3323634 0.14546899
## Young LPS 1.8571096 3.8942283 2.1040103 0.3064862 0.07702133
## Aged LPS 0.6895258 0.7539913 1.0288392 0.3637754 0.10399706