Overview

Here we present a reproducible protocol for analysing bacterial proteomes developed during the PhD research of Camila P. Perico and described in our study currently in preparation.

The workflow consists of:

obtaining bacterial proteomes;
identifying an Evolutionarily-relevant Protein Core (EPC);
vectorizing protein sequences with rSWeeP;
performing dimensionality reduction with Principal Component Analysis (PCA);
calculating pairwise Pearson correlation distances and other distances and reconstructing phylogenetic trees; and
Visualizing Phylogenetic Trees

This tutorial provides a fully reproducible implementation in R. To reproduce the study results, we have made the complete data available on the Zenodo platform via this link [DOI: 10.5281/zenodo.13999354]. This link provides access to the SWeeP matrix derived from the EPCs of the complete proteomes available on NCBI.

1. Obtaining bacterial proteomes

A total of 17,724 complete bacterial proteomes (translated CDSs) were obtained from the NCBI RefSeq on October 17, 2023 (available at <ftp.ncbi.nlm.nih.gov/refseq/release>)

2. Identifying an Evolutionarily-relevant Protein Core (EPC)

Evolutionarily-relevant Protein Core (EPC) were obtained by applying the methodology described in Aryel Oliveira’s thesis (Oliveira, 2022) Pipeline and supporting files used to generate the EPC: Zenodo link [DOI: 10.5281/zenodo.14850984].

3. Vectorizing protein sequences with rSWeeP

To facilitate testing, we provide a small example dataset in this repository GitHub link. Download the file ‘rSWeePdata_sample17k.zip’, unzip it and provide the path='pathto/rSWeePdata_sample17k/inputs/' to the folder.

library(rSWeeP)
# supply here the path to MLSA data
folder = 'pathto/rSWeePdata_sample17k/inputs/'

As a first step, vectorize the sequences using SWeeP. The study used the following masks and output length:

  mask1 = c(1,1,0,1,1,0,1,0,0,0)
  mask2 = c(0,0,1,0,0,1,0,1,1,1)
  psz = 7500

But for a lighter test use the following mask and parameters:

library(rSWeeP)
mask1 = c(1,1,0,1)
mask2 = c(0,1,1,1)
seqtype = 'AA' # amino acid type
psz =  750 # output length
ncores = 2 # number of availabel cores

swA = SWeePlite(folder, seqtype=seqtype, bin=F, mask=mask1, psz=psz, norm='logNeg', lowRAMmode=T, ncores=ncores, verbose = FALSE)

## Starting projection. Please wait.

swB = SWeePlite(folder, seqtype=seqtype, bin=F, mask=mask2, psz=psz, norm='logNeg', lowRAMmode=T, ncores=ncores, verbose = FALSE)

## Starting projection. Please wait.

Concatenate the matrices and save the vectors

sw = swA
sw$info$timeElapsed =  swA$info$timeElapsed + swB$info$timeElapsed
sw$info$mask =  list()
sw$info$mask[[1]] = swA$info$mask
sw$info$mask[[2]] = swB$info$mask

sw$info$headers = c(swA$info$headers,swB$info$headers)
sw$proj = rbind(swA$proj,swB$proj)

# saving
saveRDS(sw,'SWeePVectorsConcatenated.rds')

4. Performing dimensionality reduction with Principal Component Analysis (PCA);

Here, you can use the SWeeP vectors concatenated generated above or use the complete vectors available at Zenodo.

In our study, we found that the first PCA component is highly correlated with sequence length and is therefore not evolutionarily informative; as a result, it was discarded. The phylogenies also yielded better results when PC-1 was removed. If you wish to study sequences of similar length, this step may not be necessary.

PCA = prcomp(sw$proj)

# Discarding PC-1
PCA = PCA$x[,-1]

5. Calculating pairwise Pearson correlation distances and other distances

Load the metadata.

The metadata includes information from various taxonomic levels and it is necessary to calculate the quality metrics for the phylogenies. Provide the metadata address.

metadatafile = "pathto/rSWeePdata_sample17k/Metadata_sample17k.csv"

Calculate the distance matrices and generate the correspondent phylogenetic trees

In the study, phylogenetic reconstruction was performed with ccphylo, which is considerably more efficient for very large datasets. Here we present an R implementation for generating trees using the Neighbor-joining method, in order to facilitate reproducibility. For large-scale analyses, we recommend using ccphylo.

The Euclidean and Pearson matrices

library(ape)

tr = list()

# EUCLIDEAN
mdist1   = dist(PCA,method='euclidean')
tr[[1]] = nj(mdist1)

# PEARSON
mdist2   = cor(t(PCA))
tr[[2]] = nj(mdist2)

The \(L_k\) distance (\(k=1/3\))

For this, use the following functions

