Introduction

This code examines the association of five genes with age (KRT14, CCL21, CD1a, CD3E, CXCL12) adjusted by KRT8. We used univariate analysis of variance (ANOVA) to evaluate the association.

The associated file “GeneExpressionFileCreation.rmd” replaces expression values that are below the limit of detection (LOD) for each gene with 1/2 the LOD. It also cleans up names and id’s so that it can be merged with the subject characteristics file (e.g., age and percent thymic epithelium levels)

Read in subject characteristics data

require(dplyr)
require(plyr)
require(rms)
require(readr)

raw.data <- read_csv("Quantgene2.csv")

## data formats

final.data <- as.data.frame(raw.data) %>% 
  mutate(age.ch=as.factor(ifelse(`Age of Thymus Tissue Donor (in Years)`<=17, "Young","Adult"))) %>% # categorizing age to two groups in age.ch variable
  mutate(age.ch3=as.factor(ifelse(`Age of Thymus Tissue Donor (in Years)`<=1, "Postnatal", 
                                  ifelse(`Age of Thymus Tissue Donor (in Years)`>1 & `Age of Thymus Tissue Donor (in Years)` <=17,"Child-Adolescent","Adult")))) %>% 
  # categorizing age to two or three groups in age.ch and age.ch3 variable
  rename(replace = c("Total area of relevant tissue on slide:" = "relevant.tissue")) %>%
  rename(replace = c("Total area of lymphoid tissue on slide:" = "lymphoid.tissue")) %>%
  rename(replace = c("Total area of thymic epithelium on slide:" = "thymic.epithelium")) %>%
  rename(replace = c("Total active cortical area (thymopoiesis):" = "cortical")) %>%
  rename(replace=c("DHVI Thymus Tissue ID Number (T#)"="id")) %>%
  rename(replace =c("Age of Thymus Tissue Donor (in Years)"="Age"))

names(final.data) <- gsub(" \\(human\\)","",names(final.data))
names(final.data) <- gsub(" ","",names(final.data))

#create \% thymic epithelial and \% cortical 
area.to.total <- final.data %>%
  mutate(lymphoid.tissue.p=100*lymphoid.tissue/relevant.tissue,
         thymic.epithelium.p=100*thymic.epithelium/relevant.tissue,
         cortical.p=100*cortical/relevant.tissue) %>% 
  select(id, Age,age.ch,age.ch3,Gender,
         lymphoid.tissue.p, thymic.epithelium.p, cortical.p, 
         everything(), 
         -relevant.tissue, -lymphoid.tissue,-thymic.epithelium, -cortical)

area.to.total$age.ch3 <- factor(area.to.total$age.ch3, levels=c("Postnatal", "Child-Adolescent", "Adult"))

Read in gene expression data

new.expression <-read.csv("new.expression.csv", stringsAsFactors=F)

### P10 = P1, 9086= 908 missing 244 -- can I assume 244 was a typo in the original clinical   #####

