#  San Román et al. A geometric morphometric protocol to correct postmorten body arching in fishes.
#  email: carla.sanroman@uam.es

#  Supplementary File 1: custom code
  

library(RRPP)
library(Matrix)
library(rgl)
library(geomorph)
library(Morpho)
library(ggplot2)


#Download Data: Supplementary File 2: OriginalFossilFishes

Raw <- readland.tps("FileName",specID = "ID" )
Fit.original<-procSym(Raw, use.lm =c(1:11,13,15:17))

##-  INDEX OF CURVATURE: IC  -##

lm<-matrix(c(13,25,13,20,20,21,21,22,22,23,23,24,24,25), ncol=2, byrow=TRUE,
           dimnames = list(c("standardlength","col1","col2","col3","col4","col5","col6"),c("start", "end")))
Fit.IC.geomorph<-gpagen(Raw)
lineardists <- interlmkdist(Fit.IC.geomorph$coords, lm)
dist<-as.data.frame(lineardists)
dist$curvelength<-rowSums(dist[,2:7])
indcurvature<-dist[,"standardlength"]/dist[,"curvelength"]

summary(1-indcurvature); hist(1-indcurvature) #Explore IC



  # ---------------------------- #
  #       DATA PREPARATION       #
  # ---------------------------- #

#---A. Original data

# Slab and Counter-slab AVERAGE
## DATA (without Lm of the backbone: Data: OriginalFossilFishesNOBackbone)

Fish.Nobackbone <- read.delim("FileName")
Fish.Nobackbone<-Fish.Nobackbone[order(Fish.Nobackbone$Slab), ]

duplicate<-c()
for (i in 1:dim(Fish.Nobackbone)[1]){
  if (Fish.Nobackbone[i,"Slab"]==Fish.Nobackbone[(i+1),"Slab"]){
    duplicate<-c(duplicate,i,(i+1))
  }
  
}
vslab<-c()
for (i in 1:dim(Fish.Nobackbone)[1]){
  if (Fish.Nobackbone[i,"Slab"]==Fish.Nobackbone[(i+1),"Slab"]) {
    v<-c(Fish.Nobackbone[i,"Slab"],Fish.Nobackbone[i,"Taxa"],as.numeric(mean(Fish.Nobackbone[i:(i+1),4])),as.numeric(colMeans(Fish.Nobackbone[i:(i+1),5:dim(Fish.Nobackbone)[2]],)))
    vslab<-rbind(vslab,v)
  }
}
Noduplicate<-Fish.Nobackbone[-duplicate,-1]
dfslab<-data.frame(vslab)
dfslab[,3:dim(dfslab)[2]] <- lapply(dfslab[,3:dim(dfslab)[2]], as.numeric)

colnames(dfslab)[1]<-"Slab"
colnames(dfslab)[2]<-"Taxa"
colnames(dfslab)[3]<-"IC"
colnames(dfslab)[4:dim(dfslab)[2]]<-paste("Coord",1:30,sep="")

df<-merge(Noduplicate,dfslab,all=T)
dfGordRub<-df[with(df, order(Taxa,Slab)), ] ##Final Dataframe
coordGordRub<-dfGordRub[,4:dim(dfGordRub)[2]]

coord.original.mean <- arrayspecs(coordGordRub, p = 15, k = 2)

fit.original.mean.geomorph<- gpagen(coord.original.mean)
fit.original.mean.geomorph.PC<-gm.prcomp(fit.original.mean.geomorph$coords)
fit.original.mean<- procSym(coord.original.mean)


#---B. Regression-unbending method

Reg.delete.bent <- procD.lm(f1=fit.original.mean.geomorph$coords ~ dfGordRub$IC ,logsz = FALSE,iter=999)

plot.GLS.curve.free<-plot(Reg.delete.bent, type="regression",reg.type="RegScore", predictor = dfGordRub$IC)
RegScores.curve.free<-plot.GLS.curve.free$RegScore

curve.free.residuals <- Reg.delete.bent$residuals
curve.free.residuals <- arrayspecs(curve.free.residuals,dim(fit.original.mean.geomorph$coords)[1], 2)
curve.free.shape <- curve.free.residuals +  array(fit.original.mean.geomorph$consensus, dim(curve.free.residuals))

coord.curvefree.final <- curve.free.shape

