# GLAMOUR: Graph Learning over Macromolecule Representations

Somesh Mohapatra, Joyce An & Rafael Gómez-Bombarelli
Department of Materials Science and Engineering
Massachusetts Institute of Technology
Cambridge, MA 02139, USA
<EMAIL_ADDRESS>Corresponding author<EMAIL_ADDRESS>
###### Abstract

The near-infinite chemical diversity of natural and artificial macromolecules
arises from the vast range of possible component monomers, linkages, and
polymers topologies. This enormous macromolecular diversity contributes to the
ubiquity and indispensability of macromolecules but hinders the development of
general machine learning methods with macromolecules as input. To address
this, we developed GLAMOUR, a framework for chemistry-informed graph
representation of macromolecules that enables quantifying structural
similarity, and interpretable supervised learning for macromolecules.

##  1 Introduction

Machine learning (ML) applications to individual macromolecule classes, such
as proteins, have been very successful but typically rely on sequence-based
representations that effectively capture structures with linear architectures
(Schissel et al., 2021; Alley et al., 2019). However, these methods do not
extend to macromolecules with complex topologies, such as glycans and
biohybrid sequence-defined polymers, which exhibit non-linear structure and
higher levels of monomer and linkage diversity. Graphs are a natural and more
general macromolecule representation, with nodes representing monomers and
edges representing connecting bonds. ML over graph representations has
achieved state-of-the-art results across several fields (Hamilton et al.,
2017) (see Section A for a survey of the field) and in chemistry and life
sciences, graph neural networks (GNNs) have become the modern workhorse for
molecular property prediction (Yang et al., 2019; Jin et al., 2020).
Attribution methods are compatible with GNNs to highlight how input features
are relevant to model predictions of target properties, such as molecular
olfactory descriptors (Sanchez-Lengeling et al., 2020; 2019).

In this work, we developed a graph representation for macromolecules. The
representation generalizes ideas of macromolecule similarity from sequence
alignment to structural similarity between complex topologies. Using chemical
similarity between monomers through cheminformatic fingerprints and exact
graph edit distances (GED) or graph kernels to compare topologies, it allows
for quantification of the chemical and structural similarity of two arbitrary
macromolecule topologies. Further, the representation was coupled to
supervised GNN models to learn structure-property relationships in glycans and
anti-microbial peptides. The graph representation handles linear, branched and
cyclic topologies along with any monomer composition. Attribution methods
highlight the regions of the macromolecules and the substructures within the
monomers that are most responsible for the predicted properties.

##  2 Results and Discussion

We developed a generalized text file system to convert a macromolecule
structure into a machine-readable format (Figure 1A, Section C). The text file
has three sections – SMILES, MONOMERS, and BONDS, inspired by the PDB file
format (Berman, 2000). Under SMILES, monomer and bond names followed by the
stereochemical SMILES are noted. MONOMERS enumerates indices of all nodes
numbered from 1 to n, where n is the total number of monomers, followed by the
monomer names. Similarly, BONDS lists indices of connections between monomer
indices, followed by bond names. In this way, we are able to incorporate
complexity from the level of individual atoms to the full macromolecular
structure.

Figure 1: Macromolecules represented as graphs enable similarity computation
and interpretable machine learning. A. Macromolecular structures are converted
into text files. The text files are parsed into network graphs, where monomers
are nodes and bonds are edges. Each text file enumerates SMILES for monomers
and bonds, node indices corresponding to monomers, and pairs of node indices
for bonds. Node and edge attributes are fingerprints of respective molecules.
B. Pair-wise similarity matrix obtained for the library of macromolecules, and
dimensionality reduction followed by clustering of similarity vectors
corresponding to each macromolecule. C. Graph neural networks learn over
macromolecule graphs for various tasks. D. Graph attribution provides insight
into the model’s decision-making process by highlighting the relative
importance of different nodes/monomers (denoted as node size) in making a
specific prediction.

In our experiments, the macromolecules were then processed as NetworkX graphs
(Hagberg et al., 2008), with monomers as nodes and bonds as edges (Section D).
In line with earlier work where fingerprint-based monomer representations
worked effectively for macromolecule property prediction, the monomer and bond
molecules were featurized using standard Extended Connectivity FingerPrints
(Schissel et al., 2021). The fingerprints capture the atomic connectivity of
the monomer/bond molecule as a series of bits using circular atom
neighborhoods, for each constituent node or edge of the macromolecule graph.
This representation enables the encoding of the raw chain connectivity of
macromolecules with explicit featurization of the stereochemistry and
topology. Moreover, it provides a single framework to represent both natural
and synthetic, linear and non-linear macromolecules.

To compute the similarity between 2 or more macromolecule graphs, we used GED
and graph kernel (Figure 1B, Section E). GED (Abu-Aisheh et al., 2015)
computes the similarity between two graphs by assigning node and edge
substitution scores, similar to local sequence alignment methods, such as
BLAST (Altschul et al., 1990). Instead of evolutionary statistics-based
substitution matrices like BLOSUM62, we used Tanimoto chemical similarity
matrices that compute the similarity between molecular fingerprints. Tanimoto
similarity also extends to unnatural monomers. Since computing exact graph
edit distances is an NP-hard problem, we used propagation attribute graph
kernels to obtain approximate similarity matrices for large datasets (Neumann
et al., 2016; Siglidis et al., 2020). This graph kernel captures local monomer
node information and propagates this information along the bond edges, making
it an ideal choice for macromolecule graphs. We have demonstrated the
similarity computation for a linear glycan with six additional glycans of
different topology and/or monomer chemistry, as well as with itself, using
Tanimoto chemical similarity matrix (Figure 2A, B).

Because quantifying graph similarity enables unsupervised learning over
similarities rather than raw features, we used dimensionality reduction
methods, such as uniform manifold approximation and projection (UMAP) (McInnes
et al., 2018) and t-stochastic neighbor embeddings (t-SNE) (Van Der Maaten &
Hinton, 2008), over the similarity matrices for unsupervised learning (Section
F). For glycans with immunogenicity labels, we observed that the non-
immunogenic and immunogenic glycans are in nearly distinct regions (Figure
2C).

Figure 2: Analysis using GLAMOUR provides deeper insights into immunogenicity
of glycans. A. Pairwise chemical similarity between glycan monomers using
Tanimoto distance on ECFP fingerprints is typically low. B. Similarities
between a linear glycan consisting of six glucose monomers with itself and six
additional glycans with modified topology and monomer chemistry were computed
using exact graph edit distances, combined with the Tanimoto substitution
matrices (A). C. 2-component UMAP over the similarity matrix, colored by
immunogenicity, for 1313 glycans with immunogenicity labels. D. Mean ROC-AUC
curve for classification of immunogenic and non-immunogenic glycans, averaged
across all 25 Attentive FP models (5 random seeds, with top 5 sets of
hyperparameters), with the standard deviation shaded in light blue too
insignificant to be visible in the graph. E. Attribution analysis using
Attentive FP model architecture and integrated gradients attribution method
for a representative glycan shows that xylose (Xyl) contributes the most to
immunogenicity. In xylose, the substructure centered on the anomeric carbon,
C1, is shown to be the key contributor. A similar substructure analysis has
been shown for fucose. In the glycan graph, the size of the node corresponds
to the importance, and the color corresponds to the monomer. In the
fingerprint, only the indices with positive contribution to immunogenicity
have been visualized, with the color corresponding to the normalized
importance

For supervised learning, we trained five different GNN model architectures to
classify glycans by immunogenicity and taxonomy levels, and predict the anti-
microbial activity of peptides (Section G). Each model architecture was
trained over fingerprint and one-hot node and edge attributes on 60%,
validated on 20%, and tested on held-out 20% data set, all determined via
random splits. Hyperparameter optimization was carried out for 1000 iterations
using SigOpt (Clark & Hayes, 2019). To report final metrics on validation and
test datasets, for each model architecture, we averaged values obtained from
25 individual models – top 5 hyperparameter sets, each trained with 5 random
seeds. For glycans, our models outperformed metrics for classification for
four out of eight tasks and achieved comparable results for the rest four,
against results reported in the literature (Figure 2D, Table A.5) (Burkholz et
al., 2021).

Attribution analysis on the pre-trained models finds features – both nodes and
substructures – that are key to the model’s decision-making (Section H). To
avoid spurious hypotheses (McCloskey et al., 2019) that may occur in NN
attribution on molecular models, we chose the combination of GNN architecture
and attribution method, in this case Attentive FP - integrated gradients, that
produced the most consistent results across all immunogenic glycans over 25
replicates (Figure A.45). N-glycolylneuraminic acid, xylose, and fucose were
found to be the key monomers responsible for immunogenicity, consistent with
experimental findings due to the low prevalence of these monomers in human
glycans (Planinc et al., 2016) (Table A.7). For a representative glycan, we
observed that xylose, followed by galectin, were the monomers that contributed
most significantly to immunogenicity (Figure 2E). In addition to the
importance of individual nodes, attribution identifies the critical
substructures in the monomers that contributed most to immunogenicity, such as
the substructure centered on the anomeric carbon of xylose. These critical
features help understand the fundamental structure-function relationships that
underpin immunogenicity, postulate very explicit hypotheses that can be
validated in the lab, and may help in further design of immunogenic or non-
immunogenic scaffolds.

##  3 Conclusion

GLAMOUR provides a generalized method for representing macromolecules as
hierarchical graphs with molecular fingerprints to capture chemical
information which can be used to compute structural similarity between
macromolecules with different composition and topology, and perform
unsupervised and supervised learning. The unsupervised learning enables
visualization of the complex landscape of different classes of macromolecules
and understanding of the subtle differences between similar macromolecules.
The attribution analysis helps in cracking open the black-box of supervised
GNN models, which in turn can help elucidate fundamental design principles of
otherwise opaque structure-property relationships and assist with hypothesis
generation for future experimental studies. The codebase and data sets have
been open-sourced (Mohapatra, 2021). We therefore expect that this toolkit
will be used by both experimentalists and computational practitioners in
chemistry, biology and materials science for a variety of macromolecule
property prediction tasks, ultimately driving the design of improved
macromolecules.

## Data availability

Data used for training of the model is available at
https://github.com/learningmatter-mit/GLAMOUR, and archived in the Zenodo
repository (Mohapatra, 2021).

## Code availability

All the code used for model training and analysis is available at
https://github.com/learningmatter-mit/GLAMOUR, and archived in Zenodo
repository (Mohapatra, 2021).

## Acknowledgements

SM acknowledges funding from Novo Nordisk and MIT’s Jameel Clinic, JA
acknowledges funding from MIT’s UROP office. The authors thank Dr Daria Kim
and Sarah Antilla from the Department of Chemistry, MIT for valuable comments.

