Code from: Phenotypic plasticity in visual opsin gene expression: a meta-analysis in teleost fish
César Bertinetti
January 23 2025
The datasets needed to reproduce this analysis are deposited on Bertinetti et al. 2025. Data from: Phenotypic plasticity in visual opsin gene expression: a meta-analysis in teleost fish [Dataset]
Following R version and packages were used to run this code:
R.version$platform## [1] "aarch64-apple-darwin20"R.version$version.string## [1] "R version 4.4.1 (2024-06-14)"installed.packages()[names(sessionInfo()$otherPkgs), "Version"]## rotl geiger phytools maps phangorn ggdist
## "3.1.0" "2.0.11" "2.3-0" "3.4.2" "2.12.1" "3.3.2"
## metafor numDeriv metadat ape phylosignal phylobase
## "4.7-58" "2016.8-1.1" "1.2-0" "5.8-1" "1.3.1" "0.8.12"
## visreg bbmle MASS partR2 ggpubr rsq
## "2.7.0" "1.0.25.1" "7.3-61" "0.9.2" "0.6.0" "2.6"
## lme4 Matrix gridExtra ggplot2 nlme agricolae
## "1.1-35.5" "1.7-0" "2.3" "3.5.1" "3.1-166" "1.3-7"
## corrplot Hmisc car carData zoo MESS
## "0.95" "5.2-1" "3.1-3" "3.0-5" "1.8-12" "0.5.12"
## RColorBrewer plotrix tidyr dplyr matrixStats stringr
## "1.1-3" "3.8-4" "1.3.1" "1.1.4" "1.4.1" "1.5.1"
## readr
## "2.1.5"## HERE USE setwd() to set your working directory where the data files will be locatedOpsin Gene Expression
genes <- c("sws1","sws2","rh1","rh2","lws") # gene names
## READ OPSIN DATA WITH Mean, standard deviation and samples from each study
opsin <- read.table("Opsin_Data.csv", header=T, sep =",", dec=".",stringsAsFactors = T)
## READ SPECIES ID INCLUDED TO GENERATE PHYLOTREE
taxon <- read.table("Species_ID.csv", header=T, sep =",", dec=".",stringsAsFactors = F,quote="\"")
### CHECK NUMBER OF SPECIES
length(unique(taxon$Species))## [1] 37
Estimate standardized mean differences based on: Viechtbauer, W.
(2010). Conducting meta-analyses in R with the metafor package. Journal
of Statistical Software, 36(3), 1–48. https://doi.org/10.18637/jss.v036.i03
opsin_meta <- escalc(measure="SMD", m1i=Mean1, sd1i=SD1, n1i=N1,
m2i=Mean2, sd2i=SD2, n2i=N2, di=Dval, data=opsin,slab = Author)
## READ GENERATED PHYLOGENETIC TREE
tree <- read.tree("plast.tre")
### Transform tree into an ultrametric one
tree <- compute.brlen(tree)
### compute phylogenetic correlation matrix
A <- vcv(tree, corr=TRUE)
### MERGED PHYLOGENTIC TREE AND OPSIN PLASTICITY DATA
opsin_master <- merge(opsin_meta, taxon[,-c(2:3,5)], by="Species",all=T)
opsin_master$search_string <- sub(" ", "_", opsin_master$search_string)
colnames(opsin_master)[23] <- "Tree_ID"
# sum(is.na(opsin_master$Tree_ID)) ### Check no species IDs are mising values!
### USE ABSOLUTE SMD AS A MEASURE OF PLASTICITY REGARDLESS OF DIRECTION OF CHANGE
opsin_master$yi <- abs(opsin_master$yi)
Run mixed-effect meta-analysis models including moderators and
random effects.
model = ~ Treatment * Type
## model = ~ 0 + Treatment:Type can be used to leave out model intercept
## the following model takes a while to run
## multi-level phylogenetic meta-analytic model
overall_int <- rma.mv(yi, vi,
mods = model,
random = list(~ 1 | Study/Sample, ~ 1 | Species, ~ 1 | Tree_ID),
R=list(Tree_ID=A), data=opsin_master, control=list(rel.tol=1e-8))
overall_int##
## Multivariate Meta-Analysis Model (k = 573; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.1855 0.4307 36 no Study no
## sigma^2.2 0.0000 0.0004 63 no Study/Sample no
## sigma^2.3 0.0412 0.2029 37 no Species no
## sigma^2.4 0.0000 0.0000 37 no Tree_ID yes
##
## Test for Residual Heterogeneity:
## QE(df = 569) = 1570.6069, p-val < .0001
##
## Test of Moderators (coefficients 2:4):
## QM(df = 3) = 25.7791, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## intrcpt 1.8197 0.2147 8.4748 <.0001 1.3989 2.2406
## Treatmentlight -1.0663 0.2207 -4.8308 <.0001 -1.4990 -0.6337
## Typereared -0.4550 0.2217 -2.0521 0.0402 -0.8896 -0.0204
## Treatmentlight:Typereared 0.5694 0.2314 2.4601 0.0139 0.1158 1.0230
##
## intrcpt ***
## Treatmentlight ***
## Typereared *
## Treatmentlight:Typereared *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(overall_int, btt=3:4) ## test interaction among moderators##
## Test of Moderators (coefficients 3:4):
## QM(df = 2) = 6.9721, p-val = 0.0306
This can now be redone to test heterogeneity depending on
experimental treatments or moderators. Requires creating smaller
phylogenetic trees for eachd dataset.
## LIGHT referred as External Stimuli
kp_sp <- unique(subset(opsin_master, Treatment=="light")$Tree_ID)
over_phy <- vcv(compute.brlen(drop.tip(tree, setdiff(tree$tip.label, kp_sp))), corr=T)
## Run model, takes some time
overall_l <- rma.mv(yi, vi, subset=(Treatment=="light"),
random = list(~ 1 | Study/Sample, ~ 1 | Species, ~ 1 | Tree_ID),
R=list(Tree_ID=over_phy), data=opsin_master, control=list(rel.tol=1e-8),slab = paste(Author, Sample, Age_dph, sep=", "))
overall_l##
## Multivariate Meta-Analysis Model (k = 489; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.1916 0.4377 31 no Study no
## sigma^2.2 0.0051 0.0716 53 no Study/Sample no
## sigma^2.3 0.0333 0.1826 34 no Species no
## sigma^2.4 0.0000 0.0009 34 no Tree_ID yes
##
## Test for Heterogeneity:
## Q(df = 488) = 1351.7993, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 0.8389 0.0961 8.7273 <.0001 0.6505 1.0273 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
This can be expanded to test moderator (e.g, light spectra)
effect on each opsin.
mod_light_opsin = ~ 0 + Type:Opsin_Class
overall_l_ops_mod <- rma.mv(yi, vi, subset=(Treatment=="light"), mods=mod_light_opsin,
random = list(~ 1 | Study/Sample, ~ 1 | Species, ~ 1 | Tree_ID),
R=list(Tree_ID=over_phy), data=opsin_master, control=list(rel.tol=1e-8),slab = paste(Author, Sample, Age_dph, sep=", "))
overall_l_ops_mod##
## Multivariate Meta-Analysis Model (k = 489; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.1911 0.4371 31 no Study no
## sigma^2.2 0.0000 0.0003 53 no Study/Sample no
## sigma^2.3 0.0407 0.2018 34 no Species no
## sigma^2.4 0.0000 0.0002 34 no Tree_ID yes
##
## Test for Residual Heterogeneity:
## QE(df = 479) = 1281.4666, p-val < .0001
##
## Test of Moderators (coefficients 1:10):
## QM(df = 10) = 124.9478, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## Typeinduced:Opsin_Classlws 0.7379 0.1223 6.0347 <.0001 0.4982 0.9776
## Typereared:Opsin_Classlws 0.8205 0.1072 7.6506 <.0001 0.6103 1.0307
## Typeinduced:Opsin_Classrh1 0.5634 0.5247 1.0737 0.2830 -0.4651 1.5918
## Typereared:Opsin_Classrh1 0.4367 0.1792 2.4363 0.0148 0.0854 0.7880
## Typeinduced:Opsin_Classrh2 0.7883 0.1274 6.1877 <.0001 0.5386 1.0381
## Typereared:Opsin_Classrh2 0.8429 0.1053 8.0076 <.0001 0.6366 1.0492
## Typeinduced:Opsin_Classsws1 0.7856 0.1415 5.5504 <.0001 0.5082 1.0630
## Typereared:Opsin_Classsws1 1.2816 0.1187 10.8004 <.0001 1.0490 1.5141
## Typeinduced:Opsin_Classsws2 0.7341 0.1281 5.7320 <.0001 0.4831 0.9852
## Typereared:Opsin_Classsws2 0.7946 0.1080 7.3604 <.0001 0.5830 1.0062
##
## Typeinduced:Opsin_Classlws ***
## Typereared:Opsin_Classlws ***
## Typeinduced:Opsin_Classrh1
## Typereared:Opsin_Classrh1 *
## Typeinduced:Opsin_Classrh2 ***
## Typereared:Opsin_Classrh2 ***
## Typeinduced:Opsin_Classsws1 ***
## Typereared:Opsin_Classsws1 ***
## Typeinduced:Opsin_Classsws2 ***
## Typereared:Opsin_Classsws2 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(overall_l_ops_mod, X=c(1,-1,0,0,0,0,0,0,0,0)) # for example for lws red-sensitive opsin##
## Hypothesis:
## 1: Typeinduced:Opsin_Classlws - Typereared:Opsin_Classlws = 0
##
## Results:
## estimate se zval pval
## 1: -0.0826 0.0922 -0.8962 0.3701
##
## Test of Hypothesis:
## QM(df = 1) = 0.8032, p-val = 0.3701
Or this can also be used to test overall effect sized on each
opsin
sws1 <- subset(opsin_master, Opsin_Class == "sws1") # Subset opsin of interest
## Drop unneeded species by selecting the ones present in the subset
kp_sp <- unique(subset(opsin_master, Opsin_Class=="sws1")$Tree_ID)
sws1_phy <- vcv(compute.brlen(drop.tip(tree, setdiff(tree$tip.label, kp_sp))), corr=T)
## Run model - Overall
res.sws1 <- rma.mv(yi, vi,
random = list(~ 1 | Study/Sample, ~ 1 | Species, ~ 1 | Tree_ID),
R=list(Tree_ID=sws1_phy), data=sws1, control=list(rel.tol=1e-8),slab = paste(Author, Sample, Age_dph, sep=", "))
res.sws1##
## Multivariate Meta-Analysis Model (k = 88; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.2386 0.4885 28 no Study no
## sigma^2.2 0.0777 0.2787 48 no Study/Sample no
## sigma^2.3 0.2897 0.5382 29 no Species no
## sigma^2.4 0.0000 0.0010 29 no Tree_ID yes
##
## Test for Heterogeneity:
## Q(df = 87) = 383.1683, p-val < .0001
##
## Model Results:
##
## estimate se zval pval ci.lb ci.ub
## 1.2854 0.1755 7.3250 <.0001 0.9414 1.6293 ***
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1PHYLOGENETIC SIGNAL
This mostly based on the phylosignal package:
Keck F., Rimet F., Bouchez A. & Franc A. (2016) phylosignal: an R package to measure, test, and explore the phylogenetic signal. Ecology and Evolution 6 (9), 2774–2780. http://onlinelibrary.wiley.com/doi/10.1002/ece3.2051/epdf
## JUST A BASIC FUNCTION
momdis = function(x) {
c(m = mean(x), s = sd(x), iqr = IQR(x), n = length(x), se = std.error(x), max = max(x))
}
### SUMMARISE PLASTICITY VALUES FOR EACH SPECIES
mean_plast <- aggregate(yi ~ Species + Tree_ID, data = opsin_master, momdis)
mean_plast <- do.call(data.frame, mean_plast)
p4 <- data.frame(mean_plast[,3])
row.names(p4) <- mean_plast$Tree_ID
## COMBINE TREE AND MEAN PLASTICITY
p4d <- phylo4d(tree, p4)
#Assessing the signal depth with correlograms
plast.crlg <- phyloCorrelogram(p4d, trait = "mean_plast...3.")
plot(plast.crlg)