Revision 8575912e5aa73fde6ab54bba4fbaeadc7d8267f4 authored by Amin Haghani on 20 November 2024, 18:45:15 UTC, committed by GitHub on 20 November 2024, 18:45:15 UTC
1 parent 276d6b4
Codes_for_figures.Rmd
---
title: "Codes for Figures"
author: "Zhe Fei"
date: "2023-06-12"
output: html_document
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(tidyverse)
```
## Figure 1: Dog AROCM
```{r}
datN12all = read.csv("/Users/feiz/Dropbox/N12.2018-9238DogBloodOstrander/SampleSheetAgeN12final.csv")
datN12 = datN12all%>%
filter(ConfidenceInAgeEstimate>= 90 &
CanBeUsedForAgingStudies == "yes")
colnames(datN12)
length(unique(datN12$Breed.Index))
tabDogBreed = sort(table(datN12$DogBreedShort))
table(datN12$Age/datN12$LifespanMedianHorvath/1.5 >1)
table(datN12$Age/datN12$LifespanUpperHorvath >1)
View(datN12%>%filter(Age> 1.5*LifespanMedianHorvath))
range(datN12$Age/datN12$LifespanUpperHorvath)
range(datN12$Age/datN12$LifespanMedianHorvath)
# LifespanMedianHorvath
# LifespanUpperHorvath
# Weight.kg.avg
```
```{r}
m1 = match(datN12$Basename, colnames(dat0sesame))
summary(m1)
```
### fit slopes
```{r}
source("Codes_AROCM/0_fns_v2.R")
vnum = "dog_v5"
cex.main=1.5
cut1 = freq = 3
j=2
# for (j in 1:length(statesPRC2))
head(statesPRC2)
```
```{r}
{
cgid = cg_list[[j]]
len1 = statesPRC2[j]
dogslopes = data.frame(
matrix(NA, nrow = length(tabDogBreed), ncol = 10)
)
for(nb in 1:length(tabDogBreed)){
idx = which(datN12$DogBreedShort == names(tabDogBreed)[nb])
dogslopes[nb,1] = names(tabDogBreed)[nb]
dogslopes[nb,2] = length(idx)
dogslopes[nb,3] = datN12$LifespanMedianHorvath[idx[1]]
dogslopes[nb,4] = datN12$LifespanUpperHorvath[idx[1]]
dogslopes[nb,5] = datN12$Weight.kg.avg[idx[1]]
basename = datN12$Basename[idx]
age1 = datN12$Age[idx]
m1 = match(basename, colnames(dat0sesame))
cgmean = colMeans(dat0sesame[cgid,m1])
dogslopes[nb,6:10] = fit1 = AROCM(cgmean, age1)
}
colnames(dogslopes) = c("DogBreed", "Freq", "LifespanMedianHorvath",
"LifespanUpperHorvath", "Weight.kg.avg",
names(fit1))
}
```
```{r}
dogslopes$SD_Rage =
dogslopes$SD_Age/dogslopes$LifespanMedianHorvath/2
dogslopes$Breed.Index =
datN12$Breed.Index[match(dogslopes$DogBreed, datN12$DogBreedShort)]
tmpd = datN12%>%filter(Breed.Index == 67)
tmpd = datN12%>%filter(Breed.Index == 11)
range(tmpd$Age/tmpd$LifespanMedianHorvath/2)
```
```{r}
colnames(dogslopes)
dogplot = dogslopes%>%filter(R2> 0.2 & AROCM < 0.5) ## R2> 0.2 &
ggplot(dogplot,
aes(x = LifespanMedianHorvath, y = AROCM,
label = Breed.Index)) +
geom_text()
```
- check outliers
```{r}
nb = 90
idx = which(datN12$Breed.Index == nb)
basename = datN12$Basename[idx]
age1 = datN12$Age[idx]
m1 = match(basename, colnames(dat0sesame))
cgmean = colMeans(dat0sesame[cgid,m1])
plot(age1, cgmean)
```
### Adjusted AROCM
- AROCM vs. SD_Rage
```{r}
dogplot = dogslopes%>%filter(R2> 0.2 & AROCM < 0.5) ## R2> 0.2 &
ggplot(dogplot,
aes(x = SD_Rage, y = AROCM,
label = Breed.Index)) +
geom_text()
```
```{r}
pvec = c(-20:20)*0.05
CV = QCOD = Cor_SD_Rage = NULL
for (p in pvec) {
adjROC = dogplot$AROCM * (dogplot$SD_Rage)^(1 - p)
adjCor = dogplot$Cor / (dogplot$SD_Rage)^(p)
CV = c(CV, sd(adjCor)/mean(adjCor))
QCOD = c(QCOD, qcod(adjCor) )
Cor_SD_Rage = c(Cor_SD_Rage, cor(adjCor, dogplot$SD_Rage) )
}
adjpmat = data.frame(
pvec,
CV,
QCOD,
Cor_SD_Rage
)
write.csv(adjpmat, paste0("DogSlopes/adjROC_select_power_",
len1,"_",vnum, ".csv"))
```
- Optimal p and Adj.ROC, adj.cor
```{r}
popt = 0.1
dogslopes$adjROC = dogslopes$AROCM * (dogslopes$SD_Rage)^(1 - popt)
dogslopes$adjCor = dogslopes$Cor / (dogslopes$SD_Rage)^(popt)
# dogplot
```
### Plot slopes
- AROCM vs. Lifespan on log scale
```{r}
colnames(dogslopes)
dogplot = dogslopes%>%filter(R2> 0.2 & AROCM < 0.5)
cor1 = with(dogplot, cor.test(log(LifespanMedianHorvath), log(AROCM)))
p1 = ggplot(dogplot,
aes(log(LifespanMedianHorvath), log(AROCM))) +
geom_point() + ggtitle(paste0("Cor=",round(cor1$estimate,2),
", p=",signif(cor1$p.value,2)))
p1
```
#### AROCM vs. Lifespan
- regression:
```{r}
dogslopes$LifespanMedianScaled =
2 * dogslopes$LifespanMedianHorvath
colnames(dogslopes)
dogplot = dogslopes%>%filter(R2> 0.2 & AROCM < 0.5)
colnames(dogplot)
x = 1/dogplot$LifespanMedianScaled
y = dogplot$adjROC
lm0 = lm(y~x - 1)
summary(lm0)
plot(x,y)
```
#### AROCM vs. Lifespan
```{r}
cor1 = with(dogplot, cor.test(LifespanMedianScaled, AROCM))
p1_AROCM = ggplot(dogplot,
aes(LifespanMedianScaled, AROCM, label = Breed.Index) ) +
geom_text() + xlab("Lifespan") +
ggtitle(paste0("a. Cor=",round(cor1$estimate,2),
", p=",signif(cor1$p.value,2))) +
theme_bw()
p1_AROCM
```
#### AROCM vs. a/Lifespan
```{r}
cor1 = with(dogplot, cor.test(1/LifespanMedianScaled, AROCM))
x = 1/dogplot$LifespanMedianScaled
y = dogplot$AROCM
coef1 = coef(lm(y~x-1))
dogplot$X = coef1*x
xr = range(c(dogplot$X, dogplot$AROCM))
p1_AROCM_a = ggplot(dogplot,
aes(X, AROCM, label = Breed.Index) ) +
geom_text() + xlab(paste0("a/Lifespan, a=",round(coef1,2))) +
ggtitle(paste0("c. Cor=",round(cor1$estimate,2),
", p=",signif(cor1$p.value,2))) +
xlim(xr) + ylim(xr) +
geom_abline(intercept = 0,slope = 1, linetype=3) +
theme_bw()
p1_AROCM_a
```
#### Adj.AROCM vs. lifespan
```{r}
colnames(dogplot)
cor1 = with(dogplot, cor.test(LifespanMedianScaled, adjROC))
p1_adjroc = ggplot(dogplot,
aes(LifespanMedianScaled, adjROC, label = Breed.Index) ) +
geom_text() + xlab("Lifespan") + ylab("Adj.AROCM") +
ggtitle(paste0("b. Cor=",round(cor1$estimate,2),
", p=",signif(cor1$p.value,2))) +
theme_bw()
p1_adjroc
```
#### ADj. AROCM vs. a/Lifespan
```{r}
cor1 = with(dogplot, cor.test(1/LifespanMedianScaled, adjROC))
x = 1/dogplot$LifespanMedianScaled
y = dogplot$adjROC
coef1 = coef(lm(y~x-1))
dogplot$X = coef1*x
xr = range(c(dogplot$X, dogplot$adjROC))
p1_adjROC_a = ggplot(dogplot,
aes(X, adjROC, label = Breed.Index) ) +
geom_text() + ylab("Adj.AROCM") +
xlab(paste0("a/Lifespan, a=",round(coef1,2))) +
ggtitle(paste0("d. Cor=",round(cor1$estimate,2),
", p=",signif(cor1$p.value,2))) +
xlim(xr) + ylim(xr) +
geom_abline(intercept = 0,slope = 1, linetype=3) +
theme_bw()
p1_adjROC_a
```
#### AROCM vs. weight
```{r}
colnames(dogplot)
cor1 = with(dogplot, cor.test(Weight.kg.avg, AROCM))
p1_AROCM_wt = ggplot(dogplot,
aes(Weight.kg.avg, AROCM, label = Breed.Index) ) +
geom_text() + xlab("Weight") + ylab("AROCM") +
ggtitle(paste0("e. Cor=",round(cor1$estimate,2),
", p=",signif(cor1$p.value,2)))+
theme_bw()
p1_AROCM_wt
cor1 = with(dogplot, cor.test(Weight.kg.avg, adjROC))
p1_adjroc_wt = ggplot(dogplot,
aes(Weight.kg.avg, adjROC, label = Breed.Index) ) +
geom_text() + xlab("Weight") + ylab("Adj.AROCM") +
ggtitle(paste0("f. Cor=",round(cor1$estimate,2),
", p=",signif(cor1$p.value,2))) +
theme_bw()
p1_adjroc_wt
```
#### Dog main figure 1
```{r}
vnum = "v6"
pdf(paste0("DogSlopes/dogslopes_main_",len1,"_",vnum, ".pdf"),
onefile = TRUE,width = 6, height = 9)
grid.arrange(p1_AROCM,p1_adjroc,
p1_AROCM_a,p1_adjROC_a,
p1_AROCM_wt, p1_adjroc_wt,
nrow=3)
dev.off()
write.csv(dogslopes,
paste0("DogSlopes/dogslopes_",len1,"_",vnum, ".csv"))
dogBreedAROCM = dogplot%>%select(
DogBreed, AROCM, adjROC, LifespanMedianScaled
)
colnames(dogBreedAROCM)[3] = c("AdjAROCM")
write.csv(dogBreedAROCM, paste0("DogSlopes/dogBreedAROCM_",len1,"_",vnum, ".csv"))
```
## Figure 2: Dog Cor Bar Plot
### barplot Cor and Adj.Cor
```{r}
colnames(dogplot)
dogplot = dogplot[order(dogplot$LifespanMedianScaled),]
dogplot$rank = 1:nrow(dogplot)
dogplot$breedlab = paste0(dogplot$Breed.Index,". ",dogplot$DogBreed)
```
```{r}
cor.test(log(dogplot$LifespanMedianScaled), log(dogplot$Cor))
cor.test(log(dogplot$LifespanMedianScaled), log(dogplot$adjCor))
```
```{r}
pbar_cor = ggplot(dogplot) +
geom_bar(aes(y = Cor, x = rank),
stat="identity", position="dodge",
width = 1,
fill="grey", col="black",na.rm=TRUE) +
coord_flip() + theme_bw() +
geom_hline(yintercept = mean(dogplot$Cor), col="red") +
labs(y = "Cor(Age,Methyl)", title = "") +
scale_x_continuous("Breed", # limits = range(dogplot$rank),
breaks = dogplot$rank,
labels = dogplot$breedlab) +
# scale_y_continuous(expand = c(0,0)) +
theme(axis.text=element_text(size=8),
axis.ticks.x = element_blank(),
plot.title = element_text(size = 11, hjust = 0.5)) +
theme(legend.position = "none")
```
```{r}
pbar_adjcor = ggplot(dogplot) +
geom_bar(aes(y = adjCor, x = rank),
stat="identity", position="dodge",
width = 1,
fill="grey", col="black",na.rm=TRUE) +
coord_flip() + theme_bw() +
geom_hline(yintercept = mean(dogplot$adjCor), col="red") +
labs(y = "Adj.Cor(Age,Methyl)", title = "") +
theme(axis.text=element_text(size=8),
axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
plot.title = element_text(size = 11, hjust = 0.5)) +
theme(legend.position = "none")
pbar_adjcor
```
```{r}
pdf(paste0("DogSlopes/BarplotCor_",vnum,".pdf"),
width = 7, height = 10)
hlay = matrix(c(1,1,2,
1,1,2,
1,1,2),
nrow = 3 ,byrow = TRUE)
grid.arrange(pbar_cor, pbar_adjcor,
nrow=1,layout_matrix=hlay)
dev.off()
```
## Figure 3 & 4: All chromatin states and single species plots
```{r}
############ 1. summary state Slopes
stateSlopes = statesTab
dim(stateSlopes)
head(stateSlopes)
# stateSlopes$Cor.Slope.Identity=NULL
# head(statesTab)
length(statesPRC2)
```
- Different correlations for each state
```{r}
scaleCpG
colnames(SpeciesMat)[9:11]
if(scaleCpG){
idx1 = which(SpeciesMat$RemoveFei==0)
nc = c(10:12, grep("Slope",colnames(SpeciesMat)))
speciescor = combineTissues(SpeciesMat,idx1, nc)
name1 = "IdentityScaled"
}
dim(speciescor)
head(colnames(speciescor))
# name1
# idx1 = which(SpeciesMat$RemoveFei==0)
colnames(stateSlopes)
stateSlopes$CorSpearman = stateSlopes$CorLog =
stateSlopes$`Cor1/Lifespan` =
stateSlopes$`Cor1/Slope` = stateSlopes$NumberPositiveSlopes = NA
```
- correlations of each state
```{r}
pdf(paste0("../",outfolder,"/SpeciesStates_corplot.pdf"), onefile = TRUE)
par(mfrow=c(2,2))
ncut = 60
for(k in 1:length(statesPRC2)){
us = statesPRC2[k]
nc = grep(us, colnames(speciescor), fixed = TRUE)
slope1 = speciescor[,nc]
slope1 = slope1[,11] ### 10 or 11
col1 = grep(paste0(us,"OldSlope1"), colnames(SpeciesMat),
fixed = TRUE)[1]
tmpslope = SpeciesMat[SpeciesMat$RemoveFei==0,col1]
tmp1 = which(stateSlopes$statesPRC2 == us)
# print(tmp1)
x = speciescor$maxAgeCaesar
stateSlopes$NumberPositiveSlopes[tmp1] = npos = sum(tmpslope > 0 & !is.na(tmpslope))
# slope2 = slope1 - min(slope1,na.rm = T) + 1
if(npos < ncut ) {
slope2 = slope1
slope2[slope1>0] = NA
}
stateSlopes$CorLog[tmp1] =
ifelse(npos<ncut, cor(-log(x), log(-slope2),
use = "pairwise.complete"),
cor(log(x), log(slope1), use = "pairwise.complete") )
stateSlopes$CorSpearman[tmp1] =
cor(x, slope1, use = "pairwise.complete", method = "spearman")
stateSlopes$`Cor1/Lifespan`[tmp1] =
cor(1/x, slope1, use = "pairwise.complete")
#slope2 = slope1 - 1.5*min(slope1,na.rm = T)
stateSlopes$`Cor1/Slope`[tmp1] = ifelse(npos<ncut,
cor(x, 1/slope2, use = "pairwise.complete"),
cor(x[slope1>0], 1/slope1[slope1>0], use = "pairwise.complete"))
# cor(x[slope1>0], 1/slope1[slope1>0], use = "pairwise.complete")
main1 = paste0(us,"\n")
xlab = "maxAge"
ylab = "Slope"
verboseScatterplot(x, slope1, main = main1,
xlab=xlab,ylab=ylab)
verboseScatterplot(log(x), log(slope1), main = main1,
xlab=paste0("log ",xlab),ylab=paste0("log ",ylab))
slope3 = slope1
slope3[slope1<0] = NA
if(npos<ncut) y = 1/slope2 else y = 1/slope3
verboseScatterplot(x, y, main = main1,
xlab=xlab,ylab=paste0("1/",ylab))
verboseScatterplot(1/x, slope1, main = main1,
xlab=paste0("1/",xlab),ylab=ylab)
}
dev.off()
```
- all states combined
```{r}
# View(stateSlopes)
colnames(stateSlopes)
nc = 8
# stateSlopes[,nc] = -stateSlopes[,nc]
stateSlopes = stateSlopes[order(stateSlopes[,nc]),]
rbPal <- colorRampPalette(c('blue','red'))
rbPal(2)
#This adds a column of color values
# based on the y values
stateSlopes$Col <- rbPal(nrow(stateSlopes))
y = - stateSlopes[,nc]
head(sort(y),12)
tail(sort(y),12)
corl = -0.36
coru = 0.539
idx1 = which(y< corl |
(y> coru & stateSlopes$NumberPositiveSlopes>100))
idx1 = which(y< corl | (y> coru))
# head(stateSlopes$statesPRC2)
n1 = 6
n2 = 4
idx1 = c(1:n1, tail(1:nrow(stateSlopes),n2))
stateSlopes$statesPRC2[idx1]
write.csv(stateSlopes, paste0("../",outfolder,"/stateSlopes_",name1, ".csv"), row.names = 1:nrow(stateSlopes))
pdf(paste0("../",outfolder,"/stateSlopes_",name1, ".pdf"))
y = stateSlopes$Freq
x = -stateSlopes[,nc]
{
plot(x,y, pch=19, col = stateSlopes$Col, xlim=c(-1,1),
ylab = "# of CpGs", xlab = "Cor(AROCM, 1/LifeSpan)",
main = "CpG lists Correlations with lifespan")
text(x[idx1],y[idx1], labels = stateSlopes$statesPRC2[idx1],
pos=c(rep(4, n1),
rep(2, n2)),
cex=0.75)
abline(v= corl, col="blue", lty= 2)
abline(v= coru, col="red", lty= 2)
}
dev.off()
```
### Prepare Figure 1
```{r}
# View(SpeciesMat[,1:16]%>%filter(Tissue == "Blood"))
specs = SpeciesMat$MammalNumberHorvath[SpeciesMat$Tissue == "Blood" &
SpeciesMat$Freq>80]
## human skin
j=2
len1 = statesPRC2[j]
cgid = cg_list[[j]]
specs = c("9.2", "9.1", "4.25", "4.26", "5.12", "1.4","5.25","6.4", "6.2",
"3.23", "3.22", "3.21", "2.2","2.1",
"3.14", "1.1", "3.11")
specs = c("4.11.1", "5.1.1", "6.2.1", "1.4.1", "4.19.1", "2.1.1", "6.1.1", "1.1.1",
"2.1.2", "4.12.3", "9.9.1", "4.13.7", "5.9.1", "9.9.3", "4.7.1", "4.7.2", "4.13.8")
specs = c("Rattus norvegicus", "Mus musculus",
"Sus scrofa domesticus","Sus scrofa minusculus",
"Canis lupus familiaris", "Chlorocebus aethiops sabaeus",
"Ceratotherium simum simum", "Phoca vitulina",
"Delphinapterus leucas", "Equus caballus",
"Loxodonta africana", "Tursiops truncatus",
"Orcinus orca", "Elephas maximus",
"Megaptera novaeangliae", "Homo sapiens",
"Balaena mysticetus")
tiss = c(rep("Blood", 14),
rep("Skin", 3))
table(SpeciesMat$Tissue=="Skin")
specs = SpeciesMat$MammalNumberHorvath[SpeciesMat$Tissue == "Skin"]
tiss = rep("Skin", length(specs))
```
- Final Figure 1
```{r}
stateSlopes = read.csv("/Users/feiz/Documents/Box Sync/Horvath/Mammalian meth projects/Slope/MainFigures/stateSlopes_IdentityScaled.csv",row.names = 1)
vnum = "v8"
caseYO = 3 ## 1: Young - 0.1L, Old - ASM; 2: 0.1L; 3: ASM
# matratio = SpeciesMat$AvgMaturity/SpeciesMat$maxAgeCaesar
matratio = 0.1
statesfolder = "out_AROCM_0.1L"
```
```{r}
j=2
len1 = statesPRC2[j]
cgid = cg_list[[j]]
colnames(stateSlopes)
nc = 8
#pdf(paste0(outfolder,"/Figure1_",vnum, ".pdf"))
cex = 1.4
cex.main = cex.lab = 2
y = stateSlopes$Freq
x = -stateSlopes[,nc]
idx1 = which(stateSlopes[,nc]< corl |
stateSlopes[,nc]> coru)
idx1 = c(1:n1, tail(1:nrow(stateSlopes),n2))
stateSlopes$statelabel[idx1] = stateSlopes$statesPRC2[idx1]
range(stateSlopes$CorSpearman)
fig1a = ggplot(stateSlopes,aes(-CorSpearman,Freq,col=Col,label = statelabel))+
geom_point() + geom_text(hjust=0, vjust=0, size = 3) +
xlim(-0.5, 1.1) +
ylab("# of CpGs") + xlab("Spearman Cor(AROCM, 1/Lifespan)") +
ggtitle("a. Corr in 54 chromatin states")+
geom_vline(aes(xintercept=corl), colour="blue", linetype="dashed") +
geom_vline(aes(xintercept=coru), colour="red", linetype="dashed")+
theme_bw() + theme(legend.position = "")
write.csv(stateSlopes,
paste0("../",outfolder,"/stateSlopes_",vnum, ".csv"),
row.names = 1:nrow(stateSlopes))
```
### Figure 1b and 1c
```{r}
with(dat_scaledM_rage,cor(Rage, ScaledM))
rmcorr(Spec,Rage, ScaledM, dataset = dat_scaledM_rage)
cors2 = plyr::ddply(dat_scaledM_rage, "Spec",
function(mat1){
data.frame(
cor_rage = cor(mat1$Rage, mat1$ScaledM),
cor_tage = cor(mat1$tage, mat1$ScaledM),
cor_logliage = cor(mat1$logliage, mat1$ScaledM)
)
})
summary(cors2)
```
```{r}
lm_age = lm(ScaledM~ logliage, data = dat_scaledM_rage)
coef(lm_age)
p_logliage = ggplot(dat_scaledM_rage,aes(x=logliage, y = ScaledM,color=Spec)) +
# geom_text(aes(color=Spec, label = SpecNum),hjust=0, vjust=0, size = 1) +
geom_point() +
geom_smooth(method = lm,se = FALSE,linewidth=0.3,linetype=2)+
geom_abline(slope = coef(lm_age)[2], intercept = coef(lm_age)[1]) +
xlab("Logli RAge") + ggtitle("c. Median Cor=0.76") +
theme_bw() + theme(legend.position = "",axis.title.y=element_blank())
p_logli = ggplot(dat_scaledM_rage,aes(x=logliage, y = ScaledM,color=Spec)) +
geom_text(aes(color=Spec, label = SpecNum),hjust=0, vjust=0, size = 1) +
# geom_point() +
geom_smooth(method = lm,se = FALSE,linewidth=0.3,linetype=2)+
geom_abline(slope = coef(lm_age)[2], intercept = coef(lm_age)[1]) +
xlab("Logli Relative Age (N=10932)") + ggtitle("c. Median Cor=0.76, p=1e-16") +
theme_bw() + theme(legend.position = "",axis.title.y=element_blank())
lm_Rage = lm(ScaledM~Rage, data = dat_scaledM_rage)
coef(lm_Rage)
lfit = loess(ScaledM ~ Rage, data = dat_scaledM_rage)
dat_scaledM_rage$loessfit = lfit$fitted
p_rage = ggplot(dat_scaledM_rage,aes(x=Rage, y = ScaledM)) +
geom_point(aes(color=Spec)) +
## geom_text(aes(color=Spec, label = SpecNum),hjust=0, vjust=0, size = 1) +
geom_smooth(aes(x=Rage, y = ScaledM), data = dat_scaledM_rage,
formula = y~logli(x,m1=0.1),se=FALSE,linewidth=0.5,col="black") +
# geom_abline(slope = coef(lm_Rage)[2], intercept = coef(lm_Rage)[1]) +
xlab("Relative Age") + ggtitle("b. Median Cor=0.73") + theme_bw() +
theme(legend.position = "")
p_ra = ggplot(dat_scaledM_rage,aes(x=Rage, y = ScaledM)) +
## geom_point(aes(color=Spec)) +
geom_text(aes(color=Spec, label = SpecNum),hjust=0, vjust=0, size = 1) +
geom_smooth(aes(x=Rage, y = ScaledM), data = dat_scaledM_rage,
formula = y~log(x),se=FALSE,linewidth=0.5,col="black") +
# geom_abline(slope = coef(lm_Rage)[2], intercept = coef(lm_Rage)[1]) +
xlab("Relative Age (N=10932)") + ggtitle("b. Median Cor=0.73, p=1e-16") + theme_bw() +
theme(legend.position = "")
```
```{r, eval=FALSE}
with(dat_scaledM_rage,cor(tage, ScaledM))
rmcorr(Spec, tage, ScaledM, dataset = dat_scaledM_rage)
lm_tage = lm(ScaledM~tage, data = dat_scaledM_rage)
coef(lm_tage)
p_tage = ggplot(dat_scaledM_rage,aes(x=tage, y = ScaledM,
color=Spec, label = SpecNum)) +
## geom_point(size=0.5, alpha=0.5) +
geom_text(hjust=0, vjust=0, size = 1) +
geom_smooth(method = lm,se = FALSE,linewidth=0.3,linetype=2)+
geom_abline(slope = coef(lm_tage)[2], intercept = coef(lm_tage)[1]) +
xlab("Log-Log Relative Age") + ggtitle("b. Median Cor=0.74") + theme_bw() +
theme(legend.position = "")
colnames(dat_scaledM_rage)
```
### Single species
```{r}
dim(tslopes)
tslopes = as.data.frame(tslopes)
tslopes$Cor_diff = tslopes$Cor_TAge - tslopes$Cor_Age
```
```{r}
ss = c("9.9.1", "6.1.1", "1.1.1",
"12.3.10", "9.9.1", "9.1.3",
"1.4.1", "4.7.1", "4.19.1")
# another mouse tissue
## 4.19.1. 6.2.1 4.13.4. 9.1.3(Skin) 1.4.1. 4.7.1 (blood, kidney) 4.7.2
sname = c("Mouse","Horse","Human",
"Shrew","Mouse","Naked mole-rat",
"Green monkey","Pig", "Beluga whale")
tissues = c("Blood", "Blood", "Skin",
"Tail","Cerebellum","Skin",
"Cortex","Blood", "Blood")
data.frame(
ss,sname,tissues
)
```
```{r}
ss_age = list()
for (nk in 1:length(ss)) {
ylab1 = ifelse(nk==1,"ScaledM","")
idxss = with(dat_scaledM_rage, which(SpecNum == ss[nk]& Tissue == tissues[nk]))
cor1 = round(cor(dat_scaledM_rage$ScaledM[idxss], dat_scaledM_rage$Age[idxss]),2)
tit1 = paste0(letters[nk+ceiling(nk/3)*3],". ",sname[nk],"_",tissues[nk],", Cor=",cor1)
xint = 0.1*dat_scaledM_rage$Lifepspan[idxss[1]]
datplot = dat_scaledM_rage[idxss,]
# ymin = min(datplot$MeanMethyl)
# ymax = max(datplot$MeanMethyl)
ymin = min(datplot$ScaledM)
ymax = max(datplot$ScaledM)
maxage = dat_scaledM_rage$Lifepspan[idxss[1]]
maturity = dat_scaledM_rage$ASM[idxss[1]]
lm1 = lm(ScaledM ~ Age, data = datplot[datplot$Age<= xint,])
xr1 = range(datplot[datplot$Age<= xint,]$Age)
yr1 = range(lm1$fitted.values)
slopey = coef(lm1)[2]
lm2 = lm(ScaledM ~ Age, data = datplot[datplot$Age> xint,])
xr2 = range(datplot[datplot$Age> xint,]$Age)
yr2 = range(lm2$fitted.values)
slopeo = coef(lm2)[2]
lm3 = lm(ScaledM ~ Age, data = datplot)
xr3 = range(datplot$Age)
yr3 = range(lm3$fitted.values)
slopea = coef(lm3)[2]
segdat = data.frame(
x = c(xr1[1], xr2[1], xr3[1]), ## 0, xint, 0,
y = c( yr1[1], yr2[1], yr3[1]), ## ymin*0.92, ymin*0.91, ymin*0.9,
xend = c( xr1[2], xr2[2], xr3[2]), ## xint, 10*xint, 10*xint,
yend = c( yr1[2], yr2[2], yr3[2]), ## ymin*0.92, ymin*0.91, ymin*0.9,
col = c("blue","red","black"),
slopex = c(mean(xr1), mean(xr2), mean(xr3)),
slopey = c(yr1[1], yr2[1], yr3[2]),
slopet = round(c(slopey,slopeo,slopea),2)
)
ss_age[[nk]] = ggplot(datplot,
aes(x=Age, y = ScaledM)) +
geom_point(size=1, alpha=0.5) + geom_vline(xintercept = xint,linetype=2)+
## geom_text(hjust=0, vjust=0, size = 1) +
# geom_smooth(method = lm,se = FALSE,linewidth=0.5,linetype=2)+
xlab("Age") + ylab(ylab1) + ggtitle(tit1) + # ylim(ymin, ymax) +
geom_segment(aes(x=x,y=y,xend=xend,yend=yend,col=col),data=segdat) +
geom_text(aes(x=slopex, y=slopey, label = slopet, col=col), # color=col,
# label.size = 0, fill=NA,
fontface = "bold",
data=segdat) +
theme_bw() + theme(legend.position = "",plot.title = element_text(size = 10))
}
ss_age[[nk]]
```
- last 3 species as panels d,e,f
```{r}
for (nk in 7:9 ) {
ylab1 = ifelse(nk==7,"ScaledM","")
idxss = with(dat_scaledM_rage, which(SpecNum == ss[nk]& Tissue == tissues[nk]))
cor1test = cor.test(dat_scaledM_rage$ScaledM[idxss], dat_scaledM_rage$Age[idxss])
cor1 = round(cor1test$estimate,2)
pv1 = signif(cor1test$p.value,1)
tit1 = paste0(letters[nk-3],". ",
sname[nk],"_",tissues[nk],
", Cor=",cor1,", p=",pv1)
xint = 0.1*dat_scaledM_rage$Lifepspan[idxss[1]]
datplot = dat_scaledM_rage[idxss,]
# ymin = min(datplot$MeanMethyl)
# ymax = max(datplot$MeanMethyl)
ymin = min(datplot$ScaledM)
ymax = max(datplot$ScaledM)
maxage = dat_scaledM_rage$Lifepspan[idxss[1]]
maturity = dat_scaledM_rage$ASM[idxss[1]]
lm1 = lm(ScaledM ~ Age, data = datplot[datplot$Age<= xint,])
xr1 = range(datplot[datplot$Age<= xint,]$Age)
yr1 = range(lm1$fitted.values)
slopey = coef(lm1)[2]
lm2 = lm(ScaledM ~ Age, data = datplot[datplot$Age> xint,])
xr2 = range(datplot[datplot$Age> xint,]$Age)
yr2 = range(lm2$fitted.values)
slopeo = coef(lm2)[2]
lm3 = lm(ScaledM ~ Age, data = datplot)
xr3 = range(datplot$Age)
yr3 = range(lm3$fitted.values)
slopea = coef(lm3)[2]
segdat = data.frame(
x = c(xr1[1], xr2[1], xr3[1]), ## 0, xint, 0,
y = c( yr1[1], yr2[1], yr3[1]), ## ymin*0.92, ymin*0.91, ymin*0.9,
xend = c( xr1[2], xr2[2], xr3[2]), ## xint, 10*xint, 10*xint,
yend = c( yr1[2], yr2[2], yr3[2]), ## ymin*0.92, ymin*0.91, ymin*0.9,
col = c("blue","red","black"),
slopex = c(mean(xr1), mean(xr2), mean(xr3)),
slopey = c(yr1[1], yr2[1], yr3[2]),
slopet = round(c(slopey,slopeo,slopea),2)
)
ss_age[[nk]] = ggplot(datplot,
aes(x=Age, y = ScaledM)) +
geom_point(size=1, alpha=0.5) + geom_vline(xintercept = xint,linetype=2)+
## geom_text(hjust=0, vjust=0, size = 1) +
# geom_smooth(method = lm,se = FALSE,linewidth=0.5,linetype=2)+
xlab(paste0("Age (N=",length(idxss),")" )) + ylab(ylab1) + ggtitle(tit1) + # ylim(ymin, ymax) +
geom_segment(aes(x=x,y=y,xend=xend,yend=yend,col=col),data=segdat) +
geom_text(aes(x=slopex, y=slopey, label = slopet, col=col), # color=col,
# label.size = 0, fill=NA,
fontface = "bold",
data=segdat) +
theme_bw() + theme(legend.position = "",plot.title = element_text(size = 10))
}
```
```{r}
vnum = "v13"
pdf(paste0("../",outfolder,"/Fig1_",vnum,".pdf"),
width = 10.5,height = 7)
grid.arrange(fig1a, p_ra,p_logli,
ss_age[[7]],ss_age[[8]],ss_age[[9]],
# sslist[[1]],sslist[[2]],sslist[[3]],
nrow=2)
dev.off()
```
#### ss_age: relative age and logliage
```{r}
nkk = c(1:3,1:3,7:9)
for (nk in 4:length(ss)) {
ylab1 = ifelse(nk%%3==1,"ScaledM","")
idxss = with(dat_scaledM_rage, which(SpecNum == ss[nk]& Tissue == tissues[nk]))
cor1test = cor.test(dat_scaledM_rage$ScaledM[idxss], dat_scaledM_rage$Age[idxss])
cor1 = round(cor1test$estimate,2)
pv1 = signif(cor1test$p.value,1)
tit1 = paste0(letters[nkk[nk]],". ",sname[nk],"_",tissues[nk],
", Cor=",cor1, ", p=",pv1)
xint = 0.1
datplot = dat_scaledM_rage[idxss,]
# ymin = min(datplot$MeanMethyl)
# ymax = max(datplot$MeanMethyl)
ymin = min(datplot$ScaledM)
ymax = max(datplot$ScaledM)
maxage = dat_scaledM_rage$Lifepspan[idxss[1]]
maturity = dat_scaledM_rage$ASM[idxss[1]]
lm0 = lm(ScaledM ~ logliage, data = datplot)
datplot$fitted_logli = lm0$fitted.values
lm3 = lm(ScaledM ~ Rage, data = datplot)
datplot$fitted_rage = lm3$fitted.values
segdat = data.frame(
x = c(0, xint, 0),
y = c(ymin*0.92, ymin*0.91, ymin*0.9),
xend = c(xint, 10*xint, 10*xint),
yend = c(ymin*0.92, ymin*0.91, ymin*0.9),
col = c("blue","red","black"),
slopex = c(0.1, 0.7, 0.5),
slopey = c(ymin*0.97, ymin*0.95, ymin*0.85),
slopet = round(c(slopey,slopeo,slopea),2)
)
ss_age[[nk]] = ggplot(datplot,
aes(x=Rage, y = ScaledM)) +
geom_point(size=1, alpha=0.5) + geom_vline(xintercept = xint,linetype=2)+
## geom_text(hjust=0, vjust=0, size = 1) +
# geom_smooth(method = lm,se = FALSE,linewidth=0.5,linetype=2)+
geom_line(aes(x = Rage, y = fitted_rage), # se = FALSE,
linewidth=0.5,linetype=2,col="red") +
geom_smooth(aes(x = Rage, y = fitted_logli), se = FALSE,
linewidth=0.5,linetype=2,col="blue") +
xlab(paste0("Relative Age (N=",length(idxss),")" )) +
ylab(ylab1) + ggtitle(tit1) + # ylim(ymin, ymax) +
# geom_segment(aes(x=x,y=y,xend=xend,yend=yend,col=col),data=segdat) +
# geom_label(aes(x=slopex, y=slopey, col=col, label = slopet),data=segdat,
# label.size = NA, fill=NA) +
theme_bw() + theme(legend.position = "",plot.title = element_text(size = 10))
}
ss_age[[nk]]
```
#### sslist
```{r}
sslist = list()
nkk = c(1:6, 10:12)
for (nk in 1:length(ss)) {
ylab1 = ifelse(nk%%3==1,"ScaledM","")
idxss = with(dat_scaledM_rage, which(SpecNum == ss[nk]& Tissue == tissues[nk]))
cor1test = cor.test(dat_scaledM_rage$ScaledM[idxss], dat_scaledM_rage$Age[idxss])
cor1 = round(cor1test$estimate,2)
pv1 = signif(cor1test$p.value,1)
tit1 = paste0(letters[nkk[nk]],". ",sname[nk],"_",tissues[nk],
", Cor=",cor1, ", p=",pv1)
datplot = dat_scaledM_rage%>%
filter(SpecNum == ss[nk]&Tissue == tissues[nk])
lm0 = lm(ScaledM ~ logliage, data = datplot)
gamma1 = round(coef(lm0)[2],2)
pv2 = signif(summary(lm0)$coefficients[2,4],1)
sslist[[nk]] = ggplot(datplot,
aes(x=logliage, y = ScaledM)) +
geom_point(size=1, alpha=0.5) +
## geom_text(hjust=0, vjust=0, size = 1) +
geom_smooth(method = lm,se = FALSE,linewidth=0.5,linetype=2)+
# geom_abline(slope = coef(lm_tage)[2], intercept = coef(lm_tage)[1]) +
xlab(paste0("Logli Relative Age (N=",length(idxss),")")) +
ylab(ylab1) +
theme_bw() +
labs(title = tit1,
subtitle = bquote(gamma[1] ~ "=" ~ .(gamma1) ~ ", p="~ .(pv2)) ) +
theme(plot.title=element_text(size = titlesize),
legend.position = "",
plot.subtitle = element_text(size = 1.5*titlesize,
margin=margin(b=-20)))
}
sslist[[nk]]
```
#### Fig1_ctd
```{r}
pdf(paste0("../",outfolder,"/Fig1_ctd_v4.pdf"),
width = 9.6,height = 12)
grid.arrange(ss_age[[4]],ss_age[[5]],ss_age[[6]],
sslist[[4]],sslist[[5]],sslist[[6]],
ss_age[[7]],ss_age[[8]],ss_age[[9]],
sslist[[7]],sslist[[8]],sslist[[9]],
nrow=4)
dev.off()
```
```{r}
pdf(paste0("../",outfolder,"/SingleSpecies_1.pdf"),
width = 9,height = 6)
grid.arrange(ss_age[[1]],ss_age[[2]],ss_age[[3]],
sslist[[1]],sslist[[2]],sslist[[3]],
nrow=2)
dev.off()
pdf(paste0("../",outfolder,"/SingleSpecies_2.pdf"),
width = 9,height = 6)
grid.arrange(ss_age[[4]],ss_age[[5]],ss_age[[6]],
sslist[[4]],sslist[[5]],sslist[[6]],
nrow=2)
dev.off()
pdf(paste0("../",outfolder,"/SingleSpecies_3.pdf"),
width = 9,height = 6)
grid.arrange(ss_age[[7]],ss_age[[8]],ss_age[[9]],
sslist[[7]],sslist[[8]],sslist[[9]],
nrow=2)
dev.off()
```
```{r}
vnum = "v8"
pdf(paste0("../",outfolder,"/Fig1_",vnum,".pdf"),
width = 9,height = 9)
grid.arrange(fig1a, p_ra,p_logli,
ss_age[[1]],ss_age[[2]],ss_age[[3]],
sslist[[1]],sslist[[2]],sslist[[3]],
nrow=3)
dev.off()
```
```{r}
vnum = "v9"
pdf(paste0("../",outfolder,"/Fig1_",vnum,".pdf"),
width = 9,height = 9)
grid.arrange(fig1a, p_rage,p_logliage,
ss_age[[1]],ss_age[[2]],ss_age[[3]],
sslist[[1]],sslist[[2]],sslist[[3]],
nrow=3)
dev.off()
```
## Figure 5: Young and Old AROCM
```{r}
sd_ragemat = ragemat%>%select(contains("sd_rage"))
colnames(sd_ragemat)
k=1
l=11
sdrage = sd_ragemat[,l]
```
```{r}
setwd("/Users/feiz/Documents/Box Sync/Horvath/Mammalian meth projects/Slope")
textnum = c("1.1.1","4.11.1","9.9.1","9.9.3","12.3.10",
"5.1.1", "1.4.1","6.1.1", "2.1.1", "8.4.5" )
statesfolder = "out_AROCM_0.1L"
kk=1
plist = list()
titlesize = 10
for (j in c(2,1,3)) {
len1 = stateSlopes$statesPRC2[j]
# slopes_mat = SpeciesMat%>%select(contains(len1),maxAgeCaesar,SpeciesLatinName)
matslopes = read.csv(paste0(statesfolder,"/states/Slopes_IdentityScaled",len1, "v11.csv"), row.names = 1)
slopes_mat = matslopes%>%select(contains("Slope"),
maxAgeCaesar,SpeciesLatinName)
source("slope_strata_v2.R")
##### species level
idx1 = which(SpeciesMat$RemoveFei==0 &
SpeciesMat$ageRangeL< matratio &
!is.na(sdrage) & sdrage>= 0.06 )
idx1 = which(!is.na(slopes_mat[,k]) & slopes_mat[,k]<10 &
slopes_mat[,k] > -0.9 &
SpeciesMat$RemoveFei==0 &
# SpeciesMat$ageRangeL< matratio &
# SpeciesMat$ageRangeU >= min(prop,0.6) &
!is.na(sdrage) & sdrage> 0.03 )
# colnames(slopes_mat)
nc = 1:14
specmat = combineTissues(slopes_mat, idx1, nc)
# colnames(specmat)
prop1 = props[k]
prop2 = oldprops[l-10]
tit1 = paste0("(0,",prop1,"*Lifespan)")
tit2 = paste0("(0.1*Lifespan,Lifespan)")
x = specmat[,k+1]
y = specmat[,l+1]
### Log scale
min1 = min(c(x,y), na.rm = TRUE)
min1 = ifelse(min1<0, 2*min1, 0)
print(min1)
x = log(x-min1)
y = log(y-min1)
xr = range(c(x,y), na.rm = TRUE)+c(-0.5,2)
# specmat$textidx = (x>2 | x< -2 | y>0 | x-y>3 | x-y<0)
specmat$textidx = specmat$MammalNumberHorvath %in% textnum
specmat$x = x
specmat$y = y
cor1test = with(specmat,cor.test(x,y,use = "pairwise.complete"))
cor1 = round(cor1test$estimate, 3)
pv1 = signif(cor1test$p.value,1)
coef1 = coef(lm(y~x, data = specmat))
mytitle = paste0(letters[kk],". ",len1,', N=',
length(y))
p1 = ggplot(specmat,
aes_string(y = "y", x = "x",
label = "MammalNumberHorvath")) +
geom_point(color="darkgrey") + theme_bw() +
scale_x_continuous(name='Log Young AROCM', limits = xr) + ## , limits = xlim
scale_y_continuous(name='Log Old AROCM', limits = xr) +
labs(title = mytitle, subtitle =
paste0("Cor=",cor1, ", p=",pv1 )) +
theme(plot.title=element_text(size = titlesize),
plot.subtitle = element_text(size = 1.5*titlesize,
margin=margin(b=-20))) + # margin=margin(t=40,b=-30)
geom_text(aes(label=ifelse(textidx,
MammalNumberHorvath,'')),
hjust=0,vjust=0) +
geom_abline(slope=1, intercept = 0, linetype=3)
plist[[kk]] = p1
kk=kk+1
}
kk
length(plist)
plist[[3]]
```
```{r}
vnum = "v10"
ggsave(paste0(outfolder,"/Figure2_",vnum, "_species.png"),
grid.arrange(plist[[1]],plist[[2]],plist[[3]],
nrow=1),
width = 24,
height = 8, units = "cm")
```
### Figure 2 original scale
```{r}
# setwd("/Users/feiz/Documents/Box Sync/Horvath/Mammalian meth projects/Slope")
textnum = c("1.1.1","4.11.1","9.9.1","9.9.3","12.3.10",
"5.1.1", "1.4.1","6.1.1", "2.1.1", "8.4.5" )
# statesfolder = "out_AROCM"
titlesize = 10
for (j in c(2,1,3)) {
len1 = stateSlopes$statesPRC2[j]
# slopes_mat = SpeciesMat%>%select(contains(len1),maxAgeCaesar,SpeciesLatinName)
matslopes = read.csv(paste0("../",statesfolder,"/states/Slopes_IdentityScaled",len1, "v11.csv"), row.names = 1)
slopes_mat = matslopes%>%select(contains("Slope"),
maxAgeCaesar,SpeciesLatinName)
source("../slope_strata_v2.R")
##### species level
idx1 = which(SpeciesMat$RemoveFei==0 &
SpeciesMat$ageRangeL< matratio &
!is.na(sdrage) & sdrage>= 0.06 )
idx1 = which(!is.na(slopes_mat[,k]) & slopes_mat[,k]<10 &
slopes_mat[,k] > -0.9 &
SpeciesMat$RemoveFei==0 &
# SpeciesMat$ageRangeL< matratio &
# SpeciesMat$ageRangeU >= min(prop,0.6) &
!is.na(sdrage) & sdrage> 0.03 )
# colnames(slopes_mat)
nc = 1:14
specmat = combineTissues(slopes_mat, idx1, nc)
# colnames(specmat)
prop1 = props[k]
prop2 = oldprops[l-10]
tit1 = paste0("(0,",prop1,"*Lifespan)")
tit2 = paste0("(0.1*Lifespan,Lifespan)")
x = specmat[,k+1]
y = specmat[,l+1]
### Log scale
min1 = min(c(x,y), na.rm = TRUE)
min1 = ifelse(min1<0, 2*min1, 0)
print(min1)
# x = log(x-min1)
# y = log(y-min1)
xr = range(c(x,y), na.rm = TRUE)+c(-0.5,2)
# specmat$textidx = (x>2 | x< -2 | y>0 | x-y>3 | x-y<0)
specmat$textidx = specmat$MammalNumberHorvath %in% textnum
specmat$x = x
specmat$y = y
cor1test = with(specmat,cor.test(x,y,use = "pairwise.complete"))
cor1 = round(cor1test$estimate, 3)
pv1 = signif(cor1test$p.value,1)
coef1 = coef(lm(y~x, data = specmat))
mytitle = paste0(letters[kk],". ",len1,', N=',
length(y)) ## ,", Cor=",cor1) letters[kk],". ",
# xr = range(c(specmat$x, specmat$y))*c(1/2, 5)
p1 = ggplot(specmat,
aes_string(y = "y", x = "x",
label = "MammalNumberHorvath")) +
geom_point(color="darkgrey") + theme_bw() +
scale_x_continuous(name='Young AROCM', limits = xr) + ## , limits = xlim
scale_y_continuous(name='Old AROCM', limits = xr) +
labs(title = mytitle, subtitle =
paste0("Cor=",cor1, ", p=",pv1 )) +
theme(plot.title=element_text(size = titlesize),
plot.subtitle = element_text(size = 1.5*titlesize,
margin=margin(b=-20))) + # margin=margin(t=40,b=-30)
geom_text(aes(label=ifelse(textidx,
MammalNumberHorvath,'')),
hjust=0,vjust=0) +
geom_abline(slope=1, intercept = 0, linetype=3)
plist[[kk]] = p1
kk=kk+1
}
kk
length(plist)
plist[[6]]
```
```{r}
vnum = "v11"
ggsave(paste0("../",outfolder,"/Figure2_",vnum, "_species.png"),
grid.arrange(plist[[1]],plist[[2]],plist[[3]],
plist[[4]],plist[[5]],plist[[6]],
nrow=2),
width = 24,
height = 16, units = "cm")
```
## Figure 6: Adj.AROCM
```{r}
# SpeciesMat = read_csv(paste0(outfolder,"/ALL_states_IdentityScaled_SpeciesTissue_13Slopes_v10.csv"))
SpeciesMat = SpeciesMat[,-1]
dim(SpeciesMat)
SpeciesMat = SpeciesMat%>%
arrange(maxAgeCaesar)
write.csv(SpeciesMat, "SupplementTable1.csv")
```
```{r, eval=FALSE}
infolder = "out_AROCM_0.1L"
SpeciesMat = read.csv("/Users/feiz/Documents/Box Sync/Horvath/Mammalian meth projects/Slope/out_AROCM_0.1L/datSamp_BivProm2+_Slopes_logliAge.csv", row.names = 1, check.names = FALSE)
head(colnames(SpeciesMat),22)
```
```{r}
j=2
len1 = stateSlopes$statesPRC2[j]
# source("slope_strata_v2.R")
# dim(slopes_mat)
slopes_mat = SpeciesMat%>%select(contains(len1),
maxAgeCaesar,SpeciesLatinName,Tissue)
colnames(slopes_mat)
slopes_mat = slopes_mat[,-(1:5)]
source("../slope_strata_v2.R")
```
### ragemat
```{r}
ragemat = matslopes%>%select(contains("rage"))
sd_ragemat = ragemat%>%select(contains("sd_rage"))
colnames(sd_ragemat)
# k=1
l=11
sdrage = sd_ragemat[,l]
```
### Main plot (nj=12)
```{r}
## outfolder = "Fig6_AROCM"
colnames(SpeciesMat)[1:16]
colnames(SpeciesMat)[10:11] = c("GT", "ASM")
dim(SpeciesMat)
dim(slopes_mat)
colnames(slopes_mat)
identical(SpeciesMat$SpeciesLatinName, slopes_mat$SpeciesLatinName)
x = SpeciesMat$maxAgeCaesar
y = SpeciesMat$`IdentityScaledBivProm2+YoungSlope1`
plot(log(x), log(y))
```
```{r}
nj = 12 ##Lifespan 10:GT 11:ASM
colnames(SpeciesMat)[nj] = "Lifespan"
colnames(slopes_mat)[1:13] =
sapply(strsplit(colnames(slopes_mat)[1:13],"+",fixed = TRUE), function(v1){v1[2]})
colnames(slopes_mat)
slopes_mat[, 14] = SpeciesMat[,nj]
colnames(slopes_mat)[14] = colnames(SpeciesMat)[nj]
```
#### Including outliers
```{r}
idx1 = which(SpeciesMat$Freq>= 3 &
SpeciesMat$ageRangeL != SpeciesMat$ageRangeU)
colnames(slopes_mat)
ncs = 1:14
speciesSlopes_all = combineTissues(slopes_mat, idx1, ncs)
dim(speciesSlopes_all)
p_out = list()
colnames(speciesSlopes_all)
```
```{r}
speciesSlopes_all$x = x = speciesSlopes_all$Lifespan
speciesSlopes_all$y = y = speciesSlopes_all[,11]
plot(x,y)
plot(x,1/y)
plot(log(x), log(y))
```
- a. AROCM vs. Lifespan
```{r}
textnum = c("1.1.1","4.11.1","9.9.1","12.3.10", "9.9.15", "9.5.5",
"1.7.1", "1.4.1","11.1.3", "2.1.1", "8.4.5" )
speciesSlopes_all$x = x = speciesSlopes_all$Lifespan
speciesSlopes_all$y = y = speciesSlopes_all[,11]
xlim = range(x)+c(0,20)
cor1test = with(speciesSlopes_all,cor.test(x,y))
cor1 = round(cor1test$estimate, 2)
pv1 = signif(cor1test$p.value, 1)
# coef1 = coef(lm(y~x, data = speciesSlopes_all))
speciesSlopes_all$textidx =
speciesSlopes_all$MammalNumberHorvath %in% textnum
mytitle = paste0('a. N=', length(y),
", Cor=",cor1,", p=",pv1)
p_out[[1]] = ggplot(speciesSlopes_all,
aes_string(y = "y", x = "x",
label = "MammalNumberHorvath")) +
geom_point(color="darkgrey") + theme_bw() +
scale_x_continuous(name='Lifespan', limits = xlim)+
scale_y_continuous(name=paste0('AROCM'))+
ggtitle(mytitle) +
theme(plot.title=element_text(size = titlesize)) +
geom_text(aes(label=ifelse(textidx,
MammalNumberHorvath,'')),
hjust=0,vjust=0)
p_out[[1]]
```
- b. log AROCM vs. log Lifespan
```{r}
plotslopes = speciesSlopes_all%>%filter(y>0)
plotslopes$x = x = log(plotslopes$Lifespan)
plotslopes$y = y = log(plotslopes[,11])
xlim = c(0,6)
cor1test = with(plotslopes,cor.test(x,y))
cor1 = round(cor1test$estimate, 2)
pv1 = signif(cor1test$p.value, 1)
coef1 = coef(lm(y~x, data = plotslopes))
plotslopes$textidx =
plotslopes$MammalNumberHorvath %in% textnum
mytitle = paste0('b. N=', length(y),
", Cor=",cor1,", p=", pv1)
p_out[[2]] = ggplot(plotslopes,
aes_string(y = "y", x = "x",
label = "MammalNumberHorvath")) +
geom_point(color="darkgrey") + theme_bw() +
scale_x_continuous(name='Log Lifespan', limits = xlim)+
scale_y_continuous(name=paste0('Log AROCM'))+
ggtitle(mytitle) +
theme(plot.title=element_text(size = titlesize)) +
geom_abline(slope=coef1[2], intercept = coef1[1], linetype=3) +
geom_text(aes(label=ifelse(textidx,
MammalNumberHorvath,'')),
hjust=0,vjust=0)
p_out[[2]]
```
- c. 1/AROCM vs. Lifespan
```{r}
plotslopes = speciesSlopes_all
plotslopes$x = x = 1/ plotslopes$Lifespan
plotslopes$y = y = plotslopes[,11]
xlim = range(x)*c(1/2, 1.4)
cor1test = with(plotslopes,cor.test(x,y))
cor1 = round(cor1test$estimate, 2)
pv1 = signif(cor1test$p.value, 1)
coef1 = coef(lm(y~x, data = plotslopes))
textnum = c("1.1.1","4.11.1","9.9.1","12.3.10", "9.9.15", "9.5.5",
"1.7.1", "1.4.1","11.1.3", "2.1.1", "8.4.5" )
plotslopes$textidx =
plotslopes$MammalNumberHorvath %in% textnum
mytitle = paste0('c. N=', length(y),
", Cor=",cor1,", p=", pv1)
p_out[[3]] = ggplot(plotslopes,
aes_string(y = "y", x = "x",
label = "MammalNumberHorvath")) +
geom_point(color="darkgrey") + theme_bw() +
scale_x_continuous(name='1/Lifespan', limits = xlim)+
scale_y_continuous(name=paste0('AROCM'))+
ggtitle(mytitle) +
theme(plot.title=element_text(size = titlesize)) +
geom_abline(slope=coef1[2], intercept = coef1[1], linetype=3) +
geom_text(aes(label=ifelse(textidx,
MammalNumberHorvath,'')),
hjust=0,vjust=0)
p_out[[3]]
```
#### Species-tissue level
- a. AROCM vs. Lifespan
```{r}
colnames(slopes_mat)
plotslopes = slopes_mat%>%filter(!is.na(YoungSlope1))
tmp1 = match(plotslopes$SpeciesLatinName, SpeciesMat$SpeciesLatinName)
plotslopes$MammalNumberHorvath = SpeciesMat$MammalNumberHorvath[tmp1]
plotslopes$x = x = plotslopes$Lifespan
plotslopes$y = y = plotslopes[,10]
xlim = range(x) + c(0,20)
cor1test = with(plotslopes,cor.test(x,y))
cor1 = round(cor1test$estimate, 2)
pv1 = signif(cor1test$p.value, 1)
plotslopes$textidx =
plotslopes$MammalNumberHorvath %in% textnum
mytitle = paste0('d. N=', length(y),
", Cor=",cor1,", p=",pv1)
p_out[[4]] = ggplot(plotslopes,
aes_string(y = "y", x = "x",col="Tissue",
label = "MammalNumberHorvath")) +
geom_point(alpha=0.5) + theme_bw() +
scale_x_continuous(name='Lifespan', limits = xlim)+
scale_y_continuous(name=paste0('AROCM'))+
ggtitle(mytitle) +
theme(legend.position = "",
plot.title=element_text(size = titlesize)) +
geom_text(aes(label=ifelse(textidx,
MammalNumberHorvath,'')),
hjust=0,vjust=0)
p_out[[4]]
```
- b. log AROCM vs. log Lifespan
```{r}
plotslopes = slopes_mat%>%filter(!is.na(YoungSlope1) &
YoungSlope1 >0)
tmp1 = match(plotslopes$SpeciesLatinName, SpeciesMat$SpeciesLatinName)
plotslopes$MammalNumberHorvath = SpeciesMat$MammalNumberHorvath[tmp1]
plotslopes$x = x = log(plotslopes$Lifespan)
plotslopes$y = y = log(plotslopes[,10])
xlim = c(0,6)
cor1test = with(plotslopes,cor.test(x,y))
cor1 = round(cor1test$estimate, 2)
pv1 = signif(cor1test$p.value, 1)
coef1 = coef(lm(y~x, data = plotslopes))
textnum = c("1.1.1","4.11.1","9.9.1","12.3.10", ## "9.9.3",
"5.1.1", "1.4.3","6.1.1", "2.1.1", "8.4.5" )
plotslopes$textidx =
plotslopes$MammalNumberHorvath %in% textnum
mytitle = paste0('e. N=', length(y),
", Cor=",cor1,", p=", pv1)
p_out[[5]] = ggplot(plotslopes,
aes_string(y = "y", x = "x", col="Tissue",
label = "MammalNumberHorvath")) +
geom_point(alpha = 0.5) + theme_bw() +
scale_x_continuous(name='Log Lifespan', limits = xlim)+
scale_y_continuous(name=paste0('Log AROCM'))+
ggtitle(mytitle) +
theme(legend.position = "",
plot.title=element_text(size = titlesize)) +
geom_abline(slope=coef1[2], intercept = coef1[1], linetype=3) +
geom_text(aes(label=ifelse(textidx,
MammalNumberHorvath,'')),
hjust=0,vjust=0)
p_out[[5]]
```
- c. 1/AROCM vs. Lifespan
```{r}
plotslopes = slopes_mat%>%filter(!is.na(YoungSlope1) )
tmp1 = match(plotslopes$SpeciesLatinName, SpeciesMat$SpeciesLatinName)
plotslopes$MammalNumberHorvath = SpeciesMat$MammalNumberHorvath[tmp1]
plotslopes$x = x = 1/plotslopes$Lifespan
plotslopes$y = y = plotslopes[,10]
xlim = range(x)*c(1/2, 1.4)
cor1test = with(plotslopes,cor.test(x,y))
cor1 = round(cor1test$estimate, 2)
pv1 = signif(cor1test$p.value, 1)
coef1 = coef(lm(y~x, data = plotslopes))
textnum = c("1.1.1","4.11.1","9.9.1","12.3.10", "9.9.15", "9.5.5",
"1.7.1", "1.4.1","10.1.3", "2.1.1", "6.1.1" )
plotslopes$textidx =
plotslopes$MammalNumberHorvath %in% textnum
mytitle = paste0('f. N=', length(y),
", Cor=",cor1, ", p=", pv1)
p_out[[6]] = ggplot(plotslopes,
aes_string(y = "y", x = "x",col="Tissue",
label = "MammalNumberHorvath")) +
geom_point(alpha = 0.5) + theme_bw() +
scale_x_continuous(name='1/Lifespan', limits = xlim)+
scale_y_continuous(name=paste0('AROCM'))+
ggtitle(mytitle) +
theme(legend.position = "",
plot.title=element_text(size = titlesize)) +
geom_abline(slope=coef1[2], intercept = coef1[1], linetype=3) +
geom_text(aes(label=ifelse(textidx,
MammalNumberHorvath,'')),
hjust=0,vjust=0)
p_out[[6]]
```
- combined plot
```{r}
pdf(paste0("SuppFig_AROCM_outliers","_v2.pdf"),
width = 11.1,height = 7.4)
grid.arrange(p_out[[1]],p_out[[2]],p_out[[3]],
p_out[[4]],p_out[[5]],p_out[[6]],
nrow=2)
dev.off()
```
#### Excluding outliers
```{r}
###### new text
textnum = c("1.1.1","4.11.1","9.9.1","12.3.10", ## "9.9.3",
"5.1.1", "1.4.1","6.1.1", "2.1.1", "8.4.5" )
pw = 0.75
SpeciesMat$RemoveFei[SpeciesMat$MammalNumberHorvath %in%
c("6.1.2")] = 1
SpeciesMat$RemoveFei[SpeciesMat$MammalNumberHorvath %in%
c("4.13.3") & SpeciesMat$Tissue=="Skin"] = 1
```
```{r}
kk=1
pAROCM = list()
titlesize = 10
# for(k in c(10,1,11)){
for(k in c(10)){
sname = ifelse(k==10, "AROCM",
ifelse(k==11, "Old AROCM", "Young AROCM"))
prop = ifelse(k<= 10, props[k], oldprops[k-10])
tit1 = ifelse(k<= 10, paste0(", (0, ",props[k],"*Lifespan)"),
paste0(", (Maturity, Lifespan)"))
if(k<= 10){
idx1 = which(!is.na(slopes_mat[,k]) & slopes_mat[,k]<10 &
slopes_mat[,k] > -0.9 &
SpeciesMat$RemoveFei==0
# SpeciesMat$ageRangeL< matratio &
# & SpeciesMat$ageRangeU >= min(prop,0.6)
)
}else {
sdrage = sd_ragemat[,k]
idx1 = which(!is.na(slopes_mat[,k]) & slopes_mat[,k]<10 &
slopes_mat[,k] > -0.1 &
SpeciesMat$RemoveFei==0 &
!is.na(sdrage) & sdrage>= 0.06 )
}
ncs = 1:14 ###c(k, grep("maxAgeCaesar",colnames(slopes_mat)))
speciesSlopes = combineTissues(slopes_mat, idx1, ncs)
x = speciesSlopes[,15]
y = speciesSlopes[,1+k]
y = logplus(y)
speciesSlopes$y = log(y)
speciesSlopes$x = log(x)
xlim = range(speciesSlopes$x)*c(0.7, 1.2)
cor1test = with(speciesSlopes,cor.test(x,y))
cor1 = round(cor1test$estimate, 2)
pv1 = signif(cor1test$p.value,1)
coef1 = coef(lm(y~x, data = speciesSlopes))
speciesSlopes$textidx =
speciesSlopes$MammalNumberHorvath %in% textnum
mytitle = paste0(letters[kk],". ",sname,', N=',
length(y),", Cor=",cor1,", p=", pv1)
p1 = ggplot(speciesSlopes,
aes_string(y = "y", x = "x",
label = "MammalNumberHorvath")) +
geom_point(color="darkgrey") + theme_bw() +
scale_x_continuous(name='Log Lifespan', limits = xlim)+
scale_y_continuous(name=paste0('Log ',sname))+
ggtitle(mytitle) +
theme(plot.title=element_text(size = titlesize)) + # margin=margin(t=40,b=-30)
geom_text(aes(label=ifelse(textidx,
MammalNumberHorvath,'')),
hjust=0,vjust=0) +
geom_abline(slope=coef1[2], intercept = coef1[1], linetype=3)
pAROCM[[kk]] = p1
kk=kk+1
##### panel b
{
b=1
tmpx = (x)^(-b)
lm1 = lm(y~ tmpx-1, weights = tmpx)
coef1 = as.numeric(coef(lm1))
x1 = coef1[1]*(x)^(-b)
a = round(coef1[1],2)
speciesSlopes$x = a/x
speciesSlopes$y = y
cor1test = with(speciesSlopes,cor.test(x,y))
cor1 = round(cor1test$estimate, 2)
pv1 = signif(cor1test$p.value,1)
coef1 = coef(lm(y~x, data = speciesSlopes))
mytitle = paste0(letters[kk],". ",sname,', N=',
length(y)) ##letters[kk],". ",
xr = range(c(speciesSlopes$x, speciesSlopes$y)) *c(1/2, 5)
p1 = ggplot(speciesSlopes,
aes_string(y = "y", x = "x",
label = "MammalNumberHorvath")) +
geom_point(color="darkgrey") + theme_bw() +
scale_x_log10(name=paste0("a/",colnames(speciesSlopes)[15],
", a = ",a), limits = xr) + ## , limits = xlim
scale_y_log10(name=paste0(sname), limits = xr) +
labs(title = mytitle, subtitle =
paste0("Cor=",cor1, ", p=",pv1 )) +
theme(plot.title=element_text(size = titlesize),
plot.subtitle = element_text(size = 1.5*titlesize,
margin=margin(b=-20))) + # margin = margin(t=-5)
geom_text(aes(label=ifelse(textidx,
MammalNumberHorvath,'')),
hjust=0,vjust=0) +
geom_abline(slope=coef1[2], intercept = coef1[1], linetype=3)
pAROCM[[kk]] = p1
kk=kk+1
}
#### panel c & d
{
xmat = cbind(slopes_mat[,1:13]*sd_ragemat^pw,
1/SpeciesMat[,nj], slopes_mat[,14:15])
colnames(xmat)[14] = col2 = paste0("1/", colnames(SpeciesMat)[nj])
ncs = 1:15 ###c(k, grep("maxAgeCaesar",colnames(slopes_mat)))
species_xmat = combineTissues(xmat, idx1, ncs)
species_xmat$textidx = species_xmat$MammalNumberHorvath %in% textnum
## log scale panel
x = species_xmat[, 16]
y = species_xmat[,k+1]
y = logplus(y)
species_xmat$y = log(y)
species_xmat$x = log(x)
xlim = range(species_xmat$x)*c(0.7, 1.2)
cor1test = with(species_xmat,cor.test(x,y))
cor1 = round(cor1test$estimate, 2)
pv1 = signif(cor1test$p.value,1)
coef1 = coef(lm(y~x, data = species_xmat))
mytitle = paste0(letters[kk],". ",sname,', N=',
length(y),", Cor=",cor1,", p=", pv1)
p1 = ggplot(species_xmat,
aes_string(y = "y", x = "x",
label = "MammalNumberHorvath")) +
geom_point(color="darkgrey") + theme_bw() +
scale_x_continuous(name='Log Lifespan', limits = xlim)+
scale_y_continuous(name=paste0('Log Adjusted',sname))+
ggtitle(mytitle) +
theme(plot.title=element_text(size = titlesize)) + # margin=margin(t=40,b=-30)
geom_text(aes(label=ifelse(textidx,
MammalNumberHorvath,'')),
hjust=0,vjust=0) +
geom_abline(slope=coef1[2], intercept = coef1[1], linetype=3)
pAROCM[[kk]] = p1
kk=kk+1
###
col1 = switch (kk/4,
"YoungSlope1",
"YoungSlope0.1",
"OldSlope1")
# species_xmat = species_xmat%>%filter(!is.na(species_xmat[,k+1]))
# colnames( species_xmat)
y = species_xmat[,k+1]
x = species_xmat[, 15]
lm1 = lm(y~ x-1)
coef1 = as.numeric(coef(lm1))
a = round(coef1[1],2)
species_xmat[,15] = species_xmat[,15]*coef1
# cor1 = round(cor(x,y, use = 'pairwise.complete'),3)
cor1test = cor.test(x,y, use = 'pairwise.complete')
cor1 = round(cor1test$estimate, 2)
pv1 = signif(cor1test$p.value,1)
mytitle = paste0(letters[kk],". ",sname,', N=',length(y)) #,", Cor=",cor1)
# species_xmat$x = x
# species_xmat$y = y
xr = range(c(x, y)) *c(1/2, 5)
#dlist[[kk]] = species_xmat
p1 = ggplot(species_xmat,
aes_string(y = col1, x = col2,
label = "MammalNumberHorvath")) +
geom_point(color="darkgrey") + theme_bw() +
scale_x_log10(name= paste0("a",substring(col2,2),", a=",a),
limits = xr) +
scale_y_log10(name=paste0('Adjusted ',sname),
limits = xr)+
ggtitle(mytitle, subtitle =
paste0("Cor=",cor1, ", p=", pv1 )) +
theme(plot.title=element_text(size = titlesize),
plot.subtitle = element_text(size = 1.5*titlesize,
margin=margin(b=-20))) + # margin=margin(t=40,b=-30)
geom_text(aes(label=ifelse(textidx,
MammalNumberHorvath,'')),
hjust=0,vjust=0) +
geom_abline(slope=1, intercept = 0, linetype=3)
pAROCM[[kk]] = p1
kk=kk+1
}
}
length(pAROCM)
pAROCM[[1]]
```
```{r}
vnum = 'v11'
ggsave(paste0("Figure4_",colnames(SpeciesMat)[nj], "_",vnum, ".pdf"),
grid.arrange(pAROCM[[1]],pAROCM[[3]],
pAROCM[[2]],pAROCM[[4]],
nrow=2, as.table = FALSE),
width = 16,
height = 16, units = "cm")
```
## Figure 7: Barplot Mammalian
```{r}
outfolder = "Fig4_barplot_slopes by species"
slopefolder = "out0314_TissueSlopes"
########### correlations only, Figure 5
fisher_t = function(r){
0.5*log((1+r)/(1-r))
}
###### (0, lifespan)
j = 2
len1 = stateSlopes$statesPRC2[j]
# tmp1 = grep(len1, colnames(spectissuecormat), fixed = TRUE)
# colnames(spectissuecormat)[tmp1]
# slopes_mat = SpeciesMat%>%select(contains(len1),maxAgeCaesar,SpeciesLatinName)
matslopes = read.csv(paste0(slopefolder,"/states/Slopes_IdentityScaled",len1, "v10.csv"),
row.names = 1)
cormat = matslopes%>%select(contains("cor"))
sdmat = matslopes%>%select(contains("sd_rage"))
slopes_mat = matslopes%>%select(contains("Slope"),maxAgeCaesar,SpeciesLatinName)
source("slope_strata_v2.R")
k = 10
sdrage = sd_ragemat[,k]
prop = ifelse(k<= 10, props[k], oldprops[k-10])
LU = ifelse(k<= 10, paste0("(0, ",props[k],"*Lifespan)"),
paste0("(MaturityAge, Lifespan)"))
LU
idx2 = which(SpeciesMat$Freq>= 15)
if(k<= 10){
idx1 = which(!is.na(slopes_mat[,k]) & slopes_mat[,k]<10 &
slopes_mat[,k] > -0.9 &
SpeciesMat$RemoveFei==0 &
SpeciesMat$ageRangeL< matratio &
#SpeciesMat$ageRangeU >= min(prop,0.5) &
!is.na(sdrage) & sdrage> 0.03 )
}else {
idx1 = which(!is.na(slopes_mat[,k]) & slopes_mat[,k]<10 &
(slopes_mat[,k] > 0 | SpeciesMat$maxAgeCaesar >20) &
SpeciesMat$RemoveFei==0 &
!is.na(sdrage) & sdrage>= 0.06 )
}
idx1 = intersect(idx1, idx2)
x = cormat[idx1, k]
y = sdmat[idx1, k]
plot(y,x)
qcod = function(x){
qt_x = quantile(x, probs = c(0.25, 0.75))
(qt_x[2] - qt_x[1])/ (qt_x[2] + qt_x[1])
}
sd(x)/mean(x)
sd(x/y^0.3)/mean(x/y^0.3)
sd(x/y)/mean(x/y)
qcod(x)
qcod(x/y)
qcod(x/y^0.3)
qcods = NULL
for (pw in (0:10)/10) {
qcods = rbind(qcods,
c(qcod(x/y^pw), sd(x/y^pw)/mean(x/y^pw)))
}
tab1 = cbind((0:10)/10, qcods)
rownames(tab1) = NULL
colnames(tab1) = c("pw","Quartile_COD", "COV")
write.csv(round(tab1,3),paste0(outfolder,"/Quartile_COD.csv"))
```
```{r}
data = data.frame(matslopes%>%select(SpeciesLatinName,Tissue,Freq,
MammalNumberHorvath,maxAgeCaesar),
cors = cormat[, k],
sds = sdmat[, k])
# colnames(data)
ncs = 5:7
dat2 = combineTissues(data, idx1 = idx1, nc = ncs)
dat2freq = ddply(data[idx1,], "SpeciesLatinName", function(dat1){
sum(dat1$Freq)
})
colnames(dat2freq)[2] = "Freq"
dat2 = merge(dat2, dat2freq)
# colnames(dat2)
mean(data$cors, na.rm=TRUE)
median(data$cors, na.rm=TRUE)
mean(dat2$cors, na.rm=TRUE)
median(dat2$cors, na.rm=TRUE)
write.csv(data%>%arrange(maxAgeCaesar),paste0(outfolder,"/CorBySpeciesTissue_",LU,"_v2.csv"))
write.csv(dat2%>%arrange(maxAgeCaesar),paste0(outfolder,"/CorBySpecies_",LU,"_v2.csv"))
```
```{r}
############ BAR plot new Figure 4
colnames(matslopes)
cormat = matslopes%>%select(contains("cor_"))
cormat = cormat[,-(1:2)]
colnames(cormat)
sdmat = matslopes%>%select(contains("sd_rage"))
pw = 0.25
ratios = cormat/sdmat^pw
# ratios = fisher_t(ratios)
colnames(ratios) = paste0("pw_",pw,substring(colnames(ratios),4))
k = 10
# k = 11
data = data.frame(matslopes%>%
select(SpeciesLatinName,
Tissue,MammalNumberHorvath,
maxAgeCaesar),
sds = sdmat[, k],
cors = cormat[, k],
ratio = ratios[,k])
{
sdrage = sd_ragemat[,k]
prop = ifelse(k<= 10, props[k], oldprops[k-10])
LU = ifelse(k<= 10, paste0("(0, ",props[k],"*Lifespan)"),
paste0("(MaturityAge, Lifespan)"))
idx2 = which(ratios[,k]>0 & SpeciesMat$Freq>= 15)
if(k<= 10){
idx1 = which(!is.na(slopes_mat[,k]) & slopes_mat[,k]<10 &
slopes_mat[,k] > -0.9 &
SpeciesMat$RemoveFei==0 &
SpeciesMat$ageRangeL< matratio &
#SpeciesMat$ageRangeU >= min(prop,0.5) &
!is.na(sdrage) & sdrage> 0.03 )
}else {
idx1 = which(!is.na(slopes_mat[,k]) & slopes_mat[,k]<10 &
(slopes_mat[,k] > 0 | SpeciesMat$maxAgeCaesar >20) &
SpeciesMat$RemoveFei==0 &
!is.na(sdrage) & sdrage>= 0.06 )
}
idx1 = intersect(idx1, idx2)
# colnames(data)
ncs = 4:7
dat2 = combineTissues(data, idx1 = idx1, nc = ncs)
dat2$SpeciesLatinName[dat2$SpeciesLatinName=="1.3"] = "Lemuridae"
# colnames(dat2)
dat2 = dat2[order(dat2$maxAgeCaesar),]
dat2$rank = 1:nrow(dat2)
# dat2 = dat2%>% mutate(rank = fct_reorder(rank, rank))
# colnames(dat2)
# invmat = apply(dat2[,3:5], 2, invfn)
# pmatch(data$SpeciesLatinName, dat2$SpeciesLatinName, duplicates.ok = TRUE)
m1 = match(data$SpeciesLatinName, dat2$SpeciesLatinName)
# dat2$rank[m1]
data$rank = NA
data$rank[!is.na(m1)] = dat2$rank[m1[!is.na(m1)]]
data$Freq = SpeciesMat$Freq
}
```
```{r}
plist = hlist = list()
#for(k in c(10,11,1)){
for(kk in 1:2){
yname = ifelse(kk==1,
paste0("Cor(Age, Methyl)"),
paste0("Adjusted Cor = Cor/SD^",pw))
data$lower = data[,kk+5] - 2/sqrt(data$Freq - 3)
data$upper = data[,kk+5] + 2/sqrt(data$Freq - 3)
##### Ratio
r1 = colnames(data)[kk+5]
h1 = ggplot(data[idx1,], aes(!!ensym(r1)))+
geom_histogram()
hlist[[kk]] = h1
{
tit1 = ifelse(kk==1, paste0("a. N = ",length(idx1)),
paste0("b. Age Interval ",LU))
tmpm = match(dat2$MammalNumberHorvath, dat1$MammalNumberHorvath)
if(kk == 1) speclabsfull = paste(dat1$MammalNumberHorvath[tmpm],
dat1$SpeciesCommonName[tmpm],sep=" ") else
speclabsfull = NULL
speclabs = dat1$MammalNumberHorvath[tmpm]
# med1 = median(dat2%>%select(all_of(r1)))
p1 = ggplot(dat2) +
geom_bar(aes(y = !!ensym(r1), x = rank), stat="identity", position="dodge",
width = 1,
fill="grey", col="black",na.rm=TRUE) +
coord_flip() + theme_bw() +
geom_hline(data = data,
aes(yintercept = mean(!!ensym(r1), na.rm=TRUE), col="red")) +
labs(y = yname, title = tit1) +
scale_x_continuous(ifelse(kk==1,"Species",""),
breaks = dat2$rank, labels = speclabsfull) +
scale_y_continuous(expand = c(0,0)) +
# scale_fill_continuous(name = "Species", labels = speclabsfull) +
theme(axis.text=element_text(size=8),
axis.ticks.x = element_blank(),
plot.title = element_text(size = 11, hjust = 0.5)) +
geom_point(data = data[idx1,],
aes(rank, !!ensym(r1), col=SpeciesMat$col.tissue[idx1]), na.rm=TRUE)+
theme(legend.position = "none")
plist[[kk]] = p1
# print(p1)
}
}
# grid.arrange(grobs = hlist,ncol=2)
```
```{r}
with(dat2, cor.test(log(cors), log(maxAgeCaesar)))
with(dat2, cor.test(log(ratio), log(maxAgeCaesar)))
with(dat2, plot(log(cors), log(maxAgeCaesar)))
with(dat2, plot(log(ratio), log(maxAgeCaesar)))
cor1 = with(dat2, cor.test(log(cors), log(maxAgeCaesar)))
tit1 = paste0("c. Cor=",round(cor1$estimate,2),", p=",signif(cor1$p.value,2))
plist[[3]] = ggplot(dat2, aes(x = log(maxAgeCaesar),
y = log(cors),
label = MammalNumberHorvath))+
geom_text()+xlab("Log(Lifespan)") + ylab("Log(Cor(Age,Methyl))") +
xlim(0,6) + ggtitle(tit1) + theme_bw()
cor1 = with(dat2, cor.test(log(ratio), log(maxAgeCaesar)))
tit1 = paste0("d. Cor=",round(cor1$estimate,2),", p=",signif(cor1$p.value,2))
plist[[4]] = ggplot(dat2, aes(x = log(maxAgeCaesar),
y = log(ratio),
label = MammalNumberHorvath))+
geom_text()+xlab("Log(Lifespan)") + ylab("Log(Adj.Cor)") +
xlim(0,6) + ggtitle(tit1) + theme_bw()
```
```{r}
# pdf(paste0(outfolder,"/suppFigure1_",len1,"_k=",k,"_v2.pdf"),
pdf(paste0(outfolder,"/Figure4_",len1,"_",invname,"_v6.pdf"),
width = 7,height = 9
)
hlay = matrix(c(1,1,1,1,2,2,
1,1,1,1,2,2,
1,1,1,1,2,2,
1,1,1,1,2,2,
1,1,1,1,2,2,
1,1,1,1,2,2,
1,1,1,1,2,2,
1,1,1,1,2,2,
1,1,1,1,2,2,
1,1,1,1,2,2,
3,3,3,4,4,4,
3,3,3,4,4,4,
3,3,3,4,4,4), nrow = 13,byrow = TRUE)
grid.arrange(plist[[1]], plist[[2]],
plist[[3]], plist[[4]],
nrow=2,layout_matrix=hlay)
dev.off()
```
## Figure S5 - S7
```{r}
titlesize = 10
length(pAROCM)
kk = 5
kkseq = c(1:4, 3, 1, 4, 2,
7, 5, 8, 6)
# for(k in c(10,1,11)){
for(k in c(1,11)){
sname = ifelse(k==10, "AROCM",
ifelse(k==11, "Old AROCM", "Young AROCM"))
prop = ifelse(k<= 10, props[k], oldprops[k-10])
tit1 = ifelse(k<= 10, paste0(", (0, ",props[k],"*Lifespan)"),
paste0(", (Maturity, Lifespan)"))
if(k<= 10){
idx1 = which(!is.na(slopes_mat[,k]) & slopes_mat[,k]<10 &
slopes_mat[,k] > -0.9 &
SpeciesMat$RemoveFei==0
# SpeciesMat$ageRangeL< matratio &
# & SpeciesMat$ageRangeU >= min(prop,0.6)
)
}else {
sdrage = sd_ragemat[,k]
idx1 = which(!is.na(slopes_mat[,k]) & slopes_mat[,k]<10 &
slopes_mat[,k] > -0.1 &
SpeciesMat$RemoveFei==0 &
!is.na(sdrage) & sdrage>= 0.06 )
}
ncs = 1:14 ###c(k, grep("maxAgeCaesar",colnames(slopes_mat)))
speciesSlopes = combineTissues(slopes_mat, idx1, ncs)
x = speciesSlopes[,15]
y = speciesSlopes[,1+k]
y = logplus(y)
speciesSlopes$y = log(y)
speciesSlopes$x = log(x)
xlim = range(speciesSlopes$x)*c(0.7, 1.2)
cor1test = with(speciesSlopes,cor.test(x,y))
cor1 = round(cor1test$estimate, 2)
pv1 = signif(cor1test$p.value,1)
coef1 = coef(lm(y~x, data = speciesSlopes))
speciesSlopes$textidx =
speciesSlopes$MammalNumberHorvath %in% textnum
mytitle = paste0(letters[kkseq[kk]],". ",sname,
", Cor=",cor1,", p=", pv1)
p1 = ggplot(speciesSlopes,
aes_string(y = "y", x = "x",
label = "MammalNumberHorvath")) +
geom_point(color="darkgrey") + theme_bw() +
scale_x_continuous(name='Log Lifespan', limits = xlim)+
scale_y_continuous(name=paste0('Log ',sname))+
ggtitle(mytitle) +
theme(plot.title=element_text(size = titlesize)) + # margin=margin(t=40,b=-30)
geom_text(aes(label=ifelse(textidx,
MammalNumberHorvath,'')),
hjust=0,vjust=0) +
geom_abline(slope=coef1[2], intercept = coef1[1], linetype=3)
pAROCM[[kk]] = p1
kk=kk+1
##### panel b
{
b=1
tmpx = (x)^(-b)
lm1 = lm(y~ tmpx-1, weights = tmpx)
coef1 = as.numeric(coef(lm1))
x1 = coef1[1]*(x)^(-b)
a = round(coef1[1],2)
speciesSlopes$x = a/x
speciesSlopes$y = y
cor1test = with(speciesSlopes,cor.test(x,y))
cor1 = round(cor1test$estimate, 2)
pv1 = signif(cor1test$p.value,1)
coef1 = coef(lm(y~x, data = speciesSlopes))
mytitle = paste0(letters[kkseq[kk]],". ",sname,', N=',
length(y))
xr = range(c(speciesSlopes$x, speciesSlopes$y)) *c(1/2, 5)
p1 = ggplot(speciesSlopes,
aes_string(y = "y", x = "x",
label = "MammalNumberHorvath")) +
geom_point(color="darkgrey") + theme_bw() +
scale_x_log10(name=paste0("a/",colnames(speciesSlopes)[15],
", a = ",a), limits = xr) + ## , limits = xlim
scale_y_log10(name=paste0(sname), limits = xr) +
labs(title = mytitle, subtitle =
paste0("Cor=",cor1, ", p=",pv1 )) +
theme(plot.title=element_text(size = titlesize),
plot.subtitle = element_text(size = 1.5*titlesize,
margin=margin(b=-20))) + # margin = margin(t=-5)
geom_text(aes(label=ifelse(textidx,
MammalNumberHorvath,'')),
hjust=0,vjust=0) +
geom_abline(slope=coef1[2], intercept = coef1[1], linetype=3)
pAROCM[[kk]] = p1
kk=kk+1
}
#### panel c & d
{
xmat = cbind(slopes_mat[,1:13]*sd_ragemat^pw,
1/SpeciesMat[,nj], slopes_mat[,14:15])
colnames(xmat)[14] = col2 = paste0("1/", colnames(SpeciesMat)[nj])
ncs = 1:15 ###c(k, grep("maxAgeCaesar",colnames(slopes_mat)))
species_xmat = combineTissues(xmat, idx1, ncs)
species_xmat$textidx = species_xmat$MammalNumberHorvath %in% textnum
## log scale panel
x = species_xmat[, 16]
y = species_xmat[,k+1]
y = logplus(y)
species_xmat$y = log(y)
species_xmat$x = log(x)
xlim = range(species_xmat$x)*c(0.7, 1.2)
cor1test = with(species_xmat,cor.test(x,y))
cor1 = round(cor1test$estimate, 2)
pv1 = signif(cor1test$p.value,1)
coef1 = coef(lm(y~x, data = species_xmat))
mytitle = paste0(letters[kkseq[kk]],". ",sname,
", Cor=",cor1,", p=", pv1)
p1 = ggplot(species_xmat,
aes_string(y = "y", x = "x",
label = "MammalNumberHorvath")) +
geom_point(color="darkgrey") + theme_bw() +
scale_x_continuous(name='Log Lifespan', limits = xlim)+
scale_y_continuous(name=paste0('Log Adjusted',sname))+
ggtitle(mytitle) +
theme(plot.title=element_text(size = titlesize)) + # margin=margin(t=40,b=-30)
geom_text(aes(label=ifelse(textidx,
MammalNumberHorvath,'')),
hjust=0,vjust=0) +
geom_abline(slope=coef1[2], intercept = coef1[1], linetype=3)
pAROCM[[kk]] = p1
kk=kk+1
###
col1 = switch (kk/4,
"YoungSlope1",
"YoungSlope0.1",
"OldSlope1")
# species_xmat = species_xmat%>%filter(!is.na(species_xmat[,k+1]))
# colnames( species_xmat)
y = species_xmat[,k+1]
x = species_xmat[, 15]
lm1 = lm(y~ x-1)
coef1 = as.numeric(coef(lm1))
a = round(coef1[1],2)
species_xmat[,15] = species_xmat[,15]*coef1
# cor1 = round(cor(x,y, use = 'pairwise.complete'),3)
cor1test = cor.test(x,y, use = 'pairwise.complete')
cor1 = round(cor1test$estimate, 2)
pv1 = signif(cor1test$p.value,1)
mytitle = paste0(letters[kkseq[kk]],". ",
sname,', N=',length(y))
# species_xmat$x = x
# species_xmat$y = y
xr = range(c(x, y)) *c(1/2, 5)
#dlist[[kk]] = species_xmat
p1 = ggplot(species_xmat,
aes_string(y = col1, x = col2,
label = "MammalNumberHorvath")) +
geom_point(color="darkgrey") + theme_bw() +
scale_x_log10(name= paste0("a",substring(col2,2),", a=",a),
limits = xr) +
scale_y_log10(name=paste0('Adjusted ',sname),
limits = xr)+
ggtitle(mytitle, subtitle =
paste0("Cor=",cor1, ", p=", pv1 )) +
theme(plot.title=element_text(size = titlesize),
plot.subtitle = element_text(size = 1.5*titlesize,
margin=margin(b=-20))) + # margin=margin(t=40,b=-30)
geom_text(aes(label=ifelse(textidx,
MammalNumberHorvath,'')),
hjust=0,vjust=0) +
geom_abline(slope=1, intercept = 0, linetype=3)
pAROCM[[kk]] = p1
kk=kk+1
}
}
length(pAROCM)
pAROCM[[12]]
```
```{r}
vnum = 'v10'
ggsave(paste0("FigS6_",colnames(SpeciesMat)[nj], "_",vnum, ".pdf"),
grid.arrange(pAROCM[[6]],pAROCM[[5]],
pAROCM[[10]],pAROCM[[9]],
pAROCM[[8]],pAROCM[[7]],
pAROCM[[12]],pAROCM[[11]],
nrow=4, as.table = FALSE),
width = 16,
height = 24, units = "cm")
```
### ASM plot (nj=11)
```{r}
nj = 11
colnames(SpeciesMat)[nj]
slopes_mat[, 14] = SpeciesMat[,nj]
colnames(slopes_mat)[14] = colnames(SpeciesMat)[nj]
```
```{r}
###### new text
textnum = c("1.1.1","4.11.1","9.9.1","9.9.3","12.3.10",
"5.1.1", "1.4.1","6.1.1", "2.1.1", "8.4.5" )
pw = 0.75
kk=1
plist = list()
titlesize = 10
for(k in c(10,1,11)){
sname = ifelse(k==10, "AROCM",
ifelse(k==11, "Old AROCM", "Young AROCM"))
prop = ifelse(k<= 10, props[k], oldprops[k-10])
tit1 = ifelse(k<= 10, paste0(", (0, ",props[k],"*Lifespan)"),
paste0(", (Maturity, Lifespan)"))
if(k<= 10){
idx1 = which(!is.na(slopes_mat[,k]) & slopes_mat[,k]<10 &
slopes_mat[,k] > -0.9 &
SpeciesMat$RemoveFei==0 &
SpeciesMat$ageRangeL< matratio &
SpeciesMat$ageRangeU >= min(prop,0.6) )
# if(k> 5) idx1 = setdiff(c(idx1, 191), 192) else idx1 = setdiff(c(idx1, 192), 191)
# if(k>= 5) idx1 = setdiff(idx1, 191)
}else {
sdrage = sd_ragemat[,k]
idx1 = which(!is.na(slopes_mat[,k]) & slopes_mat[,k]<10 &
slopes_mat[,k] > -0.1 &
SpeciesMat$RemoveFei==0 &
!is.na(sdrage) & sdrage>= 0.06 )
}
ncs = 1:14 ###c(k, grep("maxAgeCaesar",colnames(slopes_mat)))
speciesSlopes = combineTissues(slopes_mat, idx1, ncs)
x = speciesSlopes[,15]
y = speciesSlopes[,1+k]
y = logplus(y)
speciesSlopes$y = log(y)
speciesSlopes$x = log(x)
xlim = range(speciesSlopes$x)*1.2
cor1 = round(with(speciesSlopes,cor(x,y)), 3)
coef1 = coef(lm(y~x, data = speciesSlopes))
speciesSlopes$textidx = speciesSlopes$MammalNumberHorvath %in% textnum
mytitle = paste0(letters[kk],". ",sname,', N=',
length(y),", Cor=",cor1) ##letters[kk],". ",
p1 = ggplot(speciesSlopes,
aes_string(y = "y", x = "x",
label = "MammalNumberHorvath")) +
geom_point(color="darkgrey") + theme_bw() +
scale_x_continuous(name='Log Lifespan', limits = xlim)+
scale_y_continuous(name=paste0('Log ',sname))+
ggtitle(mytitle) +
theme(plot.title=element_text(size = titlesize)) + # margin=margin(t=40,b=-30)
geom_text(aes(label=ifelse(textidx,
MammalNumberHorvath,'')),
hjust=0,vjust=0) +
geom_abline(slope=coef1[2], intercept = coef1[1], linetype=3)
plist[[kk]] = p1
kk=kk+1
##### panel b
{
b=1
tmpx = (x)^(-b)
lm1 = lm(y~ tmpx-1, weights = tmpx)
coef1 = as.numeric(coef(lm1))
x1 = coef1[1]*(x)^(-b)
a = round(coef1[1],2)
speciesSlopes$x = a/speciesSlopes[,15]
speciesSlopes$y = y
cor1 = round(with(speciesSlopes,cor(x,y)), 3)
coef1 = coef(lm(y~x, data = speciesSlopes))
mytitle = paste0(letters[kk],". ",sname,', N=',
length(y)) ##letters[kk],". ",
xr = range(c(speciesSlopes$x, speciesSlopes$y)) *c(1/2, 5)
p1 = ggplot(speciesSlopes,
aes_string(y = "y", x = "x",
label = "MammalNumberHorvath")) +
geom_point(color="darkgrey") + theme_bw() +
scale_x_log10(name=paste0("a/",colnames(speciesSlopes)[15],
", a = ",a), limits = xr) + ## , limits = xlim
scale_y_log10(name=paste0(sname), limits = xr) +
labs(title = mytitle, subtitle = paste0("Cor=",cor1 )) +
theme(plot.title=element_text(size = titlesize),
plot.subtitle = element_text(size = 1.5*titlesize,
margin=margin(b=-20))) + # margin = margin(t=-5)
geom_text(aes(label=ifelse(textidx,
MammalNumberHorvath,'')),
hjust=0,vjust=0) +
geom_abline(slope=coef1[2], intercept = coef1[1], linetype=3)
plist[[kk]] = p1
kk=kk+1
}
#### panel c
{
xmat = cbind(slopes_mat[,1:13]*sd_ragemat^pw,
1/SpeciesMat[,nj], slopes_mat[,14:15])
colnames(xmat)[14] = col2 = paste0("1/", colnames(SpeciesMat)[nj])
ncs = 1:15 ###c(k, grep("maxAgeCaesar",colnames(slopes_mat)))
species_xmat = combineTissues(xmat, idx1, ncs)
species_xmat$textidx = species_xmat$MammalNumberHorvath %in% textnum
col1 = switch (kk/3,
"YoungSlope1",
"YoungSlope0.1",
"OldSlope1")
# species_xmat = species_xmat%>%filter(!is.na(species_xmat[,k+1]))
# colnames( species_xmat)
y = species_xmat[,k+1]
x = species_xmat[, 15]
lm1 = lm(y~ x-1)
coef1 = as.numeric(coef(lm1))
a = round(coef1[1],2)
species_xmat[,15] = species_xmat[,15]*coef1
cor1 = round(cor(x,y, use = 'pairwise.complete'),3)
mytitle = paste0(letters[kk],". ",sname,', N=',length(y)) #,", Cor=",cor1)
#dlist[[kk]] = species_xmat
p1 = ggplot(species_xmat,
aes_string(y = col1, x = col2,
label = "MammalNumberHorvath")) +
geom_point(color="darkgrey") + theme_bw() +
scale_x_log10(name= paste0("a",substring(col2,2),", a=",a))+ ## , limits = c(0.001,3)
scale_y_log10(name=paste0('Adjusted ',sname))+
ggtitle(mytitle, subtitle = paste0("Cor=",cor1 )) +
theme(plot.title=element_text(size = titlesize),
plot.subtitle = element_text(size = 1.5*titlesize,
margin=margin(b=-20))) + # margin=margin(t=40,b=-30)
geom_text(aes(label=ifelse(textidx,
MammalNumberHorvath,'')),
hjust=0,vjust=0) +
geom_abline(slope=1, intercept = 0, linetype=3)
plist[[kk]] = p1
kk=kk+1
}
}
length(plist)
```
```{r}
vnum = 'v4'
ggsave(paste0(outfolder,"/Figure6_",colnames(SpeciesMat)[nj], "_",vnum, ".pdf"),
grid.arrange(plist[[1]],plist[[2]],plist[[3]],
plist[[4]],plist[[5]],plist[[6]],
plist[[7]],plist[[8]],plist[[9]],
nrow=3, as.table = FALSE),
width = 24,
height = 24, units = "cm")
```
### ASM and GT combined plot
```{r}
###### new text
textnum = c("1.1.1","4.11.1","9.9.1","9.9.3","12.3.10",
"5.1.1", "1.4.1","6.1.1", "2.1.1", "8.4.5" )
pw = 0.8
kk=1
pGTASM = list()
titlesize = 10
for (nj in c(10:11)) {
njname = colnames(SpeciesMat)[nj]
slopes_mat[, 14] = SpeciesMat[,nj]
colnames(slopes_mat)[14] = colnames(SpeciesMat)[nj]
k= 1
sname = ifelse(k==10, "AROCM",
ifelse(k==11, "Old AROCM", "Young AROCM"))
prop = ifelse(k<= 10, props[k], oldprops[k-10])
tit1 = ifelse(k<= 10, paste0(", (0, ",props[k],"*Lifespan)"),
paste0(", (Maturity, Lifespan)"))
if(k<= 10){
idx1 = which(!is.na(slopes_mat[,k]) & slopes_mat[,k]<10 &
slopes_mat[,k] > -0.9 &
SpeciesMat$RemoveFei==0 &
SpeciesMat$ageRangeL< matratio &
SpeciesMat$ageRangeU >= min(prop,0.6) )
# if(k> 5) idx1 = setdiff(c(idx1, 191), 192) else idx1 = setdiff(c(idx1, 192), 191)
# if(k>= 5) idx1 = setdiff(idx1, 191)
}else {
sdrage = sd_ragemat[,k]
idx1 = which(!is.na(slopes_mat[,k]) & slopes_mat[,k]<10 &
slopes_mat[,k] > -0.1 &
SpeciesMat$RemoveFei==0 &
!is.na(sdrage) & sdrage>= 0.06 )
}
ncs = 1:14 ###c(k, grep("maxAgeCaesar",colnames(slopes_mat)))
speciesSlopes = combineTissues(slopes_mat, idx1, ncs)
x = speciesSlopes[,15]
y = speciesSlopes[,1+k]
y = logplus(y)
speciesSlopes$y = log(y)
speciesSlopes$x = log(x)
xlim = range(speciesSlopes$x) + c(-1,1)
cor1test = with(speciesSlopes,cor.test(x,y))
cor1 = round(cor1test$estimate, 2)
pv1 = signif(cor1test$p.value,1)
coef1 = coef(lm(y~x, data = speciesSlopes))
speciesSlopes$textidx =
speciesSlopes$MammalNumberHorvath %in% textnum
mytitle = paste0(letters[kk],". ",sname,', N=',
length(y),", Cor=",cor1,", p=",pv1)
p1 = ggplot(speciesSlopes,
aes_string(y = "y", x = "x",
label = "MammalNumberHorvath")) +
geom_point(color="darkgrey") + theme_bw() +
scale_x_continuous(name=paste0('Log ',njname), limits = xlim)+
scale_y_continuous(name=paste0('Log ',sname))+
ggtitle(mytitle) +
theme(plot.title=element_text(size = titlesize)) + # margin=margin(t=40,b=-30)
geom_text(aes(label=ifelse(textidx,
MammalNumberHorvath,'')),
hjust=0,vjust=0) +
geom_abline(slope=coef1[2], intercept = coef1[1], linetype=3)
pGTASM[[kk]] = p1
kk=kk+1
##### panel b
{
b=1
tmpx = (x)^(-b)
lm1 = lm(y~ tmpx-1, weights = tmpx)
coef1 = as.numeric(coef(lm1))
x1 = coef1[1]*(x)^(-b)
a = round(coef1[1],2)
speciesSlopes$x = a/speciesSlopes[,15]
speciesSlopes$y = y
cor1test = with(speciesSlopes,cor.test(x,y))
cor1 = round(cor1test$estimate, 2)
pv1 = signif(cor1test$p.value,1)
coef1 = coef(lm(y~x, data = speciesSlopes))
mytitle = paste0(letters[kk],". ",sname,', N=',
length(y))
xr = range(c(speciesSlopes$x, speciesSlopes$y)) *c(1/2, 5)
p1 = ggplot(speciesSlopes,
aes_string(y = "y", x = "x",
label = "MammalNumberHorvath")) +
geom_point(color="darkgrey") + theme_bw() +
scale_x_log10(name=paste0("a/",colnames(speciesSlopes)[15],
", a = ",a), limits = xr) + ## , limits = xlim
scale_y_log10(name=paste0(sname), limits = xr) +
labs(title = mytitle, subtitle =
paste0("Cor=",cor1,", p=",pv1 )) +
theme(plot.title=element_text(size = titlesize),
plot.subtitle = element_text(size = 1.5*titlesize,
margin=margin(b=-20))) + # margin = margin(t=-5)
geom_text(aes(label=ifelse(textidx,
MammalNumberHorvath,'')),
hjust=0,vjust=0) +
geom_abline(slope=coef1[2], intercept = coef1[1], linetype=3)
pGTASM[[kk]] = p1
kk=kk+1
}
}
length(pGTASM)
```
```{r}
vnum = 'v4'
ggsave(paste0("SuppFigure6_ASM_GT_",vnum, ".pdf"),
grid.arrange(pGTASM[[1]],pGTASM[[2]],
pGTASM[[3]],
pGTASM[[4]],
nrow=2, as.table = TRUE),
width = 18,
height = 18, units = "cm")
```
### Adult weight
```{r}
slopes_mat[, 14] = SpeciesMat$weightCaesar
colnames(slopes_mat)[14] = "Weight"
```
```{r}
###### new text
textnum = c("1.1.1","4.11.1","9.9.1","9.9.3","12.3.10",
"5.1.1", "1.4.1","6.1.1", "2.1.1", "8.4.5" )
kk=1
pWeight = list()
titlesize = 10
for(k in c(10,1,11)){
sname = ifelse(k==10, "AROCM",
ifelse(k==11, "Old AROCM", "Young AROCM"))
prop = ifelse(k<= 10, props[k], oldprops[k-10])
tit1 = ifelse(k<= 10, paste0(", (0, ",props[k],"*Lifespan)"),
paste0(", (Maturity, Lifespan)"))
if(k<= 10){
idx1 = which(!is.na(slopes_mat[,k]) & slopes_mat[,k]<10 &
slopes_mat[,k] > -0.9 &
SpeciesMat$RemoveFei==0 &
SpeciesMat$ageRangeL< matratio &
SpeciesMat$ageRangeU >= min(prop,0.6) )
# if(k> 5) idx1 = setdiff(c(idx1, 191), 192) else idx1 = setdiff(c(idx1, 192), 191)
# if(k>= 5) idx1 = setdiff(idx1, 191)
}else {
sdrage = sd_ragemat[,k]
idx1 = which(!is.na(slopes_mat[,k]) & slopes_mat[,k]<10 &
slopes_mat[,k] > -0.1 &
SpeciesMat$RemoveFei==0 &
!is.na(sdrage) & sdrage>= 0.06 )
}
# colnames(slopes_mat)
ncs = 1:14
speciesSlopes = combineTissues(slopes_mat, idx1, ncs)
x = speciesSlopes$Weight
y = speciesSlopes[,1+k]
y = logplus(y)
##### panel a
{
speciesSlopes$x = log(speciesSlopes$Weight)
speciesSlopes$y = log(y)
speciesSlopes$textidx = speciesSlopes$MammalNumberHorvath %in% textnum
xr = range(speciesSlopes$x)*c(0, 1.2)
cor1test = with(speciesSlopes,cor.test(x,y))
cor1 = round(cor1test$estimate, 2)
pv1 = signif(cor1test$p.value,1)
coef1 = coef(lm(y~x, data = speciesSlopes))
mytitle = paste0(letters[kk],". ",sname)
##letters[kk],". ", ,', N=', length(y)
p1 = ggplot(speciesSlopes,
aes_string(y = "y", x = "x")) +
geom_point(color="darkgrey") + theme_bw() +
scale_x_continuous(name='Log Weight', limits = xr) + ##
scale_y_continuous(name=paste0('Log ',sname)) +
labs(title = paste0(mytitle, ", Cor=",cor1, ", p=",pv1) ) +
theme(plot.title=element_text(size = titlesize)) +
geom_text(aes(label=ifelse(textidx,
MammalNumberHorvath,'')),
hjust=0,vjust=0)
pWeight[[kk]] = p1
kk=kk+1
}
##### panel b
{
b=1
tmpx = (x)^(-b)
lm1 = lm(y~ tmpx-1, weights = tmpx)
coef1 = as.numeric(coef(lm1))
x1 = coef1[1]*(x)^(-b)
a = round(coef1[1],2)
speciesSlopes$x = 1/speciesSlopes$Weight
speciesSlopes$y = y
speciesSlopes$textidx = speciesSlopes$MammalNumberHorvath %in% textnum
cor1test = with(speciesSlopes,cor.test(x,y))
cor1 = round(cor1test$estimate, 2)
pv1 = signif(cor1test$p.value,1)
coef1 = coef(lm(y~x, data = speciesSlopes))
mytitle = paste0(letters[kk],". ",sname,', N=',
length(y)) ##letters[kk],". ",
# xr = range(c(speciesSlopes$x, speciesSlopes$y))*c(.9, 2)
xr = range(speciesSlopes$x) + c(0, 0.1)
p1 = ggplot(speciesSlopes,
aes_string(y = "y", x = "x")) +
geom_point(color="darkgrey") + theme_bw() +
scale_x_continuous(name='1/Weight', limits = xr) + ##
scale_y_continuous(name=paste0(sname)) +
labs(title = mytitle, subtitle =
paste0("Cor=",cor1,", p=",pv1)) +
theme(plot.title=element_text(size = titlesize),
plot.subtitle = element_text(size = 1.5*titlesize,
margin=margin(b=-20))) +
geom_text(aes(label=ifelse(textidx,
MammalNumberHorvath,'')),
hjust=0,vjust=0)
pWeight[[kk]] = p1
kk=kk+1
}
}
length(pWeight)
```
```{r}
vnum = 'v4'
ggsave(paste0("SuppFig_Weight_",vnum, "_species.pdf"),
grid.arrange(pWeight[[1]],pWeight[[2]],pWeight[[3]],
pWeight[[4]],pWeight[[5]],pWeight[[6]],
nrow=2, as.table = FALSE),
width = 24,
height = 16, units = "cm")
```
## Figure S9-S10
```{r}
library(ggplot2)
xlab = "PIC of Lifespan"
ylab = "PIC of AROCM"
ARtitle = c(outer(c("AROCM","AdjAROCM"),c("","_Young","_Old"), paste0))
ARtitle
ar1 = read.csv("PGM/AROCM_BivProm2. phyloEWAS.results.csv",
row.names = 1)
cor1 = with(ar1, cor.test(Lifespan,slopes))
cor1$estimate
ggplot(ar1, aes(Lifespan,slopes)) +
geom_point() + geom_smooth(method = "lm") +
labs(title = ARtitle[1],
subtitle = paste0("cor=",round(cor1$estimate,2),
", p=",signif(cor1$p.value,2)),
x=xlab, y=ylab) +
theme_bw()
# theme(plot.title = element_text(hjust = 0.5))
```
```{r}
pname = "PGM/"
AROCM_files = list.files(path = pname, pattern = "phyloEWAS")
AROCM_files
fileorder = c(4,1,6,3,5,2)
ARtitle
```
```{r}
p_pgm = list()
for (k in 1:6) {
ar1 = read.csv(paste0(pname,AROCM_files[fileorder[k]]),
row.names = 1)
tit1 = paste0(ARtitle[k], ", N=", nrow(ar1) )
cor1 = with(ar1, cor.test(Lifespan,slopes))
# cor1$estimate
p_pgm[[k]] = ggplot(ar1, aes(Lifespan,slopes)) +
geom_point() + geom_smooth(method = "lm") +
labs(title = paste0(letters[k],". ", tit1),
subtitle = paste0("cor=",round(cor1$estimate,2),
", p=",signif(cor1$p.value,2)),
x=xlab, y=ylab) +
theme_bw()
}
p_pgm[[k]]
```
```{r}
library(gridExtra)
ggsave(
"PGM/AROCM_PIC_v2.pdf",
grid.arrange(grobs = p_pgm, ncol=3, nrow=2, as.table = FALSE),
width = 10, height = 6
)
```
## Figure S11-S12
```{r}
SpeciesMat = read.csv("/Users/feiz/Documents/Box Sync/Horvath/Mammalian meth projects/Slope/out_AROCM_0.1L/states/Slopes_mapped_v1.csv",
row.names = 1, check.names = FALSE)
len1 = "Mapped"
# source("slope_strata_v2.R")
# dim(slopes_mat)
slopes_mat = SpeciesMat%>%select(contains(len1),
maxAgeCaesar,SpeciesLatinName,Tissue)
colnames(slopes_mat)
slopes_mat = slopes_mat[,-(1:5)]
```
### ragemat
```{r}
ragemat = matslopes%>%select(contains("rage"))
sd_ragemat = ragemat%>%select(contains("sd_rage"))
colnames(sd_ragemat)
# k=1
# l=11
# sdrage = sd_ragemat[,l]
```
### Lifespan (nj=12)
```{r}
## outfolder = "Fig6_AROCM"
colnames(SpeciesMat)[1:16]
colnames(SpeciesMat)[10:12] = c("GT", "ASM","Lifespan")
dim(SpeciesMat)
dim(slopes_mat)
colnames(slopes_mat)
identical(SpeciesMat$SpeciesLatinName, slopes_mat$SpeciesLatinName)
```
```{r}
nj = 12 ## Lifespan ## 10:GT 11:ASM
colnames(SpeciesMat)[nj] = "Lifespan"
colnames(slopes_mat)[1:13] = substring(colnames(slopes_mat)[1:13],7)
```
### Choose power
```{r}
len1
matslopes = read.csv(paste0("/Users/feiz/Documents/Box Sync/Horvath/Mammalian meth projects/Slope/out_AROCM_0.1L/states/Slopes_mapped_maplist.csv"),
row.names = 1)
# colnames(matslopes)
cormat = matslopes%>%select(contains("cor_"))
colnames(cormat)
cormat = cormat[,-(1:2)]
sd_ragemat = matslopes%>%select(contains("sd_rage"))
colnames(sd_ragemat)
```
```{r}
k = 10 ## AROCM
datall = data.frame(
matslopes%>%
select(SpeciesLatinName,Tissue,Freq,
MammalNumberHorvath,maxAgeCaesar),
AROCM = matslopes$YoungSlope1,
cors = cormat[, k],
sd_rage = sd_ragemat[, k])
colnames(datall)
hist(datall$sd_rage)
```
```{r}
ggplot(datall,aes(1/maxAgeCaesar, AROCM))+
geom_point()
datfilter = datplot = datall%>%filter(sd_rage> 0.04)
ggplot(datplot,aes(1/maxAgeCaesar, AROCM))+
geom_point()
```
- species level
```{r}
library(plyr)
datplot = datall%>%filter(!is.na(AROCM))
colnames(datplot)
ncs = 5:8
dat2 = combineTissues(datplot, idx1 = 1:nrow(datplot), nc = ncs)
dat2freq = plyr::ddply(datplot, "SpeciesLatinName", function(dat1){
sum(dat1$Freq)
})
colnames(dat2freq)[2] = "Freq"
dat2 = merge(dat2, dat2freq)
ggplot(dat2,aes(1/maxAgeCaesar, AROCM))+
geom_point()
```
```{r}
# datplot = dat2
# datplot = datfilter
range(datplot$sd_rage)
pvec = c(0:40)*0.05
CV = QCOD = Cor_SD_Rage = NULL
for (p in pvec) {
adjROC = datplot$AROCM * (datplot$sd_rage)^(1 - p)
adjCor = datplot$cors / (datplot$sd_rage)^p
CV = c(CV, sd(adjCor)/mean(adjCor))
QCOD = c(QCOD, qcod(adjCor) )
Cor_SD_Rage = c(Cor_SD_Rage, cor(adjCor, datplot$sd_rage) )
}
adj_pmat = data.frame(
pvec,
CV,
QCOD,
Cor_SD_Rage
)
# View(adj_pmat)
```
```{r}
write.csv(adj_pmat,
paste0("../out_AROCM_0.1L/adjROC_select_power_",len1, ".csv"))
# adj_pmat_st = adj_pmat
```
#### supp plot other chromatin states
```{r}
colnames(adj_pmat)
mp = max(abs(adj_pmat$CV))/max(abs(adj_pmat$Cor_SD_Rage))
pw_mapped = ggplot(adj_pmat) + ## %>%filter(pvec> -0.5)
geom_line(aes(pvec,CV, color="CV")) +
geom_line(aes(pvec,QCOD, color="QCOD") ) +
geom_line(aes(pvec,Cor_SD_Rage*mp, color="Cor(Adj.Cor,SD)")) +
geom_hline(yintercept = 0,linetype=2) +
geom_vline(xintercept = 0,linetype=2) +
geom_vline(xintercept = 0.85,linetype=2,col="red") +
scale_color_manual(name='Measure',
breaks=c('CV', 'QCOD', 'Cor(Adj.Cor,SD)'),
values=c('black', 'blue', 'purple'))+
xlab("Power") + ggtitle("b. Mammalian data") +
scale_y_continuous(name = "CV/QCOD",
sec.axis = sec_axis(~ ./mp, name = "Cor")) +
theme_bw() + theme(legend.position = "")
pw_mapped
```
```{r}
pdf(paste0("../out_AROCM_0.1L/adjROC_select_power_",len1,".pdf"),
width = 4.5, height = 4.5)
pw_mapped
dev.off()
```
### Excluding outliers
```{r}
###### new text
textnum = c("1.1.1","4.11.1","9.9.1","12.3.10", ## "9.9.3",
"5.1.1", "1.4.1","6.1.1", "2.1.1", "8.4.5" )
pw = 0.5
len1 = "Mapped"
SpeciesMat$RemoveFei[SpeciesMat$MammalNumberHorvath %in%
c("6.1.2")] = 1
SpeciesMat$RemoveFei[SpeciesMat$MammalNumberHorvath %in%
c("4.13.3") & SpeciesMat$Tissue=="Skin"] = 1
```
```{r}
kk=1
p_mapped = list()
titlesize = 10
# for(k in c(10,1,11)){
for(k in c(10)){
sname = ifelse(k==10, "AROCM",
ifelse(k==11, "Old AROCM", "Young AROCM"))
# sname = len1
prop = ifelse(k<= 10, props[k], oldprops[k-10])
tit1 = ifelse(k<= 10, paste0(", (0, ",props[k],"*Lifespan)"),
paste0(", (Maturity, Lifespan)"))
if(k<= 10){
idx1 = which(!is.na(slopes_mat[,k]) & slopes_mat[,k]<10 &
# slopes_mat[,k] > -0.9 &
SpeciesMat$RemoveFei==0
# SpeciesMat$ageRangeL< matratio &
# & SpeciesMat$ageRangeU >= min(prop,0.6)
)
}else {
sdrage = sd_ragemat[,k]
idx1 = which(!is.na(slopes_mat[,k]) & slopes_mat[,k]<10 &
slopes_mat[,k] > -0.1 &
SpeciesMat$RemoveFei==0 &
!is.na(sdrage) & sdrage>= 0.06 )
}
ncs = 1:14 ###c(k, grep("maxAgeCaesar",colnames(slopes_mat)))
speciesSlopes = combineTissues(slopes_mat, idx1, ncs)
x = speciesSlopes[,15]
y = speciesSlopes[,1+k]
# y = logplus(y)
# y = y - min(y) + 1
speciesSlopes$y = y
speciesSlopes$x = x
xlim = range(speciesSlopes$x)*c(0.7, 1.2)
cor1test = with(speciesSlopes,cor.test(x,y))
cor1 = signif(cor1test$estimate, 2)
pv1 = signif(cor1test$p.value,2)
coef1 = coef(lm(y~x, data = speciesSlopes))
speciesSlopes$textidx =
speciesSlopes$MammalNumberHorvath %in% textnum
mytitle = paste0(letters[kk],". ", len1,
", Cor=",cor1, ", p=",pv1)
p1 = ggplot(speciesSlopes,
aes_string(y = "y", x = "x",
label = "MammalNumberHorvath")) +
geom_point(color="darkgrey") + theme_bw() +
scale_x_continuous(name='Lifespan', limits = xlim)+
scale_y_continuous(name=paste0("AROCM"))+
ggtitle(mytitle) +
theme(plot.title=element_text(size = titlesize)) + # margin=margin(t=40,b=-30)
geom_text(aes(label=ifelse(textidx,
MammalNumberHorvath,'')),
hjust=0,vjust=0) +
geom_hline(yintercept = 0, linetype=3)
# geom_abline(slope=coef1[2], intercept = coef1[1], linetype=3)
p_mapped[[kk]] = p1
kk=kk+1
##### panel b
{
speciesSlopes$x = 1/x
speciesSlopes$y = y
cor1test = with(speciesSlopes,cor.test(x,y))
cor1 = signif(cor1test$estimate, 2)
pv1 = signif(cor1test$p.value,2)
# coef1 = coef(lm(y~x, data = speciesSlopes))
mytitle = paste0(letters[kk],". ", len1,', N=',
length(y)) ##letters[kk],". ",
# xr = range(c(speciesSlopes$x, speciesSlopes$y)) *c(1/2, 5)
xr = range(speciesSlopes$x) + c(0,0.2)
p1 = ggplot(speciesSlopes,
aes_string(y = "y", x = "x",
label = "MammalNumberHorvath")) +
geom_point(color="darkgrey") + theme_bw() +
labs(title = mytitle, subtitle =
paste0("Cor=",cor1,", p=",pv1 )) +
xlab("1/Lifespan") + ylab("AROCM") + xlim(xr) +
theme(plot.title=element_text(size = titlesize),
plot.subtitle = element_text(size = 1.2*titlesize,
margin=margin(b=-15))) +
geom_text(aes(label=ifelse(textidx,
MammalNumberHorvath,'')),
hjust=0,vjust=0) +
geom_hline(yintercept = 0, linetype=3)
# geom_abline(slope=coef1[2], intercept = coef1[1], linetype=3)
p_mapped[[kk]] = p1
kk=kk+1
}
#### panel c & d
{
xmat = cbind(slopes_mat[,1:13]*sd_ragemat^pw,
1/SpeciesMat[,nj], slopes_mat[,14:15])
colnames(xmat)[14] = col2 = paste0("1/", colnames(SpeciesMat)[nj])
ncs = 1:15 ###c(k, grep("maxAgeCaesar",colnames(slopes_mat)))
species_xmat = combineTissues(xmat, idx1, ncs)
species_xmat$textidx = species_xmat$MammalNumberHorvath %in% textnum
## log scale panel
x = species_xmat[, 16]
y = species_xmat[,k+1]
# y = logplus(y)
# y = y - min(y) + 1
species_xmat$y = y ## log(y)
species_xmat$x = x ## log(x)
xlim = range(species_xmat$x)*c(0.95, 1.1)
cor1test = with(species_xmat, cor.test(x,y))
cor1 = signif(cor1test$estimate, 2)
pv1 = signif(cor1test$p.value,2)
coef1 = coef(lm(y~x, data = species_xmat))
mytitle = paste0(letters[kk],". ", len1, ", Cor=",cor1,
", p=",pv1) ##letters[kk],". ",
p1 = ggplot(species_xmat,
aes_string(y = "y", x = "x",
label = "MammalNumberHorvath")) +
geom_point(color="darkgrey") + theme_bw() +
scale_x_continuous(name='Lifespan', limits = xlim)+
scale_y_continuous(name=paste0('Adjusted AROCM')) +
ggtitle(mytitle) +
theme(plot.title=element_text(size = titlesize)) + # margin=margin(t=40,b=-30)
geom_hline(yintercept = 0, linetype = 3) +
geom_text(aes(label=ifelse(textidx,
MammalNumberHorvath,'')),
hjust=0,vjust=0)
p_mapped[[kk]] = p1
kk=kk+1
###
col1 = switch (kk/4,
"YoungSlope1",
"YoungSlope0.1",
"OldSlope1")
# species_xmat = species_xmat%>%filter(!is.na(species_xmat[,k+1]))
# colnames( species_xmat)
y = species_xmat[,k+1]
x = species_xmat[, 15]
lm1 = lm(y~ x-1)
coef1 = as.numeric(coef(lm1))
a = round(coef1[1],2)
species_xmat$x = species_xmat[,15]
species_xmat$y = y
cor1test = cor.test(x,y, use = 'pairwise.complete')
cor1 = signif(cor1test$estimate, 2)
pv1 = signif(cor1test$p.value,2)
mytitle = paste0(letters[kk],". ", len1,', N=',length(y)) #,", Cor=",cor1)
# species_xmat$x = x
# species_xmat$y = y
xr = range(c(x)) *c(1, 1.2)
#dlist[[kk]] = species_xmat
p1 = ggplot(species_xmat,
aes_string(y = "y", x = "x", # y = col1, x = col2,
label = "MammalNumberHorvath")) +
geom_point(color="darkgrey") + theme_bw() +
scale_x_continuous(name= paste0(col2), limits = xr) +
scale_y_continuous(name=paste0('Adjusted AROCM')
#, limits = xr
)+
ggtitle(mytitle, subtitle =
paste0("Cor=",cor1,", p=",pv1 )) +
theme(plot.title=element_text(size = titlesize),
plot.subtitle = element_text(size = 1.2*titlesize,
margin=margin(b=-15))) +
geom_text(aes(label=ifelse(textidx,
MammalNumberHorvath,'')),
hjust=0,vjust=0) +
geom_hline(yintercept = 0, linetype = 3)
p_mapped[[kk]] = p1
kk=kk+1
}
}
length(p_mapped)
p_mapped[[1]]
```
```{r}
vnum = 'v12'
# ggsave(paste0("Suppfig_",len1, "_pw=",pw, ".pdf"),
ggsave(paste0("Suppfig_AROCM_",len1, "_",vnum, ".pdf"),
grid.arrange(p_mapped[[1]],p_mapped[[3]],
p_mapped[[2]],p_mapped[[4]],
nrow=2, as.table = FALSE),
width = 16,
height = 16, units = "cm")
```
```{r}
vnum = 'v10'
# ggsave(paste0("Suppfig_",len1, "_pw=",pw, ".pdf"),
ggsave(paste0("Suppfig_AROCM_",len1, "_",vnum, ".pdf"),
grid.arrange(p_mapped[[1]],p_mapped[[2]],
nrow=1, as.table = FALSE),
width = 16,
height = 8, units = "cm")
```
Computing file changes ...