sort(new.expression$id)
##  [1] "1013" "1080" "158"  "2000" "2001" "2002" "2003" "222"  "240"  "263" 
## [11] "264"  "265"  "278"  "309"  "321"  "349"  "359"  "545"  "562"  "577" 
## [21] "587"  "588"  "591"  "598"  "628"  "630"  "639"  "680"  "697"  "707" 
## [31] "717"  "908"  "921"  "924"  "926"  "960"  "997"  "P1"   "P11"  "P12" 
## [41] "P13"  "P14"  "P15"  "P16"  "P17"  "P18"  "P19"  "P2"   "P20"  "P21" 
## [51] "P22"  "P23"  "P24"  "P25"  "P4"   "P5"   "P6"   "P7"   "P8"   "P9"
sort(area.to.total$id)
##  [1] "1013" "1080" "158"  "2000" "2001" "2002" "2003" "222"  "240"  "244" 
## [11] "265"  "278"  "309"  "321"  "349"  "359"  "545"  "562"  "577"  "587" 
## [21] "588"  "591"  "598"  "628"  "630"  "639"  "680"  "697"  "707"  "717" 
## [31] "908"  "921"  "924"  "926"  "960"  "997"  "P1"   "P11"  "P13"  "P17" 
## [41] "P2"   "P20"  "P21"  "P22"  "P24"  "P4"   "P8"
names(area.to.total)
##  [1] "id"                  "Age"                 "age.ch"             
##  [4] "age.ch3"             "Gender"              "lymphoid.tissue.p"  
##  [7] "thymic.epithelium.p" "cortical.p"          "VCAM1"              
## [10] "IL12a"               "PPARG"               "KRT14"              
## [13] "LIF"                 "AIRE"                "NCAM1"              
## [16] "CCL25"               "PSMB11"              "FLT3LG"             
## [19] "FGF7"                "BMP4"                "IL6"                
## [22] "ADIPOQ"              "KRT8"                "IL22"               
## [25] "CD1a"                "GAPDH"               "IL23a"              
## [28] "KIT"                 "CCL19"               "ICAM1"              
## [31] "DLL4"                "IL7"                 "FOXP3"              
## [34] "ATF3"                "S1PR1"               "ESR1"               
## [37] "CCL21"               "AR"                  "IN3C1"              
## [40] "GUSB"                "HPRT1"               "RORC"               
## [43] "TNFSF11"             "CD3E"                "WNTF4"              
## [46] "ITGAX"               "KITLG"               "NOTCH1"             
## [49] "IL2"                 "FOXN1"               "CCND1"              
## [52] "CXCL12"              "GREM1"               "CD80"               
## [55] "JAG1"                "TNFRSF11b"           "LTB"                
## [58] "ITH3"                "CYP21a2"             "SELP"
dim(area.to.total)
## [1] 47 60
area.to.total2 <-subset(area.to.total, select = 1:8)

updated.data <-merge(area.to.total2, new.expression, by="id")

updated.data$agenew <-ifelse(updated.data$Age >= 60,"60 +", updated.data$age.ch)
knitr::kable(table(updated.data$agenew))
Var1 Freq
1 16
2 25
60 + 5
knitr::kable(table(updated.data$Gender))
Var1 Freq
Female 27
Male 18
Unknown 1

GAPDH expression decreases with age, demonstrating that it is not a good standardization gene. Also, examines univariate association with other genes.

###  Show that GAPDH expression changes with age

GAPDH <- lm(log(GAPDH) ~  age.ch, data=updated.data)
summary(GAPDH)
## 
## Call:
## lm(formula = log(GAPDH) ~ age.ch, data = updated.data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.8009 -0.7738 -0.4675  0.4242  3.9292 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   5.3965     0.3172  17.011   <2e-16 ***
## age.chYoung   0.9160     0.4303   2.129   0.0389 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.454 on 44 degrees of freedom
## Multiple R-squared:  0.09337,    Adjusted R-squared:  0.07276 
## F-statistic: 4.531 on 1 and 44 DF,  p-value: 0.03892
GAPDH.mean <- lm(log(GAPDH) ~  -1 + age.ch, data=updated.data)
summary(GAPDH.mean)
## 
## Call:
## lm(formula = log(GAPDH) ~ -1 + age.ch, data = updated.data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.8009 -0.7738 -0.4675  0.4242  3.9292 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## age.chAdult   5.3965     0.3172   17.01   <2e-16 ***
## age.chYoung   6.3124     0.2907   21.71   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.454 on 44 degrees of freedom
## Multiple R-squared:  0.9453, Adjusted R-squared:  0.9428 
## F-statistic: 380.4 on 2 and 44 DF,  p-value: < 2.2e-16
### Does KRT14 show similar difference in age?

KRT14 <- lm(log(KRT14) ~  age.ch, data=updated.data)
summary(KRT14)
## 
## Call:
## lm(formula = log(KRT14) ~ age.ch, data = updated.data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.9581 -0.8388  0.0583  0.9165  3.4747 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   6.2162     0.2927  21.236   <2e-16 ***
## age.chYoung   0.4154     0.3971   1.046    0.301    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.341 on 44 degrees of freedom
## Multiple R-squared:  0.02427,    Adjusted R-squared:  0.002094 
## F-statistic: 1.094 on 1 and 44 DF,  p-value: 0.3012
KRT14means <-lm(log(KRT14) ~  -1 + age.ch, data=updated.data)
summary(KRT14means)
## 
## Call:
## lm(formula = log(KRT14) ~ -1 + age.ch, data = updated.data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.9581 -0.8388  0.0583  0.9165  3.4747 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## age.chAdult   6.2162     0.2927   21.24   <2e-16 ***
## age.chYoung   6.6316     0.2683   24.72   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.341 on 44 degrees of freedom
## Multiple R-squared:  0.9602, Adjusted R-squared:  0.9584 
## F-statistic:   531 on 2 and 44 DF,  p-value: < 2.2e-16
##### Does KRT8 show similar difference in age?