## References

  * Abu-Aisheh et al. (2015) Zeina Abu-Aisheh, Romain Raveaux, Jean-Yves Ramel, and Patrick Martineau.  An Exact Graph Edit Distance Algorithm for Solving Pattern Recognition Problems.  _4th International Conference on Pattern Recognition Applications and Methods_ , 2015.
  * Alley et al. (2019) Ethan C Alley, Grigory Khimulya, Surojit Biswas, Mohammed AlQuraishi, and George M Church.  Unified rational protein engineering with sequence-based deep representation learning.  _Nature methods_ , 16(12):1315–1322, 2019.
  * Altschul et al. (1990) Stephen F Altschul, Warren Gish, Webb Miller, Eugene W Myers, and David J Lipman.  Basic local alignment search tool.  _Journal of molecular biology_ , 215(3):403–410, 1990.
  * Altschul et al. (1997) Stephen F Altschul, Thomas L Madden, Alejandro A Schäffer, Jinghui Zhang, Zheng Zhang, Webb Miller, and David J Lipman.  Gapped blast and psi-blast: a new generation of protein database search programs.  _Nucleic acids research_ , 25(17):3389–3402, 1997\.
  * Aoki et al. (2003) Kiyoko F Aoki, Atsuko Yamaguchi, Yasushi Okuno, Tatsuya Akutsu, Nobuhisa Ueda, Minoru Kanehisa, and Hiroshi Mamitsuka.  Efficient tree-matching methods for accurate carbohydrate database queries.  _Genome informatics_ , 14:134–143, 2003.
  * Berman (2000) H M Berman.  The Protein Data Bank.  _Nucleic Acids Research_ , 28(1):235–242, jan 2000.  ISSN 13624962.  doi: 10.1093/nar/28.1.235.  URL http://www.ncbi.nlm.nih.gov/pubmed/10592235http://www.pubmedcentral.nih.gov/articlerender.fcgi?artid=PMC102472https://academic.oup.com/nar/article-lookup/doi/10.1093/nar/28.1.235.
  * Bileschi et al. (2019) Maxwell L. Bileschi, David Belanger, Drew Bryant, Theo Sanderson, Brandon Carter, D. Sculley, Mark A. DePristo, and Lucy J. Colwell.  Using Deep Learning to Annotate the Protein Universe.  _bioRxiv:10.1101/626507_ , 2019.
  * Blumenthal & Gamper (2020) David B Blumenthal and Johann Gamper.  On the exact computation of the graph edit distance.  _Pattern Recognition Letters_ , 134:46–57, 2020.
  * Bojar et al. (2021) Daniel Bojar, Rani K. Powers, Diogo M. Camacho, and James J. Collins.  Deep-Learning Resources for Studying Glycan-Mediated Host-Microbe Interactions.  _Cell Host and Microbe_ , 29(1):132–144.e3, 2021\.
  * Boratyn et al. (2019) Grzegorz M Boratyn, Jean Thierry-Mieg, Danielle Thierry-Mieg, Ben Busby, and Thomas L Madden.  Magic-blast, an accurate rna-seq aligner for long and short reads.  _BMC bioinformatics_ , 20(1):1–19, 2019.
  * Borgwardt et al. (2020) Karsten Borgwardt, Elisabetta Ghisu, Felipe Llinares-lópez, Leslie O Bray, and Bastian Rieck.  Graph Kernels.  _arXiv:2011.03854v2_ , 2020.
  * Burkholz et al. (2021) Rebekka Burkholz, John Quackenbush, and Daniel Bojar.  Using graph convolutional neural networks to learn a representation for glycans.  _Cell Reports_ , 35(11):109251, 2021.
  * Chen et al. (2021) Lihua Chen, Ghanshyam Pilania, Rohit Batra, Tran Doan Huan, Chiho Kim, Christopher Kuenneth, and Rampi Ramprasad.  Polymer informatics: Current status and critical next steps.  _Materials Science and Engineering: R: Reports_ , 144:100595, 2021.
  * Clark & Hayes (2019) Scott Clark and Patrick Hayes.  SigOpt Web page.  https://sigopt.com, 2019.  URL https://sigopt.com.
  * Coff et al. (2020) Lachlan Coff, Jeffrey Chan, Paul A Ramsland, and Andrew J Guy.  Identifying glycan motifs using a novel subtree mining approach.  _BMC bioinformatics_ , 21(1):42, 2020.
  * Drefahl (2011) Axel Drefahl.  Curlysmiles: a chemical language to customize and annotate encodings of molecular and nanodevice structures.  _Journal of cheminformatics_ , 3(1):1–7, 2011\.
  * Eddy (2004) Sean R Eddy.  Where did the BLOSUM62 alignment score matrix come from?  _Nature biotechnology_ , 22(8):1035, 2004.
  * Gilmer et al. (2017) Justin Gilmer, Samuel S. Schoenholz, Patrick F. Riley, Oriol Vinyals, and George E. Dahl.  Neural Message Passing for Quantum Chemistry.  _arXiv:1704.01212v2_ , 2017.
  * Hagberg et al. (2008) Aric A Hagberg, Daniel A Schult, and Pieter J Swart.  Exploring Network Structure, Dynamics, and Function using NetworkX.  In Gaël Varoquaux, Travis Vaught, and Jarrod Millman (eds.), _Proceedings of the 7th Python in Science Conference_ , pp. 11–15, Pasadena, CA USA, 2008.
  * Hamilton et al. (2017) William L Hamilton, Rex Ying, and Jure Leskovec.  Representation learning on graphs: Methods and applications.  _arXiv preprint arXiv:1709.05584_ , 2017.
  * Heller et al. (2015) Stephen R Heller, Alan McNaught, Igor Pletnev, Stephen Stein, and Dmitrii Tchekhovskoi.  Inchi, the iupac international chemical identifier.  _Journal of cheminformatics_ , 7(1):1–34, 2015\.
  * Hosoda et al. (2017) Masae Hosoda, Yukie Akune, and Kiyoko F. Aoki-Kinoshita.  Development and application of an algorithm to compute weighted multiple glycan alignments.  _Bioinformatics_ , 33(9):1317–1323, 2017.
  * Jaakkola et al. (2000) Tommi Jaakkola, Mark Diekhans, and David Haussler.  A discriminative framework for detecting remote protein homologies.  _Journal of Computational Biology_ , 7(1-2):95–114, 2000.
  * Jin et al. (2020) Wengong Jin, Regina Barzilay, and Tommi Jaakkola.  Hierarchical generation of molecular graphs using structural motifs.  In _International Conference on Machine Learning_ , pp. 4839–4848. PMLR, 2020.
  * Kearnes et al. (2016) Steven Kearnes, Kevin McCloskey, Marc Berndl, Vijay Pande, and Patrick F. Riley.  Molecular graph convolutions: moving beyond fingerprints.  _Journal of Computer-Aided Molecular Design_ , 30(8):595–608, 2016.
  * Kim et al. (2018) Chiho Kim, Anand Chandrasekaran, Tran Doan Huan, Deya Das, and Rampi Ramprasad.  Polymer Genome: A Data-Powered Polymer Informatics Platform for Property Predictions.  _Journal of Physical Chemistry C_ , 122(31):17575–17585, 2018.
  * Kim et al. (2019) Sunghwan Kim, Jie Chen, Tiejun Cheng, Asta Gindulyte, Jia He, Siqian He, Qingliang Li, Benjamin A Shoemaker, Paul A Thiessen, Bo Yu, et al.  Pubchem 2019 update: improved access to chemical data.  _Nucleic acids research_ , 47(D1):D1102–D1109, 2019.
  * Kipf & Welling (2019) Thomas N. Kipf and Max Welling.  Semi-supervised classification with graph convolutional networks.  _5th International Conference on Learning Representations, ICLR 2017 - Conference Track Proceedings_ , 2019.
  * Kriege et al. (2020) NM Kriege, FD Johansson, and C Morris.  A survey on graph kernels. appl netw sci 5 (1): 6, 2020.
  * Landrum (2006) Greg Landrum.  RDKit: Open-source cheminformatics, 2006.
  * Lee et al. (2021) Chee-Kong Lee, Chengqiang Lu, Yue Yu, Qiming Sun, Chang-Yu Hsieh, Shengyu Zhang, Qi Liu, and Liang Shi.  Transfer learning with graph neural networks for optoelectronic properties of conjugated oligomers.  _The Journal of Chemical Physics_ , 154(2):024906, 2021.
  * Leslie et al. (2001) Christina Leslie, Eleazar Eskin, and William Stafford Noble.  The spectrum kernel: A string kernel for svm protein classification.  In _Biocomputing 2002_ , pp. 564–575. World Scientific, 2001.
  * Li et al. (2010) Limin Li, Wai Ki Ching, Takako Yamaguchi, and Kiyoko F. Aoki-Kinoshita.  A weighted q-gram method for glycan structure classification.  _BMC Bioinformatics_ , 11:1–6, 2010.
  * Lin et al. (2019) Tzyy-Shyang Lin, Connor W Coley, Hidenobu Mochigase, Haley K Beech, Wencong Wang, Zi Wang, Eliot Woods, Stephen L Craig, Jeremiah A Johnson, and Julia A Kalow.  BigSMILES: a structurally-based line notation for describing macromolecules.  _ACS central science_ , 5(9):1523–1531, 2019\.
  * McCloskey et al. (2019) Kevin McCloskey, Ankur Taly, Federico Monti, Michael P. Brenner, and Lucy J. Colwell.  Using attribution to decode binding mechanism in neural network models for chemistry.  _Proceedings of the National Academy of Sciences of the United States of America_ , 116(24):11624–11629, 2019.  ISSN 10916490.  doi: 10.1073/pnas.1820657116.
  * McInnes et al. (2017) Leland McInnes, John Healy, and Steve Astels.  hdbscan: Hierarchical density based clustering.  _Journal of Open Source Software_ , 2(11):205, 2017.
  * McInnes et al. (2018) Leland McInnes, John Healy, Nathaniel Saul, and Lukas Großberger.  UMAP: Uniform Manifold Approximation and Projection.  _Journal of Open Source Software_ , 3(29):861, 2018.
  * Mohapatra (2021) Somesh Mohapatra.  learningmatter-mit/glamour: v0.1, August 2021.  URL https://doi.org/10.5281/zenodo.5237237.
  * Needleman & Wunsch (1970) Saul B. Needleman and Christian D. Wunsch.  A general method applicable to the search for similarities in the amino acid sequence of two proteins.  _Journal of Molecular Biology_ , 48(3):443–453, 1970.
  * Neumann et al. (2016) Marion Neumann, Roman Garnett, Christian Bauckhage, and Kristian Kersting.  Propagation kernels: efficient graph kernels from propagated information.  _Machine Learning_ , 102(2):209–245, 2016.
  * Pedregosa et al. (2011) Fabian Pedregosa, Gaël Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, Olivier Grisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, et al.  Scikit-learn: Machine learning in python.  _the Journal of machine Learning research_ , 12:2825–2830, 2011.
  * Pirtskhalava et al. (2016) Malak Pirtskhalava, Andrei Gabrielian, Phillip Cruz, Hannah L Griggs, R Burke Squires, Darrell E Hurt, Maia Grigolava, Mindia Chubinidze, George Gogoladze, Boris Vishnepolsky, et al.  Dbaasp v. 2: an enhanced database of structure and antimicrobial/cytotoxic activity of natural and synthetic peptides.  _Nucleic acids research_ , 44(D1):D1104–D1112, 2016.
  * Planinc et al. (2016) Ana Planinc, Jonathan Bones, Bieke Dejaegher, Pierre Van Antwerpen, and Cédric Delporte.  Glycan characterization of biopharmaceuticals: updates and perspectives.  _Analytica chimica acta_ , 921:13–27, 2016.
  * Qiao et al. (2020) Bo Qiao, Somesh Mohapatra, Jeffrey Lopez, Graham M. Leverick, Ryoichi Tatara, Yoshiki Shibuya, Yivan Jiang, Arthur France-Lanord, Jeffrey C. Grossman, Rafael Gómez-Bombarelli, Jeremiah A. Johnson, and Yang Shao-Horn.  Quantitative Mapping of Molecular Substituents to Macroscopic Properties Enables Predictive Design of Oligoethylene Glycol-Based Lithium Electrolytes.  _ACS Central Science_ , 6(7):1115–1128, 2020\.
  * Riesen & Bunke (2009) Kaspar Riesen and Horst Bunke.  Approximate graph edit distance computation by means of bipartite graph matching.  _Image and Vision Computing_ , 27(7):950–959, 2009.
  * Rogers & Hahn (2010) David Rogers and Mathew Hahn.  Extended-Connectivity Fingerprints.  _J. Chem. Inf. Model._ , 50(5):742–754, 2010\.
  * Sanchez-Lengeling et al. (2019) Benjamin Sanchez-Lengeling, Jennifer N Wei, Brian K Lee, Richard C Gerkin, Alán Aspuru-Guzik, and Alexander B Wiltschko.  Machine learning for scent: learning generalizable perceptual representations of small molecules.  _arXiv preprint arXiv:1910.10685_ , 2019.
  * Sanchez-Lengeling et al. (2020) Benjamin Sanchez-Lengeling, Jennifer Wei, Brian Lee, Emily Reif, Peter Y Wang, Wei Qian, Kevin Mccloskey, Lucy Colwell, and Alexander Wiltschko.  Evaluating Attribution for Graph Neural Networks.  _Proceeedings of the Neural Information Processing Systems Conference_ , (NeurIPS), 2020.
  * Schissel et al. (2021) Carly K Schissel, Somesh Mohapatra, Justin M Wolfe, Colin M Fadzen, Kamela Bellovoda, Chia-Ling Wu, Jenna A Wood, Annika B Malmberg, Andrei Loas, Rafael Gómez-Bombarelli, et al.  Deep learning to design nuclear-targeting abiotic miniproteins.  _Nature Chemistry_ , pp. 1–9, 2021.
  * Schütt et al. (2017) Kristof T. Schütt, P.-J J. Pieter-Jan P.-J Pieter-Jan Kindermans, Huziel E. Sauceda, Stefan Chmiela, Alexandre Tkatchenko, and Klaus-Robert R. K.-R Müller.  SchNet: A continuous-filter convolutional neural network for modeling quantum interactions.  _Advances in Neural Information Processing Systems_ , 2017-Decem(1):992–1002, jun 2017.  URL http://arxiv.org/abs/1706.08566www.quantum-machine.org.
  * Seo et al. (2018) Seokjun Seo, Minsik Oh, Youngjune Park, and Sun Kim.  Deepfam: deep learning based alignment-free method for protein family modeling and prediction.  _Bioinformatics_ , 34(13):i254–i262, 2018.
  * Shrikumar et al. (2017) Avanti Shrikumar, Peyton Greenside, and Anshul Kundaje.  Learning important features through propagating activation differences.  _arXiv:1704.02685_ , 2017.
  * Siglidis et al. (2020) Giannis Siglidis, Giannis Nikolentzos, Stratis Limnios, Christos Giatsidis, and Michalis Vazirgiannis.  GraKeL : A Graph Kernel Library in Python.  _Journal of Machine Learning Research_ , 21:1–5, 2020.
  * Smith & Waterman (1981) Temple F Smith and Michael S Waterman.  Identification of common molecular subsequences.  _Journal of molecular biology_ , 147(1):195–197, 1981.
  * St John et al. (2019) Peter C. St John, Caleb Phillips, Travis W. Kemper, A. Nolan Wilson, Yanfei Guan, Michael F. Crowley, Mark R. Nimlos, and Ross E. Larsen.  Message-passing neural networks for high-throughput polymer screening.  _Journal of Chemical Physics_ , 150(23), 2019.
  * Sugiyama et al. (2018) Mahito Sugiyama, M. Elisabetta Ghisu, Felipe Llinares-López, and Karsten Borgwardt.  Graphkernels: R and Python packages for graph comparison.  _Bioinformatics_ , 34(3):530–532, 2018.
  * Sundararajan et al. (2017) Mukund Sundararajan, Ankur Taly, and Qiqi Yan.  Axiomatic attribution for deep networks.  _arXiv:1703.01365_ , 2017.
  * Van Der Maaten & Hinton (2008) Laurens Van Der Maaten and Geoffrey Hinton.  Visualizing Data using t-SNE.  _Journal of Machine Learning Research_ , 9:2579–2605, 2008\.
  * Velickovic et al. (2017) Petar Velickovic, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Lio, and Yoshua Bengio.  Graph Attention Networks.  _arXiv:1710.10903_ , 2017.
  * Wang et al. (2019) Minjie Wang, Da Zheng, Zihao Ye, Quan Gan, Mufei Li, Xiang Song, Jinjing Zhou, Chao Ma, Lingfan Yu, Yu Gai, Tianjun Xiao, Tong He, George Karypis, Jinyang Li, and Zheng Zhang.  Deep Graph Library: A Graph-Centric, Highly-Performant Package for Graph Neural Networks.  _arXiv:1909.01315_ , 2019.
  * Xiong et al. (2020) Zhaoping Xiong, Dingyan Wang, Xiaohong Liu, Feisheng Zhong, Xiaozhe Wan, Xutong Li, Zhaojun Li, Xiaomin Luo, Kaixian Chen, Hualiang Jiang, and Mingyue Zheng.  Pushing the boundaries of molecular representation for drug discovery with the graph attention mechanism.  _Journal of Medicinal Chemistry_ , 63(16):8749–8760, 2020.
  * Yamanishi et al. (2007) Yoshihiro Yamanishi, Francis Bach, and Jean-Philippe Vert.  Glycan classification with tree kernels.  _Bioinformatics_ , 23(10):1211–1216, 2007.
  * Yang et al. (2019) Kevin Yang, Kyle Swanson, Wengong Jin, Connor Coley, Philipp Eiden, Hua Gao, Angel Guzman-Perez, Timothy Hopper, Brian Kelley, Miriam Mathea, et al.  Analyzing learned molecular representations for property prediction.  _Journal of chemical information and modeling_ , 59(8):3370–3388, 2019.
  * Zeng et al. (2018) Minggang Zeng, Jatin Nitin Kumar, Zeng Zeng, Ramasamy Savitha, Vijay Ramaseshan Chandrasekhar, Kedar Hippalgaonkar, and Electronic Materials.  Graph Convolutional Neural Networks for Polymers.  _arXiv:1811.06231v1_ , 2018.
  * Zhang et al. (2012) Tianhong Zhang, Hongli Li, Hualin Xi, Robert V Stanton, and Sergio H Rotstein.  Helm: a hierarchical notation language for complex biomolecule structure representation, 2012.
  * Zhang et al. (2010) Zhenjie Zhang, Marios Hadjieleftheriou, Beng Chin Ooi, and Divesh Srivastava.  Bed-tree: an all-purpose index structure for string similarity search based on edit distance.  In _Proceedings of the 2010 ACM SIGMOD International Conference on Management of data_ , pp. 915–926, 2010.
  * Zhou et al. (2018) Jie Zhou, Ganqu Cui, Zhengyan Zhang, Cheng Yang, Zhiyuan Liu, Lifeng Wang, Changcheng Li, and Maosong Sun.  Graph Neural Networks: A Review of Methods and Applications.  _arXiv:1812.08434v4_ , 2018.

