Additional file 5

logo # Analysis of functional assignment (COG ID) from metagenomic data The process to analyse a COG database takes its compositional nature into account was:

Filtering data
ALR transformation
Procrustes analysis
Identification outliers samples
Partial least square-Discriminant analysis (PLS-DA)
Classification performance of the model

Loading libraries and datasets

library(pacman)
pacman::p_load(readxl,data.table,easyCODA,compositions,mixOmics,factoextra,dplyr)

Filtering data

#SqueezeMeta dataset for COG abundance; filter_merge17754.COG.abund.tsv

count <- fread("sqm/filter_merge17754.COG.abund.tsv")
count[1,1] <- "COG_ID"
colnames(count) <- as.character(count[1,])
count <- count[-1,]

#Statistics of reads per sample and percentage of mapped reads; 10.project_name.mappingstat

stat <- read.table("C:/Users/3ccas/OneDrive - UPV/2020_2021/Metagenomica/10.filter_merge17754.mappingstat",header=T,stringsAsFactors = F)

#Samples description

type <- read_xlsx("C:/Users/3ccas/OneDrive - UPV/2020_2021/Metabolomica/ANALISIS SANGRE_Supervisado_Cristina.xlsx")
type <- type[type$HEMBRA%in%stat$Sample,]

#Assignment the Line/Population of each sample to the "stat" variable

stat$Line <- type$LINEA[match(stat$Sample,type$HEMBRA)]
stat$Line[stat$Line=="A"] <- "Low"
stat$Line[stat$Line=="V"] <- "High"

#Dataset rearrangement 

ID.var <- count$COG_ID
count <- count[,-1]
ID.samples <- colnames(count)
count.t <- data.frame(t(count))
colnames(count.t) <- ID.var
rownames(count.t) <- ID.samples

line <- stat$Line[match(rownames(count.t),stat$Sample)]

#Remove all variables with a percentage of zeros higher than 20% within-population
n.l <- round(length(which(line=="Low"))*0.2)
n.h <- round(length(which(line=="High"))*0.2)

zero.col<-data.frame(table(which(count.t==0,arr.ind = T)[,2]))
zero.col<-as.numeric(as.character(zero.col$Var1[zero.col$Freq>(n.l+n.h)]))

index.zero<-NULL;n<-0
for(i in zero.col){
  nh<-length(which(count.t[line=="High",i]==0))
  nl<-length(which(count.t[line=="Low",i]==0))
  if(nh/nl>0.5 | nh/nl<2){
    n<-n+1
    index.zero[n]<-i
  }
}

indx.l <- which(count.t[line=="Low",]==0,arr.ind = T)
indx.h <- which(count.t[line=="High",]==0,arr.ind = T)
zeros.l <- data.frame(table(indx.l[,2]))
zeros.h <- data.frame(table(indx.h[,2]))

indx.l <- as.numeric(as.character(zeros.l$Var1[which(zeros.l$Freq>n.l)]))
indx.h <- as.numeric(as.character(zeros.h$Var1[which(zeros.h$Freq>n.h)]))

col<-unique(c(indx.h[indx.h%in%indx.l],index.zero))
count.zero <- count.t[,-col]
count.zero <- count.zero+1

write.table(count.t,"cog_raw.txt",col.names=T,row.names=F,quote=F,sep=";")
write.table(count.zero,"cog_raw0.txt",col.names=T,row.names=F,quote=F,sep=";")

Additive log-ratio (ALR) transformation

Greenacre M, Martínez-Álvaro M, Blasco A. Compositional data analysis of microbiome and any-omics datasets: a revalidation of the additive logratio transformation. bioRxiv. 2021. 10.1101/2021.05.15.444300

DOI: https://doi.org/10.1101/2021.05.15.444300

RecA

RecA (recombination protein RecA) was used as a reference variable for Kegg (K03553). This protein had been suggested as a reference variable because is present in most of the Bacteria, Archaea and Eukaryotes organisms. Moreover, it has a low copy number variation between taxa. It could be like a copy cell number estimation.

Reference:

Wu D, Jospin G, Eisen JA. Systematic identification of gene families for use as “markers” for phylogenetic and phylogeny-driven ecological studies of bacteria and Archaea and their major subgroups. PLoS One. 2013;8(10):e77033. Published 2013 Oct 17.

DOI:10.1371/journal.pone.0077033

#RecA location

indx <- grep("COG0468",names(count.zero))
recA <- log(count.zero/count.zero[,indx])

Procrustes analysis

Procrustes analysis was performed to test if the selected reference variable allowed how close they come to the exact geometry. This means that the variables maintain their relationship despite the transformation.

Procrustes analysis is based on three simple operations:

Centering
Scaling
Rotation

#Calculation of LRA (isometric matrix) and ALR matrix

LRA <- LRA(count.zero,weight = F)$rowpcoord

ALR <- ALR(count.zero,denom = indx)$LR
ALR.PCA <- PCA(ALR,weight=F)$rowpcoord

protest(LRA,ALR.PCA[,1:ncol(LRA)], permutations=0)$t0

## [1] 0.9978047

Identification of outlier samples

To identify the outlier samples a principal component analysis was performed. After the removal of animals, the centering and scaling steps must perform again.

#Autoscaling

recA.scale <- recA

for(i in 1:ncol(recA)){
  
  recA.scale[,i] <- recA[,i] - mean(recA[,i])
  recA.scale[,i] <- recA.scale[,i] / sd(recA.scale[,i])

  }

recA.scale <- recA.scale[,-indx]

#Line/Population assignment

recA.scale$Line <- stat$Line[match(rownames(recA.scale),stat$Sample)]

#Identified outlier samples

out <- grep("18070|17763|18041|17661|17901|17722",row.names(count.zero))
recA.out <- recA[-out,]

#Centering and scaling

recA.out.scale <- recA.out

for(i in 1:ncol(recA.out)){
  
  recA.out.scale[,i] <- recA.out[,i] - mean(recA.out[,i])
  recA.out.scale[,i] <- recA.out.scale[,i] / sd(recA.out.scale[,i])

  }

recA.out.scale <- recA.out.scale[,-indx]
recA.out.scale$Line <- stat$Line[match(rownames(recA.out.scale),stat$Sample)]

We removed only the most influential animals according to the additional file 4. After that, the autoscaling was computed again.

color <- c("#31A2AC","#AF1C1C","#2F2F28","#F0EFF0")
pca.alr <- prcomp(recA.out.scale[,-ncol(recA.out.scale)],scale=F)

#Write final dataset

write.table(recA.out.scale,"COG.ALR.txt",col.names=T,row.names=F,quote=F,sep=";")

#Write PCA after removal outlier samples

p <- fviz_pca_ind(pca.alr,geom=c("point","text"),col.ind=recA.out.scale$Line,axes=c(1,2),addEllipses = TRUE,ellipse.level = 0.95,pointsize=1) +
  geom_point(aes(fill=recA.out.scale$Line,color=recA.out.scale$Line,shape=recA.out.scale$Line))+
  scale_shape_manual(values=c(16,17))+
  scale_color_manual(values =color)+
  scale_fill_manual(values =color) +
  xlim(-150, 150) + ylim (-150, 150)+
  geom_hline(yintercept = 0, colour="#2F2F28", linetype="dashed") + 
  geom_vline(xintercept = 0, colour="#2F2F28", linetype="dashed") + 
  labs(x ="PC1 (18.7%)", y = "PC2 (9.7%)",title = NULL,face="bold",fill="Population",shape="Population",color="Population") + 
  theme_minimal() +
  theme_classic()+
  theme( 
    legend.position="bottom",
    panel.border = element_blank(),
    panel.grid.major.x = element_blank(),
    panel.grid.minor.x = element_blank(),
    axis.text.x = element_text( size = 9, vjust = 1.5),
    axis.text.y = element_text( size = 9, vjust = 0.7))
p

ggsave("C:/Users/3ccas/OneDrive - UPV/2020_2021/Metagenomica/Plot/pca.COG_initial.tiff",p, width = 120, height = 120, units = "mm",dpi = 600)

Partial Least Square-Discriminant Analysis (PLS-DA)