KRT8<- lm(log(KRT8) ~  age.ch, data=updated.data)
summary(KRT8)
## 
## Call:
## lm(formula = log(KRT8) ~ age.ch, data = updated.data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.8360 -0.5644 -0.2303  0.2739  3.0758 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   5.0665     0.1879   26.97   <2e-16 ***
## age.chYoung   0.2879     0.2548    1.13    0.265    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8608 on 44 degrees of freedom
## Multiple R-squared:  0.0282, Adjusted R-squared:  0.006114 
## F-statistic: 1.277 on 1 and 44 DF,  p-value: 0.2646
KRT8means <-lm(log(KRT8) ~  -1 + age.ch, data=updated.data)
summary(KRT8means)
## 
## Call:
## lm(formula = log(KRT8) ~ -1 + age.ch, data = updated.data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.8360 -0.5644 -0.2303  0.2739  3.0758 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## age.chAdult   5.0665     0.1879   26.97   <2e-16 ***
## age.chYoung   5.3545     0.1722   31.10   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8608 on 44 degrees of freedom
## Multiple R-squared:  0.9747, Adjusted R-squared:  0.9735 
## F-statistic: 847.3 on 2 and 44 DF,  p-value: < 2.2e-16
########################

KRT8 <- lm(log(KRT8/KRT14) ~ age.ch, data=updated.data)
anova(KRT8)
ABCDEFGHIJ0123456789
 
 
Df
<int>
Sum Sq
<dbl>
Mean Sq
<dbl>
F value
<dbl>
Pr(>F)
<dbl>
age.ch10.18542450.18542450.15282520.6977365
Residuals4453.38568941.2133111NANA
####################

CD1a <- lm(log(CD1a/KRT14) ~ age.ch, data=updated.data)
anova(CD1a)
ABCDEFGHIJ0123456789
 
 
Df
<int>
Sum Sq
<dbl>
Mean Sq
<dbl>
F value
<dbl>
Pr(>F)
<dbl>
age.ch12.7984222.7984222.195580.1455339
Residuals4456.0811091.274571NANA
####################

CCL21 <- lm(log(CCL.21/ KRT14) ~ age.ch, data=updated.data)
anova(CCL21)
ABCDEFGHIJ0123456789
 
 
Df
<int>
Sum Sq
<dbl>
Mean Sq
<dbl>
F value
<dbl>
Pr(>F)
<dbl>
age.ch14.7021754.7021752.4814260.1223615
Residuals4483.3777471.894949NANA
####################

CD3E <- lm(log(CD3.E/KRT14) ~ age.ch, data=updated.data)
anova(CD3E)
ABCDEFGHIJ0123456789
 
 
Df
<int>
Sum Sq
<dbl>
Mean Sq
<dbl>
F value
<dbl>
Pr(>F)
<dbl>
age.ch118.5315818.531579729.253672.474859e-06
Residuals4427.873070.6334788NANA
####################

CXCL12 <- lm(log(CXCL.12/KRT14) ~ age.ch, data=updated.data)
anova(CXCL12)
ABCDEFGHIJ0123456789
 
 
Df
<int>
Sum Sq
<dbl>
Mean Sq
<dbl>
F value
<dbl>
Pr(>F)
<dbl>
age.ch10.9294260.92942600.9985080.3231357
Residuals4440.9558530.9308148NANA

Analysis for the manuscript, genes standardized by KRT8.

########################

KRT14 <- lm(log(KRT14/KRT8) ~ age.ch, data=updated.data)
anova(KRT14)
ABCDEFGHIJ0123456789
 
 
Df
<int>
Sum Sq
<dbl>
Mean Sq
<dbl>
F value
<dbl>
Pr(>F)
<dbl>
age.ch10.18542450.18542450.15282520.6977365
Residuals4453.38568941.2133111NANA
KRT14.means <- lm(KRT14/KRT8 ~ -1 + age.ch, data=updated.data)
summary(KRT14.means)
## 
## Call:
## lm(formula = KRT14/KRT8 ~ -1 + age.ch, data = updated.data)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -5.643 -2.928 -1.246  1.486 13.529 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## age.chAdult   5.8017     0.9686   5.990 3.49e-07 ***
## age.chYoung   4.7378     0.8878   5.337 3.15e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.439 on 44 degrees of freedom
## Multiple R-squared:  0.5939, Adjusted R-squared:  0.5755 
## F-statistic: 32.18 on 2 and 44 DF,  p-value: 2.449e-09
####################