APPENDIX

##  Appendix A Related work

_Representation._ Macromolecules can be represented as strings, using
hierarchical editing language for macromolecules (HELM) (Zhang et al., 2012),
International Union of Pure and Applied Chemistry (IUPAC) international
chemical identifier (InChI) (Heller et al., 2015), CurlySMILES (Drefahl, 2011)
(where SMILES is Simplified Molecular-Input Line-Entry System) and BigSMILES
(Lin et al., 2019). Line notations do not always support all topologies,
require a fair amount of customization, and have non-canonical variants. One
exception to this is linear biological macromolecules, such as proteins and
DNA/RNA, which are represented as sequences of one/three-letter monomer codes.
In a recent attempt, glycans (non-linear biological macromolecules) were
represented as sequences, where groups of monosaccharides were clubbed into
’glycowords’ and placed in hierarchical brackets (Bojar et al., 2021).
Representing glycans presents a significant advance in notation but falls
short in terms of comparison (limited to linear glycans) and machine learning
tasks.

Hierarchical fingerprinting is another approach, which follows a hierarchy of
atomic (categorical encoding of presence/absence of contiguous atom-sets),
physicochemical (e.g., fraction of rings, molecular surface area), and
morphological descriptors (e.g., length of side chain, length of main chain)
(Kim et al., 2018). This approach is limited in its coverage of chemical
space, ability to differentiate stereochemistry, and capturing long range
through-space interactions.

_Similarity computation._ In recent times, there have been significant
advances in computation of similarity using graph edit distances (Blumenthal &
Gamper, 2020) (GED) and graph kernels (Borgwardt et al., 2020; Kriege et al.,
2020). Development of software packages, such as graphkernels (Sugiyama et
al., 2018) and GraKeL (Siglidis et al., 2020), has provided fast
implementations of graph kernels.

For linear biological macromolecules, such as proteins, DNA/RNA and linear
glycans, there are several works for computation of sequence similarities
(Bojar et al., 2021; Altschul et al., 1990; 1997; Boratyn et al., 2019).
Usually, sequence alignment is done using Smith-Waterman (Smith & Waterman,
1981) or Needleman-Wunsch (Needleman & Wunsch, 1970) algorithms, and scored
with substitution matrices, such as BLOSUM62 (Eddy, 2004) (for proteins) and
GLYSUM (Bojar et al., 2021) (for glycans). These substitution matrices are
based on evolutionary statistics thereby biasing the scoring towards the
statistical frequency of a particular monomer’s occurrence in the course of
evolution, rather than chemical similarity. Apart from sequence alignment in
linear macromolecules, edit distances (Riesen & Bunke, 2009; Zhang et al.,
2010), linear kernels (Leslie et al., 2001; Jaakkola et al., 2000) and deep
learning methods (Bileschi et al., 2019; Seo et al., 2018) have been proposed
to compute similarity. In the case of non-linear macromolecules, alignment of
glycans has been explored using q-grams (Li et al., 2010), tree matching
methods (Aoki et al., 2003; Hosoda et al., 2017; Coff et al., 2020), and using
tree kernels (Yamanishi et al., 2007). Unfortunately, the aforementioned
methods are limited to biological macromolecules, and do not extend to the
general macromolecular chemical space. Moreover, existing tools do not allow
incorporation of unnatural monomers and non-linear topologies (except
glycans).

_Machine learning._ The field of graph neural networks (GNN) has seen
substantial developments in both model architecture and attribution. The
different model architectures have demonstrated state-of-the-art results
across various domains - graph convolutional network (GCN), graph attention
network (GAT), message passing neural network (MPNN), SchNet, Weave (a variant
of GCN), and AttentiveFP (a variant of GAT) (Zhou et al., 2018; Schütt et al.,
2017; Kearnes et al., 2016; Xiong et al., 2020). Graph attribution has been
recently studied quantitatively across four metrics – accuracy, stability,
faithfulness and consistency - for a wide variety of tasks and model
architectures (Sanchez-Lengeling et al., 2020).