# L_k distance of 1/3
#This step might take a while...
## < Distance matrices functions >
lk_norm <- function(x,y,k) {
  sum(abs(x-y)^k)^(1/k)
}

mdistLk = function(PCA){
  N = dim(PCA)[1]
  mdist = matrix(nrow=N,ncol=N)
  for (i in 1:N) {
    for (j in 1:N) {
      mdist[i,j]   = lk_norm(PCA[i,],PCA[j,],(1/3)) # k = 1/3
      mdist[j,i]   = mdist[i,j]
    }
  }
  return(mdist)
}

maxmin <- function(x) {
  (x - min(x)) / (max(x) - min(x))
}

# Generate the distance matrix
mdist3   = maxmin(mdistLk(PCA))
tr[[3]] = nj(mdist3)

Jaccard and Manhattan matrices

library(vegan)

# Jaccard
mdist4    = vegan::vegdist(PCA+min(PCA),method='jaccard')
tr[[4]] = nj(mdist4)

# Manhattan
mdist5 = dist(PCA, method = "manhattan")
tr[[5]] = nj(mdist5)

Checking the metrics for the different phylogenetic trees

Metrics can be implemented to check the quality of the trees. The metrics are available:

PCCI: PhyloTaxonomic Consistency Cophenetic Index. Phylogenetic tree evaluation function, estimate of how grouped the samples of the same taxon are in the phylogenetic tree.
PMPG: Percentage of Mono or Paraphyletic Groups. Phylogenetic tree evaluation function, returns the percentage of Mono/Paraphyletic

To make it easier, use the feature below to generate scores automatically.

## < Calculate the PMPG and PCCI score for the phylogenetic tree >
TaxonBreakMetrics <- function(tr){
  library(rSWeeP)

  # c("superkingdom", "phylum", "class", "order", "family", "genus", "species")
  res = data.frame(taxon=c('2phylum','3class','4order','5family','6genus'))
  is_tip <- tr$edge[,2] <= length(tr$tip.label)
  ordered_tips <- tr$edge[is_tip, 2]

  x=tr$tip.label[ordered_tips]

  y=as.data.frame(strsplit(x,'|',fixed=T))
  names(y)=NULL

  apeM=NULL
  apeP=NULL
  PMPG_ = NULL
  PCCI_ = NULL
  for (k in 2:6) {
    print(paste('starting k=',k))
    tr2 = tr
    tr2$tip.label = unlist(y[(k),])
    PCCI_ = c(PCCI_,PCCI(tr2)$mean)
    aux = PMPG(tr2)
    PMPG_ = c(PMPG_,aux$metric)
    apeM = c(apeM,aux$percMono)
    apeP = c(apeP,aux$percPara)
  }
  res = cbind(res,apeM)
  res = cbind(res,apeP)
  res = cbind(res,PMPG_)
  res = cbind(res,PCCI_)
  x = c('mean',colMeans(res[,2:5]))
  res = rbind(res,x)

  return(res) 
}

Use this function to obtain scores at each desired taxonomic level.
Enter “i” according to the desired distance matrix (1:“Euclidean”, 2:“Pearson”, 3:“\(L_k\)”, 4:“Jaccard”, 5:“Manhattan”), and “level” according to the desired taxonomic level (‘phylum’,‘class’,‘order’,‘family’ or ‘genus’).

# Selecting the taxonomic level to apply the PCCI and PMPG metrics
level = 'family'
i = 1 # where 1 = "Euclidean"

# Extract th GCF codes from labels (file names)
GCF_order = list()
for (i in 1:5){
    GCF_order[[i]] = tr[[i]]$tip.label
}

# Obtain the PCCI and PMPG metrics
submetadata = metadata[,level]
tr[[i]]$tip.label = submetadata[match(GCF_order[[i]], metadata$GCF)]
valPCCI = PCCI(tr[[i]]) 
valPMPG = PMPG(tr[[i]])

REVER METADATA EEE FAZER ÁVORE COM UM NIVEL TAXONOMICO EMBAIXO - PODE SER ORDEM

See the values

valPCCI