The PLS-DA tries to extract the latent structures (patterns) that allow explaining a dependent variable (Y; response). With this analysis, we can reduce the dimension of the data and collect the information throughout the maximization of the covariance between the X and Y. In this case, X is the ALR matrix with all KEGG and Y is a vector with the Line/Population of each sample.

#Computing the PLS-DA model with 10 components

plsda <- plsda(recA.out.scale[,-ncol(recA.out.scale)],recA.out.scale$Line,ncomp=10,scale=F)

#4-fold cross-validation to compute the balance error rate (BER)

set.seed(30)
perf.pls <- perf(plsda, validation = "Mfold",criterion="all",folds = 4,
                 progressBar = F,nrepeat = 100,)

#Number of components with the minor balance error rate

ber <- 1
ber.n <- 0.5
c <- 20

while(ber.n < ber) {
  
  c = c+1
  ber.n <- perf.pls[["error.rate"]]$BER[c+1]
  ber <- perf.pls[["error.rate"]]$BER[c]
  
  }

comp <- which(perf.pls[["error.rate"]]$BER[21:30]==ber)

#BER for the number of components selected

err.n <- perf.pls[["error.rate"]]$BER[20+comp]
sd.n <- perf.pls[["error.rate.sd"]]$BER[20+comp]
err.total <- err.n + sd.n

#Variable important prediction (VIP). Contribution of each variable in the classification among lines/populations

vip <- data.frame(vip(plsda),stringsAsFactors = FALSE)

#A VIP higher than 1 was used as the threshold for selecting the variables with the highest contribution in the model

p <- 1
v.select <- vip[vip[,comp]>=p,]
v.ID <- row.names(v.select)

The next PLS-DA are an iterative process until the BER reaches the minimum value.

err <- 1
comp.cte <- 10
while (err.total<err | err.total<0.02) {
  
  err <- err.total
  sd.f <- sd.n
  comp.f <- comp
  
  vf<-v.ID
  vip <- data.frame(vip(plsda),stringsAsFactors = FALSE)
  v.select <- vip[vip[,comp] >= p,]
  v.ID <- row.names(v.select)
  filter <- recA.out.scale[,names(recA.out.scale) %in% c(v.ID,"Line")]
  
   if (length(v.ID)<comp.cte) {
    
    comp.cte <- length(v.ID)
  
  }
  
  plsda <- plsda(filter[,-ncol(filter)],filter$Line,ncomp=10,scale=F)
  set.seed(30)
  perf.pls <- perf(plsda, validation = "Mfold",criterion="all",folds = 4,
                   progressBar = F,nrepeat = 100,)

  ber <- 1
  ber.n <- 0.5
  c <- comp.cte*2
  while (ber.n<ber) {
    
    c = c + 1
    ber.n <- perf.pls[["error.rate"]]$BER[c+1]
    ber <- perf.pls[["error.rate"]]$BER[c]
  
    }
  
  comp <- which(perf.pls[["error.rate"]]$BER[(comp.cte*2+1):(comp.cte*3)]==ber)[1]
  err.n <- perf.pls[["error.rate"]]$BER[(comp.cte*2)+comp]
  sd.n <- perf.pls[["error.rate.sd"]]$BER[(comp.cte*2)+comp]
  err.total <- err.n + sd.n

  }

In summary, the model with the minimum BER (Balance error rate) for the classification will be the following specifications:

#Optimal number of components in the model

comp.f

## [1] 3

#Balanced error rate of the model

err - sd.f

## [1] 0.0003571429

#Standard deviation of the BER

sd.f

## [1] 0.002512595

#Number of variables included in the model

length(vf)

## [1] 164

#Relevant variables for the classification

write.table(names(recA.out[grep(paste(vf,collapse = "|"),names(recA.out))]),"C:/Users/3ccas/OneDrive - UPV/2020_2021/Metagenomica/Plot/relevant_COG_ID.txt",row.names = T,col.names = T,quote = F,sep="\t")

#Final model 

filter <- recA.out.scale[,names(recA.out.scale)%in%c(vf,"Line")]
plsda <- plsda(filter[,-ncol(filter)],filter$Line,ncomp=comp.f,scale=F)

#PLS-Plot of the final model