For macromolecular property prediction, Polymer Genome and similar works using
hierarchical fingerprints predict glass transition temperature, dielectric
point and other macromolecular properties (Kim et al., 2018; Chen et al.,
2021). There have been several attempts to extrapolate macromolecular property
by training over monomer input features (St John et al., 2019; Qiao et al.,
2020; Lee et al., 2021). GCN over macromolecule graphs with one-hot node and
edge attributes has been shown to outperform fingerprint-based models (Zeng et
al., 2018). For fingerprint-based models, the lack of a good macromolecule
representation limits the model performance, while the GCN model is limited by
the lack of chemical information in the graph representation.

##  Appendix B Dataset processing

###  B.1 Glycans

####  B.1.1 Dataset download

A dataset of 19299 glycans was accessed and downloaded from GlycoBase
(accessed on November 2, 2020) (Bojar et al., 2021). The file contained
GlycoBase ID, sequence, link (N, O, free, or none), species, and
immunogenicity information for each glycan.

For each glycan sequence string the brackets denote branches, with the point
of attachment/bonding of the branch as the monomer immediately after the
brackets. The first element within the bracket is the monomer most distant
from the point of attachment, and the last element within the bracket is the
abbreviation of the bond that connects the branch to the original main chain.
Nested brackets indicate additional sub-branches off of branches, and multiple
sets of brackets next to each other indicate several branches off of the same
monomer.

####  B.1.2 Dataset pre-processing

7 modifications and 152 deletions of glycan sequences were made before the
process for situations such as an unequal number of opening and closing
brackets and dangling branches without specified connectivity. Additional
glycan sequences were removed due to missing SMILES sequences for a number of
monomers. The original GlycoBase.csv file was curated to reflect the
modification and deletion changes, resulting in a total of 19147 glycans.

Using the curated database, we visualized the distribution of species of
origin, link, and immunogenicity of the glycans (Figures A.1, A.2, A.3).

Figure A.1: Top 50 most common glycan species of origin, sorted in descending
order of count with y-axis on a logarithmic scale. Figure A.2: Four types of
glycan links sorted in descending order of count with y-axis on a logarithmic
scale. Figure A.3: Three different types of glycan immunogenicity labels
sorted in descending order of count, with y-axis on a logarithmic scale.

####  B.1.3 SMILES compilation

The chemical composition and formula of the 959 monomers (also known as
glycoletters in GlycoBase) and 53 bonds were expressed as isomeric simplified
molecular-input line-entry system (SMILES) sequences.

In the dataset retrieved from GlycoBase, all position-specific information
about monomer modifications was removed. The supplementary information of
SweetTalk contains all the raw glycan sequences before position-specific
modification information was removed (Bojar et al., 2021). While insufficient
information is provided to directly match the raw sequences with the
corresponding edited sequences, the raw sequences can be used to look for
trends and most common modification positions for each monomer. A dictionary
was created to describe the number of times each monomer appears in all the
raw sequences. A couple of terms should be defined for consistency:

  * •

“Monosaccharides” refer to the individual glycoletters without any
modifications, such as D-glucose (D-Glu) and galactose (D-Gal). The term
“monosaccharide” is simplistic because it also covers alcohols, acids, and
other classes of molecules, but the term will be employed for the sake of
consistency and clarity.

  * •

“Modifications” refer to substitutions or additions to the original
monosaccharide sequence such as a nitrogen-linked acetyl group (NAc) and
oxygen-linked methyl group (OMe). In this dataset, monosaccharides can have
between 0-4 modifications.

  * •

“Monomers” refer to the combination of the monosaccharide and modification(s),
also known as glycoletters in GlycoBase.

Using the raw monomers dictionary, the positions of the modifications for each
monomer in the corresponding SMILES were set using a list of consistent rules
outlined below. The rules are listed in order of priority.

  1. 1.

For each monomer/glycoletter, if at least one monomer with the same
monosaccharide and set of modifications exists in the raw monomers dictionary,
assume the set of positions with the highest frequency. If no monomer exists
in the raw monomers dictionary with the same monosaccharide and set of
modifications, proceed to the following steps.

  2. 2.

If the monomer has a hydroxyl group at the 2 position, is not a furanose or
ketose, and contains a nitrogen-linked modification (N, NBut, NMe, etc.),
assign the first occurrence of a nitrogen-linked modification to the 2
position.

  3. 3.

If the monomer is a hexose or deoxy-hexose and contains a pyruvate acetal
(OPyr), assume that the acetal connects the 4 and 6 positions.

  4. 4.

If the monomer contains an O-linked phosphate (OP) or sulfate (OS)
modification, search the raw monomers dictionary for instances of the
monosaccharide with only the OP/OS modification and assume the most frequent
position.

  5. 5.

If the monomer contains a variant of the O-linked phosphate modification
(OPEtn, OPPEtn, etc.), assign the OP-variant modification the same position as
the most common position for the OP modification on the monosaccharide.

  6. 6.

If the monomer contains two modifications, search the raw glycans dictionary
for instances of the monosaccharide with only each modification separately but
not at the same time. Assign each modification the position of highest
frequency. If either monosaccharide and modification combination does not
exist in the raw monomers dictionary and the monomer is a hexose, assign
positions in the following order: 2, 4, 3, 6. If the position is already
occupied from previous steps or does not exist in a deoxy hexose, skip to the
position with next highest priority.

  7. 7.

If the monomer is a hexose or deoxy-hexose and contains more than two
modifications, assign positions in the following order: 2, 3, 4, 6. If the
position is already occupied from previous steps or does not exist in a deoxy
hexose, skip to the position with next highest priority.

  8. 8.

Assume all amino acids are connected to the monosaccharide via the oxygen on
the carboxyl. If the connectivity is not specified for a group following the
amino acid (CysAc, AlaFo), assume that the group following the amino acid is
connected to the amino acid via the amine group.

  9. 9.

For neuraminic acid (Neu), ketodeoxynononic acid (Kdn), pseudaminic acid
(Pse), legionaminic acid (Leg), and other similar ketose-based acids, assign
modification positions in the following order: 1, 4.

  10. 10.

For fructofuranose and similar furanoses for which the 1-position is not part
of the ring, assign modifications in the following order: 1, 3, 4.

  11. 11.

For alcohols, follow the same rules that apply to the oxidized form of the
alcohol. For example, for glucitol follow the same rules that apply to
glucose.

  12. 12.

Assume all rare monosaccharides (denoted as Sug) are hexoses with no specified
stereochemistry.

Because the 959 monomers contained some repeat SMILES sequences with different
monomer names, 13 redundant names were deleted so that each distinct monomer
SMILES only appears once in the SMILES compilation.

The SMILES for the 53 bond types differ in stereochemistry alone (alpha, beta,
or unspecified) but all have the same chemical composition of a glycosidic
bond. Each bond is expressed as a variation of the SMILES sequence CC(OC)CC,
which consists of the glycosidic bond C-O-C with one of the C atoms also
connected to both a methyl and ethyl group. The chiral C with the four
different attached groups is used to specify the alpha or beta stereochemistry
displayed in the bond name. The stereochemistry at the tetrahedral C is
consistently S for all alpha bonds, R for all beta bonds, and not specified
for all unspecified bonds. The 53 bond names were condensed into 3 distinct
bond types differing in stereochemistry alone.

###  B.2 Anti-microbial peptides

####  B.2.1 Dataset download

A dataset of 15864 antimicrobial peptides (AMPs), including 15450 monomers,
200 multimers, and 214 multi-peptides, from the Database of Antimicrobial
Activity and Structure of Peptides (DBAASP) was accessed and downloaded on
October 6, 2020 (Pirtskhalava et al., 2016).

Each peptide was represented in the dataset as an individual .json file
containing information about the peptide ID, name, sequence(s), unusual amino
acids, connectivity, terminal modifications, complexity, synthesis type,
target groups and objects, and target species. The term “target species” is
used loosely and actually encompasses both species and sub-classification
information such as subspecies, strain for bacteria, and forma speciales for
fungi. For each target species for any given peptide, the dataset provides the
antimicrobial peptide concentration in units of mostly either $\mu$M or
$\mu$g/mL as a function of four unique variables: activity measure, salt type,
medium, and CFU.

The dataset includes three different types of peptide complexities: monomers,
multimers, and multi-peptides. In the DBAASP dataset, monomers consist of a
single sequence, multimers between 2 and 4 separate sequences connected via
interchain bonds, and multi-peptides between 2 and 4 separate sequences
connected not via covalent bonds but instead weaker intermolecular forces.

####  B.2.2 Dataset pre-processing

The information from each .json file was inputted into a combined table and
the antimicrobial peptide concentrations for each target species was converted
into a float value by removing symbols like ‘$>$’ and ‘$<$’, taking the
average whenever the dataset provides a range, and disregarding uncertainty
values.

9 peptide monomer types included in a total of 86 peptides were removed due to
ambiguity of the molecular structure, bringing the total number of peptides to
15,778. This condensed dataset was further processed to visualize the
distribution of target species data points. For each “target species,” the
species name was separated from any sub-classifications (subspecies, strain,
forma speciales, serovar, pathovar, biovar, etc.).

Using the processed database, we visualized the distribution of species,
activity measure, salt types, mediums, units, and CFU in the target species
data points (Figures A.4-A.9).

Figure A.4: Top 50 species with greatest number of data points sorted in
descending order of count, with y-axis on a logarithmic scale. Figure A.5: Top
50 activity measures with greatest number of data points sorted in descending
order of count, with y-axis on a logarithmic scale. Figure A.6: Top 50 salt
types with greatest number of data points sorted in descending order of count,
with y-axis on a logarithmic scale. The blank label for the first column
indicates the absence of any added salt. Figure A.7: Top 50 mediums with
greatest number of data points sorted in descending order of count, with
y-axis on a logarithmic scale. Figure A.8: All peptide concentration units
sorted in descending order of count, with y-axis on a logarithmic scale.
Figure A.9: Top 50 CFUs with greatest number of data points sorted in
descending order of count, with y-axis on a logarithmic scale.

As evident from the bar plots, the vast majority of data points belong to
either _Escherichia coli_ or _Staphylococcus aereus_. To understand whether
the overall dataset could be further filtered to limit the variability in
conditions for more meaningful comparison of peptide concentrations, we
visualized the distribution of activity measure and salt type for both _E.
coli_ (Figures A.10, A.11) and _S. aereus_ (Figures A.12, A.13).

Figure A.10: Top 50 activity measures for _E. coli_ data points with greatest
number of occurrences sorted in descending order of count, with y-axis on a
logarithmic scale. Figure A.11: Top 50 salt types for _E. coli_ data points
with greatest number of occurrences sorted in descending order of count, with
y-axis on a logarithmic scale. The blank label for the first column indicates
the absence of any added salt. Figure A.12: Top 50 activity measures for _S.
aereus_ data points with greatest number of occurrences sorted in descending
order of count, with y-axis on a logarithmic scale. Figure A.13: Top 50 salt
types for _S. aereus_ data points with greatest number of occurrences sorted
in descending order of count, with y-axis on a logarithmic scale. The blank
label for the first column indicates the absence of any added salt.