fit.unbendburnary.mean.geomorph<- gpagen(coord.curvefree.final)
fit.unbendburnary.mean.geomorph.PC<-gm.prcomp(fit.unbendburnary.mean.geomorph$coords)
fit.unbendburnary.mean<- procSym(coord.curvefree.final)


#---C. Tps-unbending method

# Split sample in groups
ind.curvature.factor<-cut(indcurvature,breaks=7,labels = c(1:7))
table(ind.curvature.factor)
splitFish<-data.frame(ID=rownames(Fit.original$PCscores),
                      greenlineFig1D=dist$standardlength,redlineFig1D=dist$curvelength,IC=indcurvature,IC.factor=ind.curvature.factor)
# write.csv(splitFish,file="Filename.txt")

### Conduct unbend function in tpsUtil program (see explanation in Material and Methods) --> Supplementary File 2: TpsFossilFishes

# Slab and counter-slab average
Tpsunbent <- read.delim("FileName")
Tpsunbent<-Tpsunbent[order(Tpsunbent$Slab), ]

duplicate.tps<-c()
for (i in 1:dim(Tpsunbent)[1]){
  if (Tpsunbent[i,"Slab"]==Tpsunbent[(i+1),"Slab"]){
    duplicate.tps<-c(duplicate.tps,i,(i+1))
  }
  
}
vslab.tps<-c()
for (i in 1:dim(Tpsunbent)[1]){
  if (Tpsunbent[i,"Slab"]==Tpsunbent[(i+1),"Slab"]) {
    v<-c(Tpsunbent[i,"Slab"],Tpsunbent[i,"Taxa"],as.numeric(colMeans(Tpsunbent[i:(i+1),4:dim(Tpsunbent)[2]],)))
    vslab.tps<-rbind(vslab.tps,v)
  }
}

Tpsnoduplicate<-Tpsunbent[-duplicate.tps,-1]
dfslab.tps<-data.frame(vslab.tps)
dfslab.tps[,3:dim(dfslab.tps)[2]] <- lapply(dfslab.tps[,3:dim(dfslab.tps)[2]], as.numeric) 

colnames(dfslab.tps)[1]<-"Slab"
colnames(dfslab.tps)[2]<-"Taxa"
colnames(dfslab.tps)[3:dim(dfslab.tps)[2]]<-paste("Coord",1:30,sep="")

dfslab.tps<-merge(Tpsunbent,dfslab.tps,all=T)
dfGordRubtps<-df[with(dfslab.tps, order(Taxa,Slab)), ] ##Final Dataframe 
coordGordRubtps<-dfGordRubtps[,3:dim(dfGordRubtps)[2]]

coord.unbenttps.mean <- arrayspecs(coordGordRubtps, p = 15, k = 2)


# Include IC (index of curvature)
a<-merge(dfGordRubtps,dfGordRub, by="Slab")
a<-a[with(a, order(Taxa.x,Slab)), ]
b<-a[,c(1:32,34)]

coordGordRubtps<-b[,3:(dim(b)[2]-1)]

coord.unbendtps.mean <- arrayspecs(coordGordRubtps, p = 15, k = 2)

fit.unbendtps.mean.geomorph<- gpagen(coord.unbendtps.mean)
fit.unbendtps.mean.geomorph.PC<-gm.prcomp(fit.unbendtps.mean.geomorph$coords)
fit.unbendtps.mean<- procSym(coord.unbendtps.mean)



  # --------------------------- #
  #           ANALYSIS          #
  # --------------------------- #

# Species. note: Gord (Gordichthys conquensis); Rub (Rubiesichthys gregalis)
colgordrub<-c(rep("#FE72A9",31),rep("#b9d852",63))
symbolgordrub<-c(rep(16,31),rep(15,63))


#---A. Original data

#--1. PCA
plot(fit.original.mean$PCscores[,1],fit.original.mean$PCscores[,2],main="PCA Original data",xlab="PC1", ylab="PC2"
     ,cex=1.2,axes=FALSE,col=colgordrub, pch=symbolgordrub)
axis(1);axis(2)
barplot(fit.original.mean$Variance[1:6,2],col="plum3",ylim=c(0,40))