CD1a.krt8 <- lm(log(CD1a/KRT8) ~ age.ch, data=updated.data)
anova(CD1a.krt8)
ABCDEFGHIJ0123456789
 
 
Df
<int>
Sum Sq
<dbl>
Mean Sq
<dbl>
F value
<dbl>
Pr(>F)
<dbl>
age.ch14.4245354.42453495.2417840.02690493
Residuals4437.1399420.8440896NANA
CD1a.krt8.means <- lm(CD1a/KRT8 ~ -1 + age.ch, data=updated.data)
summary(CD1a.krt8.means)
## 
## Call:
## lm(formula = CD1a/KRT8 ~ -1 + age.ch, data = updated.data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.9927 -1.2898 -0.8416  0.8922  7.9366 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## age.chAdult   1.7287     0.5108   3.384  0.00151 ** 
## age.chYoung   3.1636     0.4682   6.758 2.59e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.341 on 44 degrees of freedom
## Multiple R-squared:  0.5649, Adjusted R-squared:  0.5451 
## F-statistic: 28.56 on 2 and 44 DF,  p-value: 1.121e-08
####################

CCL21.krt8 <- lm(log(CCL.21/ KRT8) ~ age.ch, data=updated.data)
anova(CCL21.krt8)
ABCDEFGHIJ0123456789
 
 
Df
<int>
Sum Sq
<dbl>
Mean Sq
<dbl>
F value
<dbl>
Pr(>F)
<dbl>
age.ch13.0200893.0200893.037360.08835384
Residuals4443.7498150.994314NANA
CCL21.krt8.means <- lm(CCL.21/ KRT8 ~ -1 + age.ch, data=updated.data)
summary(CCL21.krt8.means)
## 
## Call:
## lm(formula = CCL.21/KRT8 ~ -1 + age.ch, data = updated.data)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -3.984 -2.737 -0.589  0.212 38.194 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## age.chAdult    4.264      1.360   3.135  0.00306 **
## age.chYoung    1.128      1.247   0.905  0.37064   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.233 on 44 degrees of freedom
## Multiple R-squared:  0.1948, Adjusted R-squared:  0.1582 
## F-statistic: 5.322 on 2 and 44 DF,  p-value: 0.008511
####################

CD3E.krt8 <- lm(log(CD3.E/KRT8) ~ age.ch, data=updated.data)
anova(CD3E.krt8)
ABCDEFGHIJ0123456789
 
 
Df
<int>
Sum Sq
<dbl>
Mean Sq
<dbl>
F value
<dbl>
Pr(>F)
<dbl>
age.ch122.4244122.42440722.188232.493195e-05
Residuals4444.468341.010644NANA
CD3E.krt8.means <- lm(CD3.E/KRT8 ~ -1 + age.ch, data=updated.data)
summary(CD3E.krt8.means)
## 
## Call:
## lm(formula = CD3.E/KRT8 ~ -1 + age.ch, data = updated.data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.2151 -2.2259 -0.9735  1.2725  9.7558 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## age.chAdult   2.6502     0.8181   3.240  0.00228 ** 
## age.chYoung   6.9412     0.7498   9.258 6.82e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.749 on 44 degrees of freedom
## Multiple R-squared:  0.6862, Adjusted R-squared:  0.6719 
## F-statistic:  48.1 on 2 and 44 DF,  p-value: 8.463e-12
####################