For both _E. coli_ and _S. aereus_ , the majority of data points have an
activity measure of MIC and no added salt. Due to the clear majority for both
activity measure and salt type as well as the difficulty of interconversion
between the units of $\mu$M and $\mu$g/mL, the overall dataset was filtered
into two separate, smaller datasets: one containing only data points falling
into all of the categories of _E. coli_ , MIC, no salt, and $\mu$M data points
and the other containing only data points falling into all of the categories
of _S. aereus_ , MIC, no salt, and $\mu$M. Each of the two sub-datasets was
further filtered to remove any peptides for which the standard deviation of
the peptide concentrations meeting the criteria for the sub-dataset was not
zero, resulting in a total of 4,445 peptides for the _E. coli_ dataset and
3,686 peptides for the _S. aereus_ dataset.

####  B.2.3 SMILES compilation

The chemical composition and formula of the 332 amino acids, 93 N-termini
modifications, 37 C-termini modifications, and 22 bonds were expressed as
isomeric simplified molecular-input line-entry system (SMILES) sequences. The
majority of isomeric SMILES sequences were found on PubChem from the National
Library of Medicine (Kim et al., 2019). The remaining SMILES were collected
through review of literature cited in the DBAASP peptide cards of peptides
containing the corresponding monomers or bonds.

For the monomers, no additional sequences beyond those of the 20 standard
L-amino acids and unusual amino acids, C-termini modifications, and N-termini
modifications found in the Help section of the DBAASP website were included in
the SMILES compilation. 8 monomers had repeat names between the three
different types of monomers, and the names were changed so that each monomer
corresponds with a unique name. For the bonds, the sidechain-mainchain bonds
(SMB) were split into 7 subtypes determined through review of the literature.
In addition, 6 additional bond types were added for cases in which the
connection between the N or C-termini modification and terminal amino acid is
not an amide bond. The 7 SMB subtypes, 6 additional terminal bond types, and 9
intrachain bond types in the Help section of the DBAASP website resulted in a
total of 22 bonds.

##  Appendix C Text file system

###  C.1 Format

The text files to convert each macromolecule structure into machine-readable
format consist of three sections: the SMILES sequences for each unique monomer
or bond, the positions of each monomer, and the two monomer positions
connected by each bond. Each section starts with a header to indicate the
start of a new part. The first section contains the abbreviation of each
unique monomer or bond in the glycan followed by the corresponding SMILES
sequence, with each entry on a separate line. The monomers section consists of
the monomer position followed by the monomer abbreviation. The bonds section
contains the two connectivity positions in the glycan for the bond separated
by a space, followed by the bond abbreviation (Figure A.14).

Figure A.14: Standard text file format with three sections: SMILES of unique
monomers and bonds, monomer positions and abbreviations, and bond
connectivities and abbreviations.

###  C.2 Text File Parser

A text file parser converts the macromolecule information stored in the .txt
file to a NetworkX graph with monomers expressed as nodes and bonds expressed
as edges. The parser goes through the .txt file line by line, stores the
monomer information in a dictionary with keys as integer positions and values
as monomer abbreviations, and stores the bond information in a dictionary with
keys as tuples containing bond connectivities and values as bond
abbreviations. Afterwards, the reader uses NetworkX to add each key in the
monomer dictionary as a node and each key in the bond dictionary as an edge,
storing the abbreviations as attributes for the corresponding node or edge.
The resulting NetworkX graphs include both linear and highly branched
architectures (Figure A.15, A.16).

Figure A.15: NetworkX graph representations of glycans that consist of both
linear and highly branched architectures of varying molecular weights and
complexities. Figure A.16: NetworkX graph representations of peptides of all
three complexities of varying molecular weights.

##  Appendix D Node and edge attributes

###  D.1 Fingerprints – generation and optimization of hyperparameters

We used RDKit to generate stereochemical extended connectivity fingerprints
(Landrum, 2006; Rogers & Hahn, 2010). Radius and number of bits were optimized
by calculating mean and standard deviation, and visualizing the distribution
of Tanimoto similarity of all monomers in the glycans dataset (Figure
A.17A-C). We aimed to obtain fingerprints with lower number of bits and
optimal radius that could represent the monomer aptly in both similarity
computation and graph neural network models. For 64 bits, we observed that the
mean similarity was as high as 0.4 and standard deviation for similarity went
up after radius 3 indicating higher hash collision and lesser
differentiability. For 128 bits, the mean and standard deviation plateaued
around 0.3 and 0.08, respectively. Further, the similarity distribution was
qualitatively spread over a larger range, as compared to fingerprints with
larger number of bits. For fingerprint bits longer than 128, we observed lower
mean and standard deviation, and decreasing spread in the similarity
distribution. The decreasing trend indicates that the fingerprints with bits
higher than 128 are equally dissimilar, thus, if used, will lead to glycans
with equally high dissimilar scores. Hence, we chose fingerprints with 128
bits and radius 3 to generate node attributes for glycans. For edge
attributes, we had 3 types of glycosidic bonds differing by the
stereochemistry alone, hence, we chose fingerprints with 16 bits and radius 3.
We employed similar reasoning for the peptides to select 128 bits and radius 3
for both monomers and bonds.

We observed that there were 3 sets with 7 glycans which had exact fingerprints
(Figure A.17D). Since, the difference was in the number of carbon atoms in the
aliphatic chain, we used the fingerprints as is.

Figure A.17: A. Mean, B. standard deviation and C. distribution of Tanimoto
similarity of all monomers in the glycans dataset, calculated using
stereochemical extended connectivity fingerprints of different radii and bits.
D. Glycans in the same box have the exact fingerprint for radius 3 and 128
bits.

###  D.2 One-hot encoding benchmark

One-hot encodings of the 946 monomer and 3 bond types for glycans and the 461
monomer and 22 bond types for peptides were also employed as feature types for
benchmarking with featurization using molecular fingerprints. The dimensions
of the node and edge one-hot encoding features are 946 and 3, respectively,
for glycans and 461 and 22, respectively, for peptides.

##  Appendix E Similarity computation

###  E.1 Analysis of graphs

Most glycan graphs are sparse (Figure A.8). Complete graphs have density of 1,
while graphs with density $\ >$ 0.5 have been defined as dense (Borgwardt et
al., 2020).

Figure A.18: Distribution of number of A. nodes, B. edges and C. density for
glycan graphs.

###  E.2 Similarity matrix

We computed the$\ (n\times n)$ similarity matrix for all glycans with labels
on at least one taxonomic level using propagation attribute kernel in GraKeL
(Figure A.18) (Siglidis et al., 2020). Each pair of graph similarity was
computed for a maximum of 100 iterations. This resulted in 5% of the pairs
being assigned a 0 similarity (10% of all indices in the similarity matrix are
0).

It may be noted that the computation of similarity using graph kernels is way
more accessible than graph edit distances. The current computation was done on
in minutes (wall time), parallelized across 24 cores. From visual inspection
using htop, only about 30% of 2 cores were being used at any particular time,
and less than 5% of all other cores were used, although there were 24 jobs
running in parallel.

Figure A.19: 2D plot for the$\ (n\times n)$ similarity matrix of glycans. This
is not symmetric because each row is normalized by its maximum.

###  E.3 Dissimilarity computation using exact graph edit distances

We performed a grid search over combinations of possible node/edge
substitution costs and multipliers for node/edge insertion/deletion costs to
find an optimal set of values. Instead of similarity, we use the opposite –
dissimilarity – as it is more intuitive when compared to a baseline case. We
use the representative glycans in Figure 2B for this experiment, and present
them here in Figure A.20 with the names we use in Table A.1.

For node/edge substitution, the Tanimoto distance between the stereochemical
fingerprints, a value in the range of 0 and 1 (where 0 is a perfect match), is
multiplied by the substitution cost to obtain the final graph edit distance.
For node/edge insertion/deletion, a constant value which is the node/edge
insertion/deletion is added to the graph edit distance. The baseline case
where the edit distance is calculated for the macromolecule with itself is 0.

Figure A.20: Representative glycans used in the hyperparameter search for
dissimilarity computation using exact graph edit distances, and similarity
computation using graph kernels. Table A.1: Dissimilarities computed using
exact graph edit distances using different node/edge insertion/deletion and
substitution multipliers. The remarks column notes down the hyperparameters,
and baseline and other classes of glycans. Abbreviations - Hparam:
hyperparameter, L: linear, B: branched, glc: glucose, xyl: xylose, fuc:
fucose.

Remarks | Hyperparameters/Glycans | Cases
---|---|---
HParam1 | Insertion/Deletion Cost | 1 | 1 | 1 | 1 | 3 | 3 | 3 | 3 | 5 | 5 | 5 | 5 | 10 | 10 | 10 | 10
HParam2 | Substitution Cost | 1 | 3 | 5 | 10 | 1 | 3 | 5 | 10 | 1 | 3 | 5 | 10 | 1 | 3 | 5 | 10
Baseline | L-glc | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0
Topology | M-glc | 1 | 1 | 1 | 1 | 3 | 3 | 3 | 3 | 5 | 5 | 5 | 5 | 10 | 10 | 10 | 10
B1-glc | 2 | 2 | 2 | 2 | 6 | 6 | 6 | 6 | 10 | 10 | 10 | 10 | 20 | 20 | 20 | 20
B2-glc | 4 | 4 | 4 | 4 | 12 | 12 | 12 | 12 | 20 | 20 | 20 | 20 | 40 | 40 | 40 | 40
Monomer Chemistry | L-glc-xyl | 0.68 | 6 | 6 | 6 | 0.68 | 2.05 | 3.42 | 18 | 0.68 | 2.05 | 3.42 | 6.84 | 0.68 | 2.05 | 3.42 | 6.84
L-glc-xyl-fuc | 1.43 | 12 | 12 | 12 | 1.43 | 4.28 | 7.13 | 36 | 1.43 | 4.28 | 7.13 | 14.26 | 1.43 | 4.28 | 7.13 | 14.26
Topology and Chemistry | B-glc-xyl-fuc | 6.17 | 14 | 14 | 14 | 14.17 | 18.5 | 22.84 | 42 | 22.17 | 26.5 | 30.84 | 41.67 | 42.17 | 46.5 | 50.84 | 61.67

###  E.4 Similarity computation using graph kernels

We performed a grid search over hyperparameters of propagation attribute
kernel – bin width 1, 3, 10, 100, and the preserved distance metric on local
sensitive hashing ‘L1-norm’, ‘L2-norm’. We used 30 as the maximum number of
iterations for the kernel computation for this experiment. We use the
representative glycans Figure A.20, and present the results in Table A.2.

Table A.2: Similarities computed using graph kernels with different values of
bin width and distance metric. The remarks column notes down the
hyperparameters, and baseline and other classes of glycans. Abbreviations -
Hparam: hyperparameter, L: linear, B: branched, glc: glucose, xyl: xylose,
fuc: fucose.