# In "fit.original.mean.geomorph.PC$shapes$shapes.comp1$min" put "shapes.comp1" for seeing PC1 shape changes, "shapes.comp2" for seeing PC2 shape changes, etc.
deformedmin <- tps3d(fit.original.mean.geomorph$coords[,,17], fit.original.mean.geomorph$coords[,,17], fit.original.mean.geomorph.PC$shapes$shapes.comp1$min)
deformedmax <- tps3d(fit.original.mean.geomorph$coords[,,17], fit.original.mean.geomorph$coords[,,17], fit.original.mean.geomorph.PC$shapes$shapes.comp1$max)

deformGrid2d(fit.original.mean.geomorph$coords[,,17],deformedmin,ngrid =0,lines=F,lwd=2,lcol="black",col1="white",col2="black",wireframe=c(c(1:8,10:15,10,15,14,1)), margin=T, cex1=1.5)
deformGrid2d(fit.original.mean.geomorph$coords[,,17],deformedmax,ngrid =0,lines=F,lwd=2,lcol="black",col1="white",col2="black",wireframe=c(c(1:8,10:15,10,15,14,1)), margin=T, cex1=1.5)


#--2. SHAPE DATA ~ IC
Reg.bent <- procD.lm(f1=fit.original.mean.geomorph$coords ~ dfGordRub$IC ,logsz = FALSE,iter=999)
summary(Reg.bent)

plot(Reg.bent,type="regression",predictor=1-dfGordRub$IC, reg.type="RegScore",
     pch = 16, xlab = "IC", ylab="shape score",cex=1.5,col=colgordrub,axes=F)
axis(1);axis(2)


#--3. PC~ IC
Reg.bend.PC <- lm(fit.original.mean$PCscores[,1] ~ dfGordRub$IC)  #Change "PCscores[,1]" to see other PCs
summary(Reg.bend.PC)


#--4. Differences IC ~ Species (Gordichthys conquensis and Rubiesichthys gregalis)
shapiro.test(dfGordRub$IC) #Normality test
wilcox.test(dfGordRub$IC[1:31],dfGordRub$IC[32:63])



#---B. Regression- Unbending Data

#--1. PCA
plot(fit.unbendburnary.mean$PCscores[,1],fit.unbendburnary.mean$PCscores[,2],main="PCA Regression-Unbending Method",xlab="PC1", ylab="PC2" 
     ,cex=1.2,axes=FALSE,col=colgordrub, pch=symbolgordrub)
axis(1);axis(2)
barplot(fit.unbendburnary.mean$Variance[1:6,2],col="plum3",ylim=c(0,40))

#change "shapes.comp1" to see others PCs
deformedmin <- tps3d(fit.unbendburnary.mean.geomorph$coords[,,17], fit.unbendburnary.mean.geomorph$coords[,,17], fit.unbendburnary.mean.geomorph.PC$shapes$shapes.comp2$min)
deformedmax <- tps3d(fit.unbendburnary.mean.geomorph$coords[,,17], fit.unbendburnary.mean.geomorph$coords[,,17], fit.unbendburnary.mean.geomorph.PC$shapes$shapes.comp2$max)
deformGrid2d(fit.unbendburnary.mean.geomorph$coords[,,17],deformedmin,ngrid =0,lines=F,lwd=2,lcol="black",col1="white",col2="black",wireframe=c(c(1:8,10:15,10,15,14,1)), margin=T, cex1=1.5)
deformGrid2d(fit.unbendburnary.mean.geomorph$coords[,,17],deformedmax,ngrid =0,lines=F,lwd=2,lcol="black",col1="white",col2="black",wireframe=c(c(1:8,10:15,10,15,14,1)), margin=T, cex1=1.5)
fit.unbendburnary.mean$Variance



#---C. Tps- Unbending Data

#--1. PCA
fit.unbendtps.mean$Variance

plot(fit.unbendtps.mean$PCscores[,1],fit.unbendtps.mean$PCscores[,2],main="PCA Tps-Unbending Method",xlab="PC1", ylab="PC2" 
     ,cex=1.2,axes=FALSE,col=colgordrub, pch=symbolgordrub)
axis(1);axis(2)
barplot(fit.unbendtps.mean$Variance[1:6,2],col="plum3",ylim=c(0,40))