## $tab
##                     taxa      cost
## 1           Vibrionaceae 0.9375000
## 2          Moraxellaceae 0.9375000
## 3           Yersiniaceae 0.7500000
## 4         Legionellaceae 1.0000000
## 5            Frankiaceae 0.7500000
## 6         Rickettsiaceae 1.0000000
## 7         Clostridiaceae 1.0000000
## 8       Acetobacteraceae 1.0000000
## 9       Mycobacteriaceae 0.7500000
## 10    Corynebacteriaceae 0.7500000
## 11        Morganellaceae 0.7500000
## 12      Xanthomonadaceae 1.0000000
## 13 Peptostreptococcaceae 0.7500000
## 14     Salinibacteraceae 0.7500000
## 15     Streptomycetaceae 0.9375000
## 16      Streptococcaceae 0.9375000
## 17  Propionibacteriaceae 0.7500000
## 18     Helicobacteraceae 1.0000000
## 19          Rhizobiaceae 0.7500000
## 20     Planctomycetaceae 1.0000000
## 21            Thermaceae 0.7500000
## 22           Bacillaceae 0.8333333
## 23      Tsukamurellaceae 0.7500000
## 24      Caulobacteraceae 0.7500000
## 25    Bifidobacteriaceae 0.7500000
## 26      Actinomycetaceae 0.7500000
## 27   Syntrophomonadaceae 1.0000000
## 28          Trueperaceae 0.7500000
## 29         Chlamydiaceae 0.8867188
## 30      Lactobacillaceae 1.0000000
## 31        Leptospiraceae 0.7500000
## 32  Metamycoplasmataceae 1.0000000
## 33           Waddliaceae 0.7500000
## 34       Segniliparaceae 0.7500000
## 35       Brachyspiraceae 0.7500000
## 36      Burkholderiaceae 1.0000000
## 37     Xanthobacteraceae 0.7500000
## 38    Enterobacteriaceae 0.7500000
## 39          Listeriaceae 0.7500000
## 40    Blattabacteriaceae 0.7500000
## 
## $mean
## [1] 0.8430013

valPMPG

## $tab
##                     taxa  mono  para
## 1           Vibrionaceae FALSE FALSE
## 2          Moraxellaceae FALSE FALSE
## 3           Yersiniaceae FALSE FALSE
## 4         Legionellaceae  TRUE FALSE
## 5            Frankiaceae FALSE FALSE
## 6         Rickettsiaceae  TRUE FALSE
## 7         Clostridiaceae  TRUE FALSE
## 8       Acetobacteraceae  TRUE FALSE
## 9       Mycobacteriaceae FALSE FALSE
## 10    Corynebacteriaceae FALSE FALSE
## 11        Morganellaceae FALSE FALSE
## 12      Xanthomonadaceae  TRUE FALSE
## 13 Peptostreptococcaceae FALSE FALSE
## 14     Salinibacteraceae FALSE FALSE
## 15     Streptomycetaceae FALSE FALSE
## 16      Streptococcaceae FALSE  TRUE
## 17  Propionibacteriaceae FALSE FALSE
## 18     Helicobacteraceae  TRUE FALSE
## 19          Rhizobiaceae FALSE FALSE
## 20     Planctomycetaceae  TRUE FALSE
## 21            Thermaceae FALSE FALSE
## 22           Bacillaceae FALSE FALSE
## 23      Tsukamurellaceae FALSE FALSE
## 24      Caulobacteraceae FALSE FALSE
## 25    Bifidobacteriaceae FALSE FALSE
## 26      Actinomycetaceae FALSE FALSE
## 27   Syntrophomonadaceae  TRUE FALSE
## 28          Trueperaceae FALSE FALSE
## 29         Chlamydiaceae  TRUE FALSE
## 30      Lactobacillaceae  TRUE FALSE
## 31        Leptospiraceae FALSE FALSE
## 32  Metamycoplasmataceae  TRUE FALSE
## 33           Waddliaceae FALSE FALSE
## 34       Segniliparaceae FALSE FALSE
## 35       Brachyspiraceae FALSE FALSE
## 36      Burkholderiaceae  TRUE FALSE
## 37     Xanthobacteraceae FALSE FALSE
## 38    Enterobacteriaceae FALSE FALSE
## 39          Listeriaceae FALSE FALSE
## 40    Blattabacteriaceae FALSE FALSE
## 
## $percMono
## [1] 30
## 
## $percPara
## [1] 2.5
## 
## $metric
## [1] 32.5

6. Visualizing Phylogenetic Trees

You can see the phylogenies, based on euclidean distance with family taxonomic level for example, using

library(ggtree)

# Selecting the taxonomic level to apply the PCCI and PMPG metrics
level = 'family'
i = 1 # where 1 = "Euclidean"

newtr = tr[[i]]
submetadata = metadata[,level]
newtr$tip.label = submetadata[match(GCF_order[[i]], metadata$GCF)]

library(ggtree)
p<- ggtree(newtr,branch.length="none")+ geom_tiplab()+ xlim(0,50)
p

References

De Pierri, C.R., Voyceik, R., Santos de Mattos, L.G.C. et al. SWeeP: representing large biological sequences datasets in compact vectors. Sci Rep 10, 91 (2020). DOI: 10.1038/s41598-019-55627-4

Oliveira, A.M.R. Prospecção in silico de bactérias diazotróficas noduladoras a partir de genomas completos. PhD Thesis. Federal University of Parana. 2022.

Bacterial Proteomes Analysis Protocol with rSWeeP

Camila Pereira Perico

2026-09-05