Remarks | Hyperparameters/Glycans | Cases
---|---|---
HParam1 | Bin width | 1 | 1 | 3 | 3 | 10 | 10 | 100 | 100
HParam2 | Distance metric | L1 | L2 | L1 | L2 | L1 | L2 | L1 | L2
Baseline | L-glc | 230 | 230 | 230 | 230 | 230 | 230 | 230 | 230
Topology | M-glc | 61 | 61 | 61 | 61 | 61 | 61 | 61 | 61
B1-glc | 56 | 56 | 56 | 56 | 56 | 56 | 56 | 56
B2-glc | 51 | 51 | 51 | 51 | 51 | 51 | 51 | 51
Monomer Chemistry | L-glc-xyl | 50 | 50 | 50 | 50 | 50 | 55 | 56 | 61
L-glc-xyl-fuc | 39 | 39 | 39 | 39 | 45 | 44 | 51 | 61
Topology and Chemistry | B-glc-xyl-fuc | 23 | 23 | 23 | 23 | 35 | 33 | 46 | 51

##  Appendix F Dimensionality reduction

###  F.1 Hyperparameter optimization for Uniform Manifold Approximation and
Projection (UMAP)

Number of neighbors was optimized for 2-component UMAP dimensionality
reduction of similarity vectors (Figure A.21) (McInnes et al., 2018). From
visual inspection, UMAP with 128 neighbors seems to resolve into optimal size
and number of clusters. The subplot shows distinct regions for the immunogenic
and non-immunogenic glycans. We note that there is more to similarity than
graph kernel distances, and observe that it has been effectively captured
using the GNN models (Appendix G).

Figure A.21: Visualization of scatter plots for UMAP components, obtained by
dimensionality reduction of similarity vectors. The number of neighbors for
each UMAP computation has been noted in the title of the respective sub-plot.
Coloration is by immunogenicity (red: immunogenic, blue: non-immunogenic
glycans).

UMAP dimensionality reduction for higher number of components does not provide
more information. The 3 components UMAP looked similar to the 2 components,
with slightly more disentangled families (Figure A.22). We limited our visual
analysis to 3 components. To check if more components can help in finding
distinct clusters, we did dimensionality reduction for $\\{2,3,5,10,30,50\\}$
components and let HDBSCAN - an unsupervised clustering algorithm – to figure
out how many clusters are there (McInnes et al., 2017). We noted that the
number of clusters are pretty similar, and in low 400s (Figure A.23). The high
number of clusters indicates the diversity of the space, and the differences
in terms of taxonomy. As a further check to see if the distribution of glycans
in different clusters is different, we plotted the histograms of the glycans
assigned to each cluster. Across all the components, the histograms seem to be
consistent with the number of glycans in each of them (Figure A.24).

Figure A.22: Visualization of A. 2 and B. 3 components UMAP, colored by
domain. C. Glycans corresponding to individual domains are shown, with the
text noting the domain, number of glycans, and transparency ($\alpha$) of each
point on a scale of 0 to 1, where 1 is opaque Figure A.23: Number of HDBSCAN
unsupervised clusters obtained from UMAP dimensionality reduction for
different components. Figure A.24: Distribution of glycans in each cluster for
HDBSCAN unsupervised clusters obtained from UMAP dimensionality reduction for
different components. The components have been noted as titles for the sub-
plots.

###  F.2 T-SNE benchmark

We benchmarked the dimensionality reduction results obtained from UMAP against
a broad range of t-stochastic neighbor embeddings (t-SNE) models (Van Der
Maaten & Hinton, 2008). For the different models, we varied perplexity as
$\\{2,5,30,50,100\\}$, and number of steps as $\\{500,1000,5000\\}$. From the
scatter plot, colored by immunogenicity labels, we noted that dimensionality
reduction using t-SNE was not able to deduce the differences and getting the
glycans into distinct areas (Figure A.25).

Figure A.25: Visualization of scatter plots for t-SNE components, obtained by
dimensionality reduction of similarity vectors. The (perplexity, number of
steps) for each t-SNE computation has been noted in the title of the
respective sub-plot. Coloration is by immunogenicity (red: immunogenic, blue:
non-immunogenic glycans).

##  Appendix G Supervised learning with graph neural networks

###  G.1 Deep graphs library (DGL) graphs

Following featurization, NetworkX graphs were converted into undirected,
unweighted, and homogenous DGL graphs (Wang et al., 2019). For GCN and GAT
model architectures, self-loops were added to the DGL graphs to prevent silent
performance regression due to zero-in-degree nodes during training.

###  G.2 Model architectures

We performed graph classification using five distinct graph neural network
model architectures detailed below:

  * •

Weave (Kearnes et al., 2016)

  * •

Message Passing Neural Networks (MPNNs) (Gilmer et al., 2017)

  * •

Attentive FP (Xiong et al., 2020)

  * •

Graph convolutional networks (GCN) (Kipf & Welling, 2019)

  * •

Graph Attention Networks (GAT) (Velickovic et al., 2017)

While Weave, MPNN, and Attentive FP utilize both node and edge attributes in
prediction, GCN and GAT only consider node attributes. The models were trained
using implementations in the DGL LifeSci library (Wang et al., 2019). For
classification, the optimization was done by minimization of average cross-
entropy loss between batches and additional metrics such as F1 score, recall,
precision and accuracy were noted. For regression, the optimization was done
by minimization of root-mean-squared-error loss on the validation dataset.

###  G.3 Glycan graphs classification

####  G.3.1 Immunogenicity

Dataset. 1313 glycans in the database have immunogenicity labels, 631 of which
are immunogenic and 682 of which are not immunogenic. The training was
performed on 60%, validated on 20%, and tested on held-out 20% data.

Models. We classified immunogenicity using 5 model architectures combined with
2 different node and edge featurization types, for a total of 10 model
architecture-attribute pairs. For each benchmark, hyperparameter optimization
against minimization of binary cross entropy loss was performed on SigOpt, a
standardized hyperparameter optimization platform, for 1000 observations and
the 5 best sets of hyperparameters were extracted. Each model architecture and
featurization combination was trained using the 5 best sets of hyperparameters
from SigOpt using 5 distinct random seeds for splitting the dataset into
train-validation-test datasets, for a total of 25 trainings per model per
attribute type. The tables below report the metrics for the most optimal set
of hyperparameters, the values of the most optimal hyperparameters, and mean
metrics for each training across all 25 runs. All models achieve stellar
performance on all metrics, with little meaningful difference in performance
between fingerprint and one-hot encoding featurization (Figures A.26-A.28 ;
Table A.3).

Figure A.26: ROC-AUC curves for fingerprint and one-hot encoding-featurized
graphs. A. Mean ROC-AUC curve of all 25 fingerprint-featurized experiments (5
model architectures with 5 sets of hyperparameters for each architecture),
with the standard deviation shaded in light blue too insignificant to be
visible in the graph. B. Mean ROC-AUC curve of all 25 one-hot encoding-
featurized experiments. Figure A.27: Mean ROC-AUC curves for fingerprint-
featurized experiments for each of the five model architectures (Weave, MPNN,
Attentive FP, GCN, and GAT), with the standard deviation shaded in light blue
too insignificant to be visible. A standard deviation of 0.000 denotes a value
of$\ <$0.001. Figure A.28: Mean ROC-AUC curves for one-hot encoding-featurized
experiments for each of the five model architectures (Weave, MPNN, Attentive
FP, GCN, and GAT), with the standard deviation shaded in light blue too
insignificant to be visible. A standard deviation of 0.000 denotes a value
of$\ <$0.001. Table A.3: Test dataset metrics for the most optimal set of
hyperparameters, or the set of hyperparameters that results in the lowest
loss. For each metric, the mean $\mu$ and standard deviation $\sigma$ are
displayed across all 5 random seeds. “FP” denotes condensed fingerprint
featurization, and “One-hot” denotes one-hot encoding featurization.

Model | Attribute Type | ROC-AUC | F1 | Recall | Precision | Accuracy | Loss
---|---|---|---|---|---|---|---
|  | $\mu$ | $\sigma$ | $\mu$ | $\sigma$ | $\mu$ | $\sigma$ | $\mu$ | $\sigma$ | $\mu$ | $\sigma$ | $\mu$ | $\sigma$
Weave | FP | 0.990 | 0.005 | 0.956 | 0.011 | 0.955 | 0.012 | 0.957 | 0.011 | 0.956 | 0.011 | 0.133 | 0.106
One-hot | 0.992 | 0.005 | 0.962 | 0.007 | 0.962 | 0.007 | 0.962 | 0.007 | 0.962 | 0.007 | 0.105 | 0.055
MPNN | FP | 0.985 | 0.009 | 0.955 | 0.012 | 0.955 | 0.012 | 0.957 | 0.012 | 0.955 | 0.012 | 0.137 | 0.100
One-hot | 0.988 | 0.006 | 0.953 | 0.010 | 0.953 | 0.010 | 0.954 | 0.010 | 0.953 | 0.010 | 0.123 | 0.052
AttentiveFP | FP | 0.990 | 0.005 | 0.954 | 0.011 | 0.953 | 0.011 | 0.955 | 0.012 | 0.954 | 0.011 | 0.125 | 0.104
One-hot | 0.990 | 0.005 | 0.951 | 0.011 | 0.951 | 0.012 | 0.952 | 0.011 | 0.952 | 0.011 | 0.120 | 0.084
GCN | FP | 0.992 | 0.004 | 0.953 | 0.013 | 0.953 | 0.013 | 0.954 | 0.012 | 0.953 | 0.012 | 0.110 | 0.086
One-hot | 0.990 | 0.006 | 0.955 | 0.012 | 0.956 | 0.012 | 0.956 | 0.013 | 0.956 | 0.012 | 0.123 | 0.110
GAT | FP | 0.991 | 0.006 | 0.954 | 0.014 | 0.954 | 0.013 | 0.954 | 0.014 | 0.954 | 0.014 | 0.109 | 0.081
One-hot | 0.991 | 0.004 | 0.954 | 0.011 | 0.954 | 0.011 | 0.955 | 0.010 | 0.955 | 0.011 | 0.119 | 0.085

####  G.3.2 Taxonomy

Dataset. The taxonomy of the glycans was considered on eight levels: domain,
kingdom, phylum, class, order, genus, and species. First, the classification
of each species into the other seven taxonomic levels was obtained from the
supplementary tables of SweetOrigins (Bojar et al., 2021). For each glycan
with species information, taxonomic information was added in a new .csv file,
with multiple labels on the same taxonomic level separated by commas.
Afterwards, any labels with fewer than five glycans were removed, and any
species names ending with “sp” that are therefore genus labels in disguise
were filtered out to produce the final dataset for taxonomic training. The
final dataset consists of 8899 unique glycans with labels on at least one
taxonomic level encompassing 4 domains, 9 kingdoms, 33 phyla, 70 classes, 144
orders, 253 families, 400 genus, and 567 species. The training was performed
on 60%, validated on 20%, and tested on held-out 20% data.

Graph labels were generated as one-hot encodings for each taxonomic level,
with the length of each one-hot encoding tensor corresponding with the number
of unique labels in the taxonomic level. The taxonomy labels accommodate cases
in which glycans possess multiple labels within the same taxonomic level.

Models. Benchmarks of multi-label classification were performed using a
similar process as with immunogenicity classification for 5 different model
architectures, 2 different node and edge attribution types, and 8 different
taxonomic levels for a total of 80 model architecture–attribute–taxonomic
level combination. For each combination, hyperparameter optimization against
minimization of binary cross entropy loss was performed on SigOpt for 1000
observations and the 5 best sets of hyperparameters were extracted. Each model
architecture, featurization, and taxonomic level combination was trained using
the 5 best sets of hyperparameters from SigOpt using 5 distinct random seeds
for splitting the dataset into train-validation-test datasets, for a total of
25 trainings per combination. The pipeline and models achieve state-of-the-art
performance on multilabel taxonomic classification on all taxonomic levels
(Figures A.29-A.36, Table A.4). Additional supplementary tables display the
validation and test metrics for the most optimal set of hyperparameters, the
values of the most optimal hyperparameters, and mean test metrics for each
training across all 25 runs (Supplementary Tables 1-2).

Table A.4: The most optimal model architecture and node/edge attribute type
that results in the lowest loss for prediction of all taxonomic levels. For
each metric, the mean $\mu$ and standard deviation $\sigma$ are displayed
across all 5 random seeds. “FP” denotes condensed fingerprint featurization,
and “One-hot” denotes one-hot encoding featurization.

Taxonomic Level | Model, Attribute Type | ROC-AUC | F1 | Recall | Precision | Accuracy | Loss | Hamming Loss
---|---|---|---|---|---|---|---|---
$\mu$ | $\sigma$ | $\mu$ | $\sigma$ | $\mu$ | $\sigma$ | $\mu$ | $\sigma$ | $\mu$ | $\sigma$ | $\mu$ | $\sigma$ | $\mu$ | $\sigma$
Domain | Attentive FP, FP | 0.993 | 0 | 0.938 | 0.002 | 0.935 | 0.004 | 0.941 | 0.001 | 0.925 | 0.001 | 0.087 | 0.002 | 0.03 | 0.001
Kingdom | Attentive FP, FP | 0.994 | 0 | 0.912 | 0.001 | 0.891 | 0.003 | 0.934 | 0.002 | 0.892 | 0.005 | 0.056 | 0.003 | 0.019 | 0
Phylum | GCN, FP | 0.99 | 0 | 0.84 | 0.009 | 0.8 | 0.011 | 0.885 | 0.009 | 0.802 | 0.01 | 0.029 | 0.001 | 0.009 | 0
Class | GCN, FP | 0.986 | 0.001 | 0.741 | 0.008 | 0.666 | 0.01 | 0.836 | 0.006 | 0.67 | 0.009 | 0.022 | 0 | 0.007 | 0
Order | GCN, FP | 0.979 | 0.001 | 0.638 | 0.012 | 0.541 | 0.023 | 0.778 | 0.018 | 0.548 | 0.016 | 0.017 | 0 | 0.004 | 0
Family | GAT, FP | 0.975 | 0.003 | 0.557 | 0.015 | 0.456 | 0.017 | 0.715 | 0.015 | 0.469 | 0.019 | 0.013 | 0 | 0.003 | 0
Genus | GCN, One-hot | 0.963 | 0.003 | 0.513 | 0.015 | 0.397 | 0.015 | 0.723 | 0.020 | 0.412 | 0.010 | 0.009 | 0 | 0.002 | 0
Species | GCN, FP | 0.968 | 0.003 | 0.487 | 0.021 | 0.382 | 0.031 | 0.675 | 0.018 | 0.395 | 0.020 | 0.007 | 0 | 0.002 | 0

Figure A.29: ROC-AUC curves for classification on the domain level for each of
four model architectures (Weave, Attentive FP, GCN, and GAT), with the
standard deviation shaded in light blue too insignificant to be visible. A
standard deviation of 0.000 denotes a value of $\ <$0.001. A. Mean ROC-AUC
curves for fingerprint-featurized experiments, with each graph displaying the
mean and standard deviation across 5 sets of hyperparameters. B. Mean ROC-AUC
curves for one-hot encoding-featurized experiments. Figure A.30: ROC-AUC
curves for classification on the kingdom level for each of four model
architectures (Weave, Attentive FP, GCN, and GAT), with the standard deviation
shaded in light blue too insignificant to be visible. A standard deviation of
0.000 denotes a value of $\ <$0.001. A. Mean ROC-AUC curves for fingerprint-
featurized experiments, with each graph displaying the mean and standard
deviation across 5 sets of hyperparameters. B. Mean ROC-AUC curves for one-hot
encoding-featurized experiments. Figure A.31: ROC-AUC curves for
classification on the phylum level for each of four model architectures
(Weave, Attentive FP, GCN, and GAT), with the standard deviation shaded in
light blue too insignificant to be visible. A standard deviation of 0.000
denotes a value of $\ <$0.001. A. Mean ROC-AUC curves for fingerprint-
featurized experiments, with each graph displaying the mean and standard
deviation across 5 sets of hyperparameters. B. Mean ROC-AUC curves for one-hot
encoding-featurized experiments. Figure A.32: ROC-AUC curves for
classification on the class level for each of four model architectures (Weave,
Attentive FP, GCN, and GAT), with the standard deviation shaded in light blue
too insignificant to be visible. A standard deviation of 0.000 denotes a value
of $\ <$0.001. A. Mean ROC-AUC curves for fingerprint-featurized experiments,
with each graph displaying the mean and standard deviation across 5 sets of
hyperparameters. B. Mean ROC-AUC curves for one-hot encoding-featurized
experiments. Figure A.33: ROC-AUC curves for classification on the order level
for each of four model architectures (Weave, Attentive FP, GCN, and GAT), with
the standard deviation shaded in light blue too insignificant to be visible. A
standard deviation of 0.000 denotes a value of $\ <$0.001. A. Mean ROC-AUC
curves for fingerprint-featurized experiments, with each graph displaying the
mean and standard deviation across 5 sets of hyperparameters. B. Mean ROC-AUC
curves for one-hot encoding-featurized experiments. Figure A.34: ROC-AUC
curves for classification on the family level for each of four model
architectures (Weave, Attentive FP, GCN, and GAT), with the standard deviation
shaded in light blue too insignificant to be visible. A standard deviation of
0.000 denotes a value of $\ <$0.001. A. Mean ROC-AUC curves for fingerprint-
featurized experiments, with each graph displaying the mean and standard
deviation across 5 sets of hyperparameters. B. Mean ROC-AUC curves for one-hot
encoding-featurized experiments. Figure A.35: ROC-AUC curves for
classification on the genus level for each of four model architectures (Weave,
Attentive FP, GCN, and GAT), with the standard deviation shaded in light blue
too insignificant to be visible. A standard deviation of 0.000 denotes a value
of $\ <$0.001. A. Mean ROC-AUC curves for fingerprint-featurized experiments,
with each graph displaying the mean and standard deviation across 5 sets of
hyperparameters. B. Mean ROC-AUC curves for one-hot encoding-featurized
experiments. Figure A.36: ROC-AUC curves for classification on the species
level for each of four model architectures (Weave, Attentive FP, GCN, and
GAT), with the standard deviation shaded in light blue too insignificant to be
visible. A standard deviation of 0.000 denotes a value of $\ <$0.001. A. Mean
ROC-AUC curves for fingerprint-featurized experiments, with each graph
displaying the mean and standard deviation across 5 sets of hyperparameters.
B. Mean ROC-AUC curves for one-hot encoding-featurized experiments.

###  G.4 Benchmarking top GNN models against results reported in literature -
Glycans

For the top-performing models in Tables A.3 and A.4, hyperparameter
optimization on SigOpt was performed again using the same method as the
benchmarks in the work by Burkholz et al. (2021): training on 80% of the
dataset and reporting metrics on the remaining 20% used as a validation
dataset (Table A.5). The experiment on SigOpt was optimized through
minimization of the loss, and the remaining metrics (ROC-AUC, F1, recall,
precision, and accuracy) were obtained as stored metrics in the SigOpt
experiment. However, to obtain more robust metrics, we averaged across 25
models, instead of the 5 models in Burkholz et al. (2021). In 4 out of 8
tasks, the models presented in our work perform better, and are comparable for
the rest 4.

Table A.5: Validation metrics of the top-performing model architecture and
attribute combinations for immunogenicity and taxonomic levels obtained from
SigOpt using a 0.8, 0.2 train-validation split of the dataset. The subset
accuracy values are compared directly with the results for augmented models
presented in Table 1 of Burkholz et al. (2021). The higher value is bolded for
each task. In 4 out of 8 tasks, the models presented in our work perform
better, and are comparable for the rest 4.

Class | Paper | Model + Attribute Type | ROC-AUC | F1 | Recall | Precision | Accuracy | CE Loss | Hamming Loss
---|---|---|---|---|---|---|---|---|---
Immunogenicity | This Work | Weave, one-hot | 0.999 | 0.989 | 0.989 | 0.989 | 0.989 | 0.018 | -
Burkholz, et. al. | - | - | - | - | - | 0.942 | 0.149 | -
Domain | This Work | Attentive FP, fp | 0.994 | 0.946 | 0.941 | 0.951 | 0.935 | 0.081 | 0.027
Burkholz, et. al. | - | - | - | - | - | 0.939 | 0.181 | -
Kingdom | This Work | Attentive FP, fp | 0.997 | 0.936 | 0.895 | 0.922 | 0.921 | 0.043 | 0.014
Burkholz, et. al. | - | - | - | - | - | 0.907 | 0.301 | -
Phylum | This Work | GCN, fp | 0.991 | 0.845 | 0.806 | 0.889 | 0.808 | 0.029 | 0.009
Burkholz, et. al. | - | - | - | - | - | 0.841 | 0.603 | -
Class | This Work | GCN, FP | 0.987 | 0.776 | 0.713 | 0.850 | 0.718 | 0.021 | 0.006
Burkholz, et. al. | - | - | - | - | - | 0.745 | 0.976 | -
Order | This Work | GCN, fp | 0.979 | 0.638 | 0.541 | 0.778 | 0.548 | 0.017 | 0.004
Burkholz, et. al. | - | - | - | - | - | 0.582 | 1.728 | -
Family | This Work | GAT, fp | 0.984 | 0.617 | 0.525 | 0.748 | 0.542 | 0.011 | 0.003
Burkholz, et. al. | - | - | - | - | - | 0.535 | 2.051 | -
Genus | This Work | GCN, fp | 0.975 | 0.546 | 0.453 | 0.688 | 0.456 | 0.009 | 0.002
Burkholz, et. al. | - | - | - | - | - | 0.475 | 2.320 | -
Species | This Work | GCN, fp | 0.976 | 0.510 | 0.403 | 0.694 | 0.438 | 0.007 | 0.002
Burkholz, et. al. | - | - | - | - | - | 0.428 | 2.505 | -

###  G.5 Peptide graphs regression

Dataset. As explained in Appendix Section B.2, the overall peptides dataset
was filtered into two sub-datasets, one for E. coli and one for S. aereus. The
E. coli dataset contains 4,445 peptides, including 4330 monomers, 71
multimers, and 44 multi-peptides. The S. aereus dataset contains 3,686
peptides, including 3617 monomers, 61 multimers, and 8 multi-peptides. The
training for prediction of the target activity, or antimicrobial
concentration, of the peptides was performed on 60%, validated on 20%, and
tested on held-out 20% data.

Models. We performed target activity regression using 5 model architectures
combined with 2 different node and edge featurization types and 2 dataset
normalization types, for a total of 20 model combinations for each of the two
species. The 2 normalization methods were standard scaler and quantile
transform from scikit-learn (Pedregosa et al., 2011). For each benchmark,
hyperparameter optimization against minimization of root mean square error
(RMSE) was performed on SigOpt (Clark & Hayes, 2019), for 500 observations and
the 5 best sets of hyperparameters were extracted. Each model architecture and
featurization combination was trained using the 5 best sets of hyperparameters
from SigOpt using 5 distinct random seeds for splitting the dataset into
train-validation-test datasets, for a total of 25 trainings per model
combination. All models achieve stellar performance on all metrics, with
little meaningful difference in performance between fingerprint and one-hot
encoding featurization (Figures A.37-A.44, Appendix Table A.6). Additional
supplementary tables display the validation and test metrics for the most
optimal set of hyperparameters, the values of the most optimal
hyperparameters, and mean test metrics for each training across all 25 runs
(Supplementary Tables 3-4).

Table A.6: The most optimal model architecture, node/edge attribute type, and
normalization method that results in the lowest RMSE for target activity
regression. For each metric, the mean $\mu$ and standard deviation $\sigma$
are displayed across all 5 random seeds. “FP” denotes condensed fingerprint
featurization, and “One-hot” denotes one-hot encoding featurization.

Species | Model, Attribute, Normalization Type | Pearson r2 | MAE | RMSE | Spearman r | Kendall’s Tau | Loss
---|---|---|---|---|---|---|---
$\mu$ | $\sigma$ | $\mu$ | $\sigma$ | $\mu$ | $\sigma$ | $\mu$ | $\sigma$ | $\mu$ | $\sigma$ | $\mu$ | $\sigma$
_E. coli_ | Weave, FP, Quantile transform | 0.465 | 0.023 | 0.166 | 0.004 | 0.212 | 0.006 | 0.68 | 0.019 | 0.501 | 0.014 | 0.023 | 0.002
_S. aereus_ | GCN, One-hot, Quantile transform | 0.38 | 0.032 | 0.185 | 0.008 | 0.229 | 0.007 | 0.62 | 0.027 | 0.451 | 0.023 | 0.026 | 0.002

Figure A.37: Parity plots for _Escherichia coli_ validation dataset normalized
with standard scaler for regression of peptide target activity for each of
five model architectures (Attentive FP, Weave, MPNN, GCN, and GAT), all using
the same train-validation-test dataset split. A. Parity plots for fingerprint-
featurized experiments, with each graph displaying the top set of
hyperparameters. B. Parity plots for one-hot encoding-featurized experiments,
with each graph displaying the top set of hyperparameters. Figure A.38: Parity
plots for _Escherichia coli_ test dataset normalized with standard scaler for
regression of peptide target activity for each of five model architectures
(Attentive FP, Weave, MPNN, GCN, and GAT), all using the same train-
validation-test dataset split. A. Parity plots for fingerprint-featurized
experiments, with each graph displaying the top set of hyperparameters. B.
Parity plots for one-hot encoding-featurized experiments, with each graph
displaying the top set of hyperparameters. Figure A.39: Parity plots for
_Escherichia coli_ validation dataset normalized with quantile transform for
regression of peptide target activity for each of five model architectures
(Attentive FP, Weave, MPNN, GCN, and GAT). A. Parity plots for fingerprint-
featurized experiments, each graph for the top set of hyperparameters and the
same train-validation-test dataset split. B. Parity plots for one-hot
encoding-featurized experiments, each graph for the top set of hyperparameters
and the same train-validation-test dataset split. Figure A.40: Parity plots
for _Escherichia coli_ test dataset normalized with quantile transform for
regression of peptide target activity for each of five model architectures
(Attentive FP, Weave, MPNN, GCN, and GAT). A. Parity plots for fingerprint-
featurized experiments, each graph for the top set of hyperparameters with the
same train-validation-test dataset split across all plots. B. Parity plots for
one-hot encoding-featurized experiments, each graph for the top set of
hyperparameters with the same train-validation-test dataset split across all
plots. Figure A.41: Parity plots for _Staphylococcus aereus_ validation
dataset normalized with standard scaler for regression of peptide target
activity for each of five model architectures (Attentive FP, Weave, MPNN, GCN,
and GAT). A. Parity plots for fingerprint-featurized experiments, each graph
for the top set of hyperparameters with the same train-validation-test dataset
split across all plots. B. Parity plots for one-hot encoding-featurized
experiments, each graph for the top set of hyperparameters with the same
train-validation-test dataset split across all plots. Figure A.42: Parity
plots for _Staphylococcus aereus_ test dataset normalized with standard scaler
for regression of peptide target activity for each of five model architectures
(Attentive FP, Weave, MPNN, GCN, and GAT). A. Parity plots for fingerprint-
featurized experiments, each graph for the top set of hyperparameters with the
same train-validation-test dataset split across all plots. B. Parity plots for
one-hot encoding-featurized experiments, with each graph displaying the top
set of hyperparameters with the same train-validation-test dataset split
across all plots. Figure A.43: Parity plots for _Staphylococcus aereus_
validation dataset normalized with quantile transform for regression of
peptide target activity for each of five model architectures (Attentive FP,
Weave, MPNN, GCN, and GAT). A. Parity plots for fingerprint-featurized
experiments, each graph for the top set of hyperparameters with the same
train-validation-test dataset split across all plots. B. Parity plots for one-
hot encoding-featurized experiments, each graph for the top set of
hyperparameters with the same train-validation-test dataset split across all
plots. Figure A.44: Parity plots for _Staphylococcus aereus_ test dataset
normalized with quantile transform for regression of peptide target activity
for each of five model architectures (Attentive FP, Weave, MPNN, GCN, and
GAT). A. Parity plots for fingerprint-featurized experiments, each graph for
the top set of hyperparameters with the same train-validation-test dataset
split across all plots. B. Parity plots for one-hot encoding-featurized
experiments, each graph for the top set of hyperparameters with the same
train-validation-test dataset split across all plots.

##  Appendix H Attribution analysis

###  H.1 Selection of model architecture - attribution method

We used integrated gradients (IGs) (Sundararajan et al., 2017) and Input x
Grad (InpGrad) (Shrikumar et al., 2017) for the attribution analysis of the
GNNs - Weave, Attentive FP and MPNN. The model architecture selection was done
to have one of each type of architecture – Weave (graph convolution),
Attentive FP (graph attention), and MPNN (message passing).

The macromolecule is represented as a graph, $\mathcal{G}(V,E)$, where V
represents monomers at vertices/nodes, and E represents connecting bonds at
the edges.

IGs interpolate between the input graph and a baseline graph, where all
features are zero, and accumulate the gradient values for each node (Equation
1). The notation follows (Sanchez-Lengeling et al., 2020).

$\mathcal{G}_{A}=\left(\mathcal{G}-\mathcal{G}^{\prime}\right)\int_{\alpha=0}^{1}\frac{dy\left(\mathcal{G}^{\prime}+\alpha\left(\mathcal{G}-\mathcal{G}^{\prime}\right)\right)}{d\mathcal{G}}d\alpha$
(1)

InpGrad is the element-wise product of the input graph and the gradient.

$\mathcal{G}_{A}={\left(\frac{d\hat{y}}{d\mathcal{G}}\right)}^{T}\mathcal{G}$
(2)

For the attention-based GNN, AttentiveFP, in addition to IGs and InpGrad, we
evaluated attribution using attention weights, where the node attention
weights are obtained by averaging over the attention scores of the adjacent
nodes.

For each attribution method, we obtained the node weights by multiplying the
positive weights with the input fingerprint vectors -

$\mathbf{n}=\sum_{nodes}\mathcal{G}_{A}^{+}\;\mathcal{G},\;\;\;\;$ (3)

The node weights were normalized to the maximum node weight to obtain the
normalized weights.

$\mathbf{n_{norm}}=\frac{\mathbf{n}}{max(\mathbf{n})}$ (4)

Consistency of node weights, as defined in Sanchez-Lengeling et al. (2020),
was used for evaluation of different attribution methods and model
architectures.

The distribution of standard deviation for the node weights calculated using
IG and InpGrad have varying attribution consistency, both across attribution
methods and model architectures (Figure A.45). By visual inspection, Attentive
FP - IG has the smallest mean amongst all other distributions, thus the most
consistent model architecture - attribution. All further attribution analysis
has been done using node weights obtained from Attentive FP - IG.

Figure A.45: Distribution of normalized node weights for different model
architecture - attribution methods. Standard deviation of node weights
obtained from 25 different model implementations – top 5 hyperparameter sets
and 5 random seeds – were pooled across all immunogenic glycans to calculate
the distribution for each model architecture and attribution combination.
Table A.7: Top monomers contributing to glycan immunogenicity, arranged in
descending order of average node weight. Monomers which occur at least 20
times in the entire data set have been included in the table. The node weights
have been obtained using Attentive FP – IG combination.

Monomer | Mean Node Weight
---|---
NeuNGc | 0.73
Xyl | 0.65
Fuc | 0.65
Glc | 0.64
ManNAc | 0.64
Galf | 0.63
NeuNAc | 0.62
FucNAc | 0.61
QuiNAc | 0.61

###  H.2 Visualization of key substructures

To visualize the responsible substructures in the monomers, we used
$\mathcal{G}_{A}$ in Equation 1, and multiplied the weights of the with the
respective monomer fingerprint. This approach resulted in a weights vector
with the same size as the fingerprint, with the most positively to the most
negatively influencing substructure for the prediction. Using RDKit, we
visualized the chemical substructures at different fingerprint indices and
mapped it to the weights.

###  H.3 Ablation analysis

We performed ablation analysis to check if the monomer importance or node
weights changed when we assigned a zero tensor to one or more nodes (Figure
A.46). For node- and monomer type-specific ablation, we observed that the
importance of non-ablated nodes remained consistent.

Figure A.46: Ablation analysis for immunogenic glycans shows that the
importance of non-ablated monomers remains consistent, irrespective of node-
and monomer-specific ablation. A. Glycan without any ablations is visualized
as the baseline case. B. Monomer-specific ablation shows that ablation of a
particular monomer has a negligible impact on the importance/weight of other
monomers. Ablation was done by assigning zero tensors to nodes of respective
monomers. This experiment was done by subsequently mutating all nodes
corresponding to Gal, Xyl, Glc, and Fuc. C. Node-specific ablation shows that
other nodes of a particular monomer have a negligible change when one of the
nodes is assigned a zero tensor. This experiment was done by mutating all 4
nodes corresponding to Xyl, one at a time.

## Supplementary Tables

Supplementary Table 1: Validation and test metrics for the top performing sets
of model hyperparameters averaged across model implementations with 5 random
seeds for glycan classification, as well as test metrics averaged across 25
model implementations (5 hyperparameter sets and 5 random seeds)

Supplementary Table 2: GNN model hyperparameter sets for the top performing
glycan classification models

Supplementary Table 3: Validation and test metrics for the top performing sets
of model hyperparameters averaged across model implementations with 5 random
seeds for peptide regression, as well as test metrics averaged across 25 model
implementations (5 hyperparameter sets and 5 random seeds)

Supplementary Table 4: GNN model hyperparameter sets for the top performing
peptide regression models