plotIndiv(plsda,ind.names = TRUE, legend=TRUE,style = "ggplot2",rep.space = "X-variate",
             ellipse = TRUE, centroid=TRUE,title = 'PLS-DA on Line',
             X.label = 'Comp 1', Y.label = 'Comp 2',col = color[1:2],abline=TRUE,background = NULL,cex = c(5,5),
             legend.title = "Line")

#PCA-Plot of final model

pca <- prcomp(filter[,-ncol(filter)],scale=F)
fviz_pca_ind(pca,axes=c(1,2),geom = c("point","text"),col.ind=filter$Line,addEllipses = T,palette=color,ellipse.level = 0.95,pointsize = 1)

vip <- data.frame(vip(plsda),stringsAsFactors = FALSE)

write.table(vip[,1:comp.f],"C:/Users/3ccas/OneDrive - UPV/2020_2021/Metagenomica/Plot/VIP_COG.txt",sep = "\t",quote=F)
write.table(filter,"C:/Users/3ccas/OneDrive - UPV/2020_2021/Metagenomica/Plot/relevant_COG_total.txt",row.names = T,col.names = T,quote = F,sep="\t")

Quality of the model

To check the quality of the model, two tests were performed using a 4-fold Cross-validation 10,000 times specifying the Mahalanobis distance

Confusion matrix

Allow knowing the success rate for the prediction of each rabbit line/population. Percentage of false positive and false negative.

library(caTools)

cross.val <- filter
confusion.total <- matrix(ncol=2,nrow=2,0)
x.total <- NULL

for (i in 1:10000) {
  
  data.RF <- cross.val
  
  sample = sample.split(data.RF$Line, SplitRatio = .70)
  train = subset(data.RF, sample == TRUE)
  test  = subset(data.RF, sample == FALSE)
  dim(train)
  dim(test)
  
  x <- data.frame(table(test$Line))
  x.total <- rbind(x,x.total)
  
  plsda.train <- plsda(train[,-ncol(train)],train$Line,ncomp=comp.f,scale=F)
  test.predict <- predict(plsda.train,test[,-ncol(test)],dist = "mahalanobis.dist")
  prediction <- test.predict$class$mahalanobis.dist[,comp.f]
  
  confusion.mat <- get.confusion_matrix(truth = test$Line,predicted =prediction )
  confusion.total <- confusion.total+confusion.mat

  }

High <- sum(x.total$Freq[x.total$Var1=="High"])
Low <- sum(x.total$Freq[x.total$Var1=="Low"])

confusion.total[1,]<- 100*confusion.total[1,]/High
confusion.total[2,] <- 100*confusion.total[2,]/Low

confusion.total

##      predicted.as.High predicted.as.Low
## High           99.8775           0.1225
## Low             0.0040          99.9960

Permutation matrix

Allow computing the spurious rate for the randomly prediction of each rabbit line/population

library(gtools)

cross.val <- filter
permutation.total <- matrix(ncol=2,nrow=2,0)
x.total <- NULL

for (i in 1:10000) {
  
  data.RF$Line <- permute(cross.val$Line)
  
  sample = sample.split(data.RF$Line, SplitRatio = .70)
  train = subset(data.RF, sample == TRUE)
  test  = subset(data.RF, sample == FALSE)
  dim(train)
  dim(test)
  
  x <- data.frame(table(test$Line))
  x.total <- rbind(x,x.total)
  
  plsda.train <- plsda(train[,-ncol(train)],train$Line,ncomp=comp.f,scale=F)
  test.predict <- predict(plsda.train,test[,-ncol(test)],dist = "mahalanobis.dist")
  prediction <- test.predict$class$mahalanobis.dist[,comp.f]
  
  confusion.mat <- get.confusion_matrix(truth = test$Line,predicted =prediction )
  permutation.total <- permutation.total+confusion.mat
}
High<-sum(x.total$Freq[x.total$Var1=="High"])
Low<-sum(x.total$Freq[x.total$Var1=="Low"])

permutation.total[1,]<-100*permutation.total[1,]/High
permutation.total[2,]<-100*permutation.total[2,]/Low

permutation.total

##      predicted.as.High predicted.as.Low
## High            45.790           54.210
## Low             45.799           54.201