CXCL12.krt8 <- lm(log(CXCL.12/KRT8) ~ age.ch, data=updated.data)
anova(CXCL12.krt8)
ABCDEFGHIJ0123456789
 
 
Df
<int>
Sum Sq
<dbl>
Mean Sq
<dbl>
F value
<dbl>
Pr(>F)
<dbl>
age.ch10.28457740.28457741.7313260.1950563
Residuals447.23226230.1643696NANA
CXCL12.krt8.means <- lm(CXCL.12/KRT8 ~ -1 + age.ch, data=updated.data)
summary(CXCL12.krt8.means)
## 
## Call:
## lm(formula = CXCL.12/KRT8 ~ -1 + age.ch, data = updated.data)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -4.986 -2.320 -0.350  2.491  5.035 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## age.chAdult   7.5520     0.5866   12.87  < 2e-16 ***
## age.chYoung   6.2209     0.5376   11.57 6.16e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.688 on 44 degrees of freedom
## Multiple R-squared:  0.872,  Adjusted R-squared:  0.8661 
## F-statistic: 149.8 on 2 and 44 DF,  p-value: < 2.2e-16

plots that correspond to the analysis of variance.

library(beanplot)
library(tidyr)

updated.data.nogapdh <-subset(updated.data, select=-c(GAPDH))

updated.data.long <- updated.data.nogapdh %>% gather(Gene, Expression, c(KRT14, CD1a:CXCL.12))


#png("beanplot expression.png")
beanplot(log(Expression/KRT8)~age.ch*Gene, data=updated.data.long,  ll = 0.04,
         side = "both", xlab="", ylab="Log Ratio to KRT8", cex.axis=0.75,
         col = list("grey56", c("grey92", "black")),
         axes=T)

#dev.off()

library(dplyr) 
library(ggplot2)
library(ggpubr)
theme_set(theme_pubclean())

updated.data.long$Expression.adjusted <-log(updated.data.long$Expression/updated.data.long$KRT8)

#png("boxplot expression.png")
ggplot(updated.data.long, aes(Gene, Expression.adjusted)) +
  geom_boxplot(aes(color = age.ch))+ labs(y= "Log Expression Standardized to KRT8") +
  scale_color_manual(values = c("black", "grey41")) +
  stat_compare_means(aes(group = age.ch), label = "p.signif")

#dev.off()
## 
##   0 4.9   7   9  16  17  22  27  28  29  32  34  36  38  39  40  45  49  53  58 
##  20   1   1   1   1   1   1   1   1   1   1   2   1   1   2   1   1   1   1   1 
##  59  60  64  67  70  78 
##   1   1   1   1   1   1
## 
## Call:
## lm(formula = thymic.epithelium.p ~ Gender, data = updated.data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -56.727 -28.767   3.004  30.261  46.019 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     43.274      6.488   6.670 3.88e-08 ***
## GenderMale      15.897     10.258   1.550    0.129    
## GenderUnknown    9.416     34.331   0.274    0.785    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 33.71 on 43 degrees of freedom
## Multiple R-squared:  0.05306,    Adjusted R-squared:  0.009021 
## F-statistic: 1.205 on 2 and 43 DF,  p-value: 0.3097
## 
## Call:
## lm(formula = thymic.epithelium.p ~ age.zip, data = area.to.total)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -43.357  -9.181   3.026  11.506  23.109 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 76.06707    3.27195   23.25  < 2e-16 ***
## age.zip     -1.21148    0.09884  -12.26 6.12e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 16.49 on 45 degrees of freedom
## Multiple R-squared:  0.7695, Adjusted R-squared:  0.7644 
## F-statistic: 150.2 on 1 and 45 DF,  p-value: 6.123e-16
## 
## Call:
## lm(formula = thymic.epithelium.p ~ Age, data = non.zero)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -36.896  -9.717   0.108  12.674  34.090 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  61.6117     7.0463   8.744 4.48e-09 ***
## Age          -0.9154     0.1613  -5.674 6.59e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 16.36 on 25 degrees of freedom
## Multiple R-squared:  0.5629, Adjusted R-squared:  0.5454 
## F-statistic: 32.19 on 1 and 25 DF,  p-value: 6.592e-06
## 
## Call:
## lm(formula = thymic.epithelium.p ~ Age, data = area.to.total)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -43.436  -9.201   2.979  11.495  22.975 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 76.21779    3.28803   23.18  < 2e-16 ***
## Age         -1.21414    0.09932  -12.22 6.72e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 16.52 on 45 degrees of freedom
## Multiple R-squared:  0.7685, Adjusted R-squared:  0.7634 
## F-statistic: 149.4 on 1 and 45 DF,  p-value: 6.724e-16
## png 
##   2
## png 
##   2