deformedmin <- tps3d(fit.unbendtps.mean.geomorph$coords[,,17], fit.unbendtps.mean.geomorph$coords[,,17], fit.unbendtps.mean.geomorph.PC$shapes$shapes.comp2$min)
deformedmax <- tps3d(fit.unbendtps.mean.geomorph$coords[,,17], fit.unbendtps.mean.geomorph$coords[,,17], fit.unbendtps.mean.geomorph.PC$shapes$shapes.comp2$max)
deformGrid2d(fit.unbendtps.mean.geomorph$coords[,,17],deformedmin,ngrid =0,lines=F,lwd=2,lcol="black",col1="white",col2="black",wireframe=c(c(1:8,10:15,10,15,14,1)), margin=T, cex1=1.5)
deformGrid2d(fit.unbendtps.mean.geomorph$coords[,,17],deformedmax,ngrid =0,lines=F,lwd=2,lcol="black",col1="white",col2="black",wireframe=c(c(1:8,10:15,10,15,14,1)), margin=T, cex1=1.5)





  # ----------------------------- #
  #   Testing Unbending methods   #
  # ----------------------------- #

#-1. Disparity. Calcule the disparity of each group and then compare between groups

##DATA. 
 #We need to merge data and yield a single common GPA. Then use the coords resulting for subsequent analysis

#BENT sample
write.csv(coordGordRub1, 'bentcoords.csv') 
#UNBENT by TPS sample
write.csv(coordGordRubtps, 'tpscoords.csv') 
#UNBENT by regression sample
regression<-two.d.array(coord.curvefree.final)#transform array into "matrix"
write.csv(regression, 'regcoords.csv')

#Combine csv and import the combined csv in excel
#Data: BENT_TPS_REG-coords.csv

benttpscoords <- read.csv("C:/Users/Carla San Román/OneDrive - UAM/TESIS/Siensia/Unbending/data/final/benttpsregcoords.csv", sep=";")
bent.tps.reg <- arrayspecs(benttpscoords[2:dim(benttpscoords)[2]], p = 15, k = 2)
benttpsgpagen<-gpagen(bent.tps.reg)
method<-c(rep("original",94),rep("tps",94),rep("reg",94))#GROUPS: original- tps- regression data

##Dispary analysis. Calculate disparity in each sample and contrast (pairwise test) differences between samples. 
#Also gives procrustes variance of each sample

disparity.unbending <- morphol.disparity(benttpsgpagen$coords ~ method, groups = ~ method, iter=999)
summary(disparity.unbending)
disparity.unbending$Procrustes.var




#-2. Procrustes variance and variance per LM

library(tibble)

#-- Procrustes coordinates of each sample

per_lm_variance1 <- function (shape.data) #para modificar los colores de los LMs, meter los codigos que gusten en cols1
{
  variances <- rowSums(apply(shape.data, c(1, 2), var))
  cols1 <- colorRampPalette(c("gray88", "yellow", "red"))
  cols <- cols1(100)
  x = (log10(variances))
  xlims <- NULL
  tol <- 1e-06
  xlims <- range(x) + c(-tol, tol)
  nbin = 100
  breaks <- 0:nbin/nbin * (xlims[2] - xlims[1]) + xlims[1]
  whichColor <- function(p, cols, breaks) {
    i <- 1
    while (p >= breaks[i] && p > breaks[i + 1]) i <- i +
        1
    cols[i]
  }
  variancecolors <- sapply(x, whichColor, cols = cols, breaks = breaks)
  variance_table <- tibble(Per_Lm_Variance = variances, Log_Variance = x,Variance_Colors = variancecolors)
  return(variance_table)
}

original.variance<-per_lm_variance1(fit.original.mean.geomorph$coords)
tps.variance<-per_lm_variance1(fit.unbendtps.mean.geomorph$coords)
burnary.variance<-per_lm_variance1(fit.unbendburnary.mean.geomorph$coords)


#Color gradient is useful for representing LM variance in deformGrid2d function. We can use the mean

deformGrid2d(fit.unbendtps.mean.geomorph$coords[,,17],mshape(fit.unbendtps.mean.geomorph$coords),ngrid =0,lines=F,lwd=2,
             lcol="black",col1="white",col2=tps.variance$Variance_Colors,wireframe=c(c(1:8,10:15,10,15,14,1)), margin=T, cex2=1.5)