# Classifying Online Dating Profiles on Tinder using FaceNet Facial Embeddings

###### Abstract

A method to produce personalized classification models to automatically review
online dating profiles on Tinder, based on the user’s historical preference,
is proposed. The method takes advantage of a FaceNet facial classification
model to extract features which may be related to facial attractiveness. The
embeddings from a FaceNet model were used as the features to describe an
individual’s face. A user reviewed 8,545 online dating profiles. For each
reviewed online dating profile, a feature set was constructed from the profile
images which contained just one face. Two approaches are presented to go from
the set of features for each face to a set of profile features. A simple
logistic regression trained on the embeddings from just 20 profiles could
obtain a 65% validation accuracy. A point of diminishing marginal returns was
identified to occur around 80 profiles, at which the model accuracy of 73%
would only improve marginally after reviewing a significant number of
additional profiles.

Index Terms—  facial classification, facial attractiveness, online dating,
classifying dating profiles

##  1 Introduction

Online dating has become a commonplace in today’s society. The top grossing
iOS application in September 2017 was an online dating service named Tinder.
Users of online dating services are generally expected to spend a significant
amount of time filtering through the profiles of potential partners. This
study investigates whether a pattern in facial features can be used to filter
online dating profiles on Tinder.

Tinder was selected for this study because of the popularity of the
application [1]. The mobile dating application allows users to browse dating
profiles of nearby singles. Users are presented with a single profile at a
time. At first glance, the user can see the profile picture, first name, and
age of an individual on Tinder. A user must then decide whether they like or
dislike the profile on the spot. The user can view additional information such
as an optional biography or extra pictures, but no new profiles can be
reviewed until a decision is made. When two individuals have liked each other,
they are presented with a notification and the opportunity to message each
other.

A custom application was developed to interface with Tinder. The intention of
the application was to allow a user to like and dislike profiles while
recording everything about the profiles. The profile information contains the
mandatory requirements from tinder, which include the profile’s name, age, and
pictures. Every Tinder profile includes at least one image. Additionally, it
is optional to include a biography, current job, and current school. All of
this information is stored along with the like or dislike verdict in a
database.

A method for obtaining personalized classification models from a user’s
historical preferences is presented. The method takes advantage of recent
advancements in computer vision related to facial detection, facial
classification, and facial attractiveness. The results of a personalized
classification model for a user’s reviewed Tinder profiles are then presented.

##  2 Related Works

Related work focuses on literature regarding modeling facial attractiveness.
Attractiveness is known to play a key role in online dating profiles [1], [2].

There has been substantial work in literature to predict facial attractiveness
with computer models. Traditionally, facial attractiveness has been modeled
with some sort of eigenface transformation. Eisenthal et al. [3] were able to
find a correlation of facial attractiveness using an eigenface transformation
on a small data set of 92 images. These eigenface transformations are run
through a principle component analysis to determine the important facial
features [4], [5].

Convolutional neural networks (CNN) have become popular for image processing
in recent years. Gray et al. [6] showed that a CNN trained to predict facial
attractiveness could extract more meaningful features than an eigenface
transformation. Consequently, facial features related to attractiveness could
be extracted automatically without identifying landmarks using the CNN.

The recent works of Rothe et al. [7] demonstrated that a CNN can be trained on
a large data set to predict facial attractiveness. The training set contained
the faces of women who were liked or disliked by male participants.
Personalized predictions were made based on the historical preferences of
8,560 heterosexual male users. The model only considered the first profile
image from each profile. The results of the personalized predictions were
impressive, with a mean accuracy of 75% for the male users.

Zang et al. [8] used a pre-trained VGG-Face CNN descriptor to predict the
facial attractiveness of male and female users for an online dating site. They
demonstrated that features from a facial classification model could be used to
predict facial attractiveness. The model in Rothe et al. [7] was trained on a
data set related to facial attractiveness, while the model used in Zang et al.
[8] was trained on faces for the purpose of classification. These studies
highlight that a large scale facial classification model is useful to predict
facial attractiveness. The VGG-Face CNN used was created by Parkhi et al. [9]
and scores an impressive 98.95% accuracy on the Labeled Faces in the Wild
(LFW) database [10].

These works focused solely on rating individual photos, but have not
progressed to a usable model that likes or dislikes complete online dating
profiles. The work presented in this paper strives to close this gap.

##  3 Our Methodology

The methodology proposed here attempts to classify an online dating profile as
either a like or dislike. Two different approaches are proposed to consolidate
multiple facial features from the images in a profile into a single vector of
features that describes the profile. Like the related works of [7], [8], the
last layer of a CNN was used as the facial features for each face. A new
implementation of the FaceNet classification model first described by Schroff
et al. [11] is used with a slightly higher LFW score than used in Zang et al.
[8].

The detection of profile images that contain only one face per image was
automated using computer vision techniques. These faces were fed into a
FaceNet model to extract the facial features as embeddings. A set of
embeddings for reviewed online dating profiles was used to train a
personalized classification model.

The major assumptions of the purposed method are as follows: 1) An online
dating profile can be reviewed using only the profile images; 2) The face of
the individual profile can be found from the profile pictures that contain
only one face per image; 3) Images with more than one face in an online dating
profile can be ignored. Profiles that can’t be identified to a single face can
be rejected; 4) A pattern exists in the faces of individuals who were liked or
disliked by a user reviewing online dating profiles; 5) A trained FaceNet
model can be evaluated on new faces to extract the facial features of the
individuals.

###  3.1 FaceNet implementation

A Python library called facenet was used to calculate the facial embeddings of
the dating profile pictures. These embeddings are from the last layer of a
CNN, and can be thought of as the unique features that describe an
individual’s face. The facenet library was created by Sandberg as a TensorFlow
implementation of the FaceNet paper by Schroff et al. [11], with inspirations
from [9, 12, 13]. The library uses the MIT license and is available online at
https://github.com/davidsandberg/facenet.

There are pre-trained facenet models available online. The models have been
validated on the LFW database [10]. The current best model has a LFW accuracy
of 99.2% and was trained as a classification model on a subset of the MS-
Celeb-1M database [14]. The model’s architecture follows the Inception-
ResNet-v1 network as described by Szegedy et al. [15]. The facenet library
includes an implementation of Multitask CNN (MTCNN) by Zhang et al. [16] to
detect facial landmarks, which was used to create training faces as 182x182
pixel images from the MS-Celeb-1M database.

The facenet model turns a color image of a face into a vector of 128 floating
point numbers. These 128 embeddings can be used as features for classification
or clustering [11]. The facenet library includes a script to calculate the
embeddings from images of faces using a pre-trained model.

###  3.2 Classification methodology

Classification models were determined for two different approaches on the
embeddings. One approach considered all of the embeddings, from the images
containing just one face, in the dating profile. These embeddings were used to
describe the entire profile. The other approach rather considered the average
embedding values across the images. Again, only images containing exactly one
face were considered.

The first approach used the 128 embeddings from each image as the features of
the profile. The embeddings from the images of the profiles can be described
as the vectors of

$\displaystyle\bm{i}_{1}$ $\displaystyle=[x_{1},x_{2},\cdots,x_{128}]$ (1)
$\displaystyle\bm{i}_{2}$ $\displaystyle=[y_{1},y_{2},\cdots,y_{128}]$ (2)
$\displaystyle\vdots$ $\displaystyle=\vdots$ (3) $\displaystyle\bm{i}_{n}$
$\displaystyle=[z_{1},z_{2},\cdots,z_{128}]$ (4)

for $n$ number of profile images. Then a single vector of embeddings can be
constructed for the profile as

$\displaystyle\bm{i}_{p}$
$\displaystyle=[\bm{i}_{1},\bm{i}_{2},\cdots,\bm{i}_{n}]$ (5)

where $\bm{i}_{p}$ is a vector containing $128n$ values.

The second approach considered the average embedding value of the facial
images. Thus a profile with one facial image would have 128 unique embeddings.
A profile could be described as

$\displaystyle\bm{i}_{1}$ $\displaystyle=[x_{1},x_{2},\cdots,x_{128}]$ (6)
$\displaystyle\bm{i}_{2}$ $\displaystyle=[y_{1},y_{2},\cdots,y_{128}]$ (7)
$\displaystyle\vdots$ $\displaystyle=\vdots$ (8) $\displaystyle\bm{i}_{f}$
$\displaystyle=[z_{1},z_{2},\cdots,z_{128}]$ (9)

where $\bm{i}_{f}$ is the vector of embedding from the $f$ image in the
profile. Then an average embeddings could be calculated as

$\displaystyle\bm{i}_{\text{avg}}$
$\displaystyle=\begin{bmatrix}\frac{x_{1}+y_{1}+\cdots+z_{1}}{f}\\\
\frac{x_{2}+y_{2}+\cdots+z_{2}}{f}\\\ \vdots\\\
\frac{x_{128}+y_{128}+\cdots+z_{128}}{f}\\\ \end{bmatrix}$ (10)

where $\bm{i}_{\text{avg}}$ is a vector with the same size as the number of
embeddings calculated.

Calculating the facial embeddings from a user’s reviewed online dating
profiles is computationally inexpensive, as the calculation is simply a
function evaluation on a pre-trained CNN. Then, classification models were
trained using either $\bm{i}_{p}$ or $\bm{i}_{\text{avg}}$ as the input
features. Personalized classification models could be constructed based on the
preference from an individual’s historically reviewed online dating profiles.

##  4 Experimental Results

A heterosexual male used the custom application with the intention of finding
a romantic partner. The reviewing of tinder profiles went on for a month, but
stopped early because the user found a girlfriend in the process. It may be
important to mention that males may have different online dating tendencies
than females [1, 2]. The user took about one hour to review 100 profiles. In
the end, a data set was created which reviewed 8,545 tinder profiles. The user
liked a total of 2,411 profiles. Additionally, the data set contains 38,218
images from the profiles browsed. Each image has a resolution of 640x640
pixels (px).

The results were split into two categories. The first subsection presents the
results of the data set after pre-processing was performed. The data set was
transitioned from complete online dating profiles to a data set of faces for
each profile. The faces were then run through a FaceNet model to extract the
embeddings for each face. The second section then presents the results of
classifying these embeddings for the two proposed input dimensions.

###  4.1 Data set after pre-processing

The MTCNN described by [16] was used to detect and box faces from the 640x640
px profile images. Faces were selected with a minimum size of 60x60 px and a
threshold of $0.8$. Profile images that contained just one face were extracted
and re-sized. A profile that did not contain a single image with only one
face, was immediately removed. There were 24,486 images that contained only
one face in the image (according to the MTCNN). Fortunately 8,130 profiles of
the 8,545 reviewed (or 95.1%) contained at least one uniquely identifiable
face. The images containing just one face were cropped to 182x182 px images
with a margin of 44 px around the face. A face at the minimum size was
enlarged, while larger faces were reduced in size.

The MTCNN results were impressive, despite the substantial amount of noise in
the images. Noise includes everything from sunglasses, hats, and scarfs to
Snapchat filters. For example, a particular popular Snapchat filter applies
the ears, nose, and mouth of a dog to an individual’s face. The MTCNN appeared
to work well despite the noise in the data. There was a limited number of
false positives, of which a few are presented in Fig. 1. The false positives
were not removed from the training set, as the noise they provide may be
useful to construct a robust classifier. The true rate of false positives and
false negatives was not studied, as the locations of faces in the original
38,218 images were not recorded.

Fig. 1: Examples of false positives which the MTCNN identified as human faces.

The embeddings were calculated for the faces from the 8,130 profiles using the
FaceNet implementation described. The average profile reviewed had 3.01 images
of a single face, with a standard deviation of 1.34. Ten was the maximum
number of images for a profile in the new data set. Thus $\bm{i}_{p}$ was a
vector of $128\times 10$ in length. Profiles with fewer than ten images would
have zeros in place of the missing images. Essentially a profile with just one
facial image would have 128 unique embeddings and 1,152 zeros, a profile with
two facial images would have 256 unique embeddings and 1,024 zeros, and so
forth. The other input feature $\bm{i}_{\text{avg}}$ was calculated for each
profile. The supplementary material includes the two input dimensions
($\bm{i}_{p}$ and $\bm{i}_{\text{avg}}$) with binary labels to show whether
the profile was either liked or disliked.

###  4.2 Classification models

In order to build a reasonable classification model, it was important to
demonstrate how many profiles were required to be reviewed. Classification
models were trained using various fractions of the entire data, ranging from
0.125% to 95% of the 8,130 profiles. At the low end, just 10 profiles were
used to train the classification model, while the remaining 8,120 profiles
were used to validate the trained classification model. On the other spectrum,
classification models were trained using 7,723 profiles and validated on 407
profiles.

The classification models were scored on accuracy, specifically the number of
correctly classified labels over the number of profiles. The training accuracy
refers to the accuracy in the training set, while the validation accuracy
refers to the accuracy in the test set.

This data set suffers from a class imbalance, as only 28% of the total Tinder
profiles reviewed were liked. The classification models were trained assuming
a balanced class. A balanced class indicates that each profile considered had
the same weight, regardless of whether the profile was liked or disliked. The
class weight can be user dependent, as some users would value correctly liking
profiles more than incorrectly disliking profiles.

A like accuracy was introduced to represent the number of correctly labeled
liked profiles out of the total number of liked profiles in the test set.
Complementary, a dislike accuracy was used to measure the disliked profiles
predicted correctly out of the total number of disliked profiles in the test
set. A model that disliked every single profile, would have a 72% validation
accuracy, a 100% dislike accuracy, but a 0% like accuracy. The like accuracy
is the true positive rate (or recall), while the dislike accuracy is the true
negative rate (or specificity).

Various classification models were fit using either $\bm{i}_{p}$ or
$\bm{i}_{\text{avg}}$ as the input. Scikit-learn was used to fit logistic
regression and support vector machines (SVM) [17], while Keras and TensorFlow
were used to fit various neural networks [18, 19] to the embeddings.

The receiver operating characteristic (ROC) for logistic regression (Log),
neural network (NN), and SVM using radial basis function (RBF) are presented
in Fig. 2. Two different layer configurations of neural networks are presented
for each input dimension as NN 1 and NN 2. Additionally, the area under curve
(AUC) for each classification model is presented. The complete input dimension
feature of $\bm{i}_{p}$ did not appear to offer any advantages over
$\bm{i}_{\text{avg}}$ when considering AUC. A neural network had the best AUC
score of 0.83, but it was only slightly better than a logistic regression with
an AUC score of 0.82. This ROC study was performed using a random 10:1
train:test split (training on 7,317 and validation on 813 profiles).

Fig. 2: Receiver operating characteristic (ROC) and area under curve (AUC) for
various classification models using a 10:1 train:test split.

Since the AUC scores were comparable, the remaining results only consider
classification models fit to $\bm{i}_{\text{avg}}$. Models were fit using
various train-to-test ratios. The train:test split was performed at random;
however each model used the same random state for a given number of training
profiles. The ratio of likes to dislikes was not preserved in the random
splits. The training accuracy of the models is presented in Fig. 3 and the
validation accuracy for these models is presented in Fig. 4. The first data
point represents a training size of 10 profiles and a validation size of 8,120
profiles. The last data point uses 7,723 training profiles and validation on
407 profiles (a 20:1 split). The logistic regression model (Log) and neural
network (NN 2) converge to a comparable training accuracy of 0.75.
Impressively, a model can have a validation accuracy greater than 0.5 after
being trained on just 20 profiles. A reasonable model with a validation
accuracy near 0.7 was trained on just 40 profiles.

Fig. 3: Accuracy of the classification models evaluated on the training data
(number of correctly classified training profiles / total number of training
profiles) vs the number of training profiles. Fig. 4: Validation accuracy of
classification models (number of correctly classified test profiles / total
number of test profiles) as a function of the number of training profiles.

The validation accuracy can be misleading because there were significantly
more disliked profiles than liked profiles. Thus it is important to consider
the true positive and true negative rates to asses the model quality. The like
accuracy (or true positive rate) of the models is presented in Fig. 5, and the
dislike accuracy (or true negative rate) is presented in Fig. 6. It can be
noted that the neural networks had a like bias, while the SVM has a dislike
bias. The dislike bias resulted in a slightly higher validation accuracy since
there were more profiles disliked than liked. Again, after 20 profiles, a
reasonable classification model can be constructed. The logistic regression
model trained on 20 profiles had a like accuracy of 0.7 and a dislike accuracy
of 0.6.

Fig. 5: Like accuracy (true positive rate) of classification models as a
function of the number of training profiles. Fig. 6: Dislike accuracy (true
negative rate) of classification models as a function of the number of
training profiles.

The results presented thus far could be the result of the random split.
Another study is presented to better understand how the number of reviewed
profiles may influence the personalized classification model. Random splits
were performed 10,000 times for training sizes of 10, 20, 40, 81, and 406
profiles. The random split required at least one profile from each class (like
and dislike) in order to construct a classification model. A unique logistic
regression model was fit for each split and validated on the remaining test
data. Again, the training of 10 profiles was validated on 8,120 profiles, and
so forth for the other training sizes. The resulting validation accuracies
followed a skew-normal distribution, and probability density functions (PDF)
were calculated for each training size. The resulting PDFs are presented in
Fig. 7 and compared to the PDF of a completely random classifier. The
validation accuracy for a completely random classifier was simulated 10,000
times and followed a normal distribution. The results demonstrate that
training on just 10 dating profiles offers a significant advantage over a
random classifier. The variance associated with a model’s validation accuracy
is shown to reduce with the number of trained profiles. This reduction of
variance is significant when going from training on 81 profiles to 406
profiles.

Fig. 7: Probability density functions (PDF) for validation accuracy of
classifiers trained on either 10, 20, 40, 81, and 406 profiles. The PDF from a
completely random classifier is shown as a baseline comparison.

A Python script is included in the supplementary material to calculate the the
results presented here for the logistic regression model using either
$\bm{i}_{p}$ or $\bm{i}_{\text{avg}}$ as the input dimension.

##  5 Conclusion

A method was presented to build personalized classification models for online
dating profiles based on a user’s historical preference. The method could be
used to improve the user experience of online dating by reducing the time
required to filter profiles. A custom data set was collected which reviewed
over 8,000 Tinder profiles. Profile images containing just one face were run
through a FaceNet model to extract the unique features as embeddings. Two
different approaches were presented to combine these features from faces in a
profile, to a unique vector representing the features of that profile. A
classification model was then constructed either considering a 128 or 1280
input dimension. A simple logistic regression model was shown to find an
accuracy greater than 60% after being trained on just 20 profiles. The
classification methodology continuously improves as more online dating
profiles are reviewed. Additionally it was demonstrated that a classification
model trained on just 10 profiles would, on average, have a much higher
validation accuracy than a random classifier.

A Python command line application called tindetheus has been released to
reproduce the methodology presented in this paper. The application has three
major functions: 1) Build a data set as a user browses Tinder. 2) Train a
classification model to the data set. 3) Use the trained model to
automatically like new Tinder profiles. Tindetheus is available online at
https://github.com/cjekel/tindetheus/ or
https://pypi.python.org/pypi/tindetheus .

## References

  * Tyson et al. [2016] G. Tyson, V. C. Perta, H. Haddadi, and M. C. Seto, “A first look at user activity on tinder,” in _2016 IEEE/ACM International Conference on Advances in Social Networks Analysis and Mining (ASONAM)_ , aug 2016, pp. 461–466.
  * Abramova et al. [2016] O. Abramova, A. Baumann, H. Krasnova, and P. Buxmann, “Gender differences in online dating: What do we know so far? A systematic literature review,” in _Proceedings of the Annual Hawaii International Conference on System Sciences_ , vol. 2016-March, jan 2016, pp. 3858–3867.
  * Eisenthal et al. [2006] Y. Eisenthal, G. Dror, and E. Ruppin, “Facial attractiveness: Beauty and the machine,” _Neural Computation_ , vol. 18, no. 1, pp. 119–142, Jan 2006.
  * Mu [2013] Y. Mu, “Computational facial attractiveness prediction by aesthetics-aware features,” _Neurocomputing_ , vol. 99, no. Supplement C, pp. 59 – 64, 2013\. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S092523121200495X
  * Chua et al. [2013] M. Chua, Y. Akimoto, H. Aguirre, and K. Tanaka, “Asian female face classification incorporating personal attractive preference,” in _2013 International Symposium on Intelligent Signal Processing and Communication Systems_ , Nov 2013, pp. 413–418.
  * Gray et al. [2010] D. Gray, K. Yu, W. Xu, and Y. Gong, _Predicting Facial Beauty without Landmarks_. Berlin, Heidelberg: Springer Berlin Heidelberg, 2010, pp. 434–447. [Online]. Available: https://doi.org/10.1007/978-3-642-15567-3_32
  * Rothe et al. [2016] R. Rothe, R. Timofte, and L. V. Gool, “Some like it hot!; visual guidance for preference prediction,” in _2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , June 2016, pp. 5553–5561.
  * Zang et al. [2017] X. Zang, T. Yamasaki, K. Aizawa, T. Nakamoto, E. Kuwabara, S. Egami, and Y. Fuchida, “Prediction of users’ facial attractiveness on an online dating website,” in _2017 IEEE International Conference on Multimedia Expo Workshops (ICMEW)_ , jul 2017, pp. 255–260.
  * Parkhi et al. [2015] O. M. Parkhi, A. Vedaldi, and A. Zisserman, “Deep Face Recognition,” in _British Machine Vision Conference_ , 2015.
  * Huang et al. [2007] G. B. Huang, M. Ramesh, T. Berg, and E. Learned-Miller, “Labeled faces in the wild: A database for studying face recognition in unconstrained environments,” University of Massachusetts, Amherst, Tech. Rep. 07-49, oct 2007\.
  * Schroff et al. [2015] F. Schroff, D. Kalenichenko, and J. Philbin, “FaceNet: A Unified Embedding for Face Recognition and Clustering,” in _The IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , jun 2015.
  * Wen et al. [2016] Y. Wen, K. Zhang, Z. Li, and Y. Qiao, “A Discriminative Feature Learning Approach for Deep Face Recognition,” in _Computer Vision – ECCV 2016: 14th European Conference, Amsterdam, The Netherlands, October 11–14, 2016, Proceedings, Part VII_ , B. Leibe, J. Matas, N. Sebe, and M. Welling, Eds. Cham: Springer International Publishing, 2016, pp. 499–515. [Online]. Available: https://doi.org/10.1007/978-3-319-46478-7_31
  * Amos et al. [2016] B. Amos, B. Ludwiczuk, and M. Satyanarayanan, “OpenFace: A general-purpose face recognition library with mobile applications,” CMU-CS-16-118, CMU School of Computer Science, Tech. Rep., 2016.
  * Guo et al. [2016] Y. Guo, L. Zhang, Y. Hu, X. He, and J. Gao, “MS-Celeb-1M: A dataset and benchmark for large scale face recognition,” in _European Conference on Computer Vision_. Springer, 2016.
  * Szegedy et al. [2016] C. Szegedy, S. Ioffe, V. Vanhoucke, and A. Alemi, “Inception-v4, inception-resnet and the impact of residual connections on learning,” in _Proceedings of the Thirty-First AAAI Conference on Artificial Intelligence, February 4-9, 2017, San Francisco, California, USA._ , S. P. Singh and S. Markovitch, Eds. AAAI Press, 2016, pp. 4278–4284. [Online]. Available: http://arxiv.org/abs/1602.07261
  * Zhang et al. [2016] K. Zhang, Z. Zhang, Z. Li, and Y. Qiao, “Joint Face Detection and Alignment Using Multitask Cascaded Convolutional Networks,” _IEEE Signal Processing Letters_ , vol. 23, no. 10, pp. 1499–1503, oct 2016.
  * Pedregosa et al. [2011] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay, “Scikit-learn: Machine Learning in Python,” _Journal of Machine Learning Research_ , vol. 12, pp. 2825–2830, 2011.
  * Chollet and Others [2015] F. Chollet and Others, “Keras,” https://github.com/fchollet/keras, 2015\.
  * Abadi and et al. [2015] M. Abadi and et al., “TensorFlow: Large-scale machine learning on heterogeneous systems,” 2015, software available from tensorflow.org. [Online]. Available: https://www.tensorflow.org/


# Towards complete representation of bacterial contents in metagenomic samples

X.F.Xiaowen Feng  H.L.Heng Li  Department of Data Sciences, Dana-Farber Cancer
Institute, Boston, US  Department of Biomedical Informatics, Harvard Medical
School, Boston, US

###### Abstract

#### Background:

In the metagenome assembly of a microbiome community, we may think abundant
species would be easier to assemble due to their deeper coverage. However,
this conjucture is rarely tested. We often do not know how many abundant
species we are missing and do not have an approach to recover these species.

#### Results:

Here we proposed k-mer based and 16S RNA based methods to measure the
completeness of metagenome assembly. We showed that even with PacBio High-
Fidelity (HiFi) reads, abundant species are often not assembled as high strain
diversity may lead to fragmented contigs. We developed a novel algorithm to
recover abundant metagenome-assembled genomes (MAGs) by identifying circular
assembly subgraphs. Our algorithm is reference-free and complement to standard
metagenome binning. Evaluated on 14 real datasets, it rescued many abundant
species that would be missing with existing methods.

#### Conclusions:

Our work stresses the importance of metagenome completeness which is often
overlooked before. Our algorithm generates more circular MAGs and moves a step
closer to the complete representation of microbiome communities.

#### Keywords:

Metagenome, Binning, Metagenome-assembled genomes, assembly completeness

Research

## Background

De novo metagenome assembly has been an allusive promise of unbiased and
comprehensive snapshot of microbial communities of interest, independent to
isolation and cultivation [1, 2, 3]. Neither of the two aspects has been close
to realization. First, most past metagenome sequencing projects were based on
short read sequencing and produced short contigs of tens of kilobases (kb)
long, which need to be clustered to form metagenome-assembled genomes (MAGs).
Most de novogenome-complete [4, 5] MAGs still contain an average of 87
assembly gaps with median length about 1.3kb. We manually checked some and
found that these gaps either have no presence in the BLAST nr/nt database
using BLASTn, or were homologous to shared genes such as ribosomal RNA (rRNA)
operons. Second, MAGs are rarely checked for their representation-
completeness. Studies often assumed that sufficiently abundant or the most
abundant species will be reconstructed [6, 7, 8]. There is no such guarantee
despite efforts to approve it [9]. A distinct species with low coverage can be
easier to recover than abundant but highly similar species or strains.
Horizontal gene transfer (HGT) or large duplication events are even more
difficult to resolve. One major obstacle for improving the situation was that
16S rRNA sequences, which is a proxy of species definition [10, 11, 12], were
no easier to assemble than the whole genome [13] and will require amplicon
sequencing [14, 15, 16]. Therefore it is hard to cross validate between gene-
based composition inference and assembly [8, 17, 18].

Pooling short read MAGs to form a complete reference catalogs [19, 18, 20] is
also flawed. The human gut microbiome is the most genome-sequenced metagenome
context. SRA has near 83k results assigned to this category as of July 2022.
However, recent large-scale studies still found 42% of their quality MAGs
missing from major public repositories [18]. For gut microbiome of other
species, this number can be as high as 86% [21]. Moreover, taxonomic profiling
of 1004 faecal samples from TwinsUK registry found each species was observed
in a median of 2.7% samples, with 12% of species being sample-specific and 50%
species found in less than 1% samples [22], demonstrating vast diversities. To
make the situation more complicated, there is no consensus on the optimal
library size. Most sequencing projects prioritized sample size over high-
coverage, perhaps due to the observation that MAG yield per gigabases
sequenced per sample peaks at the lower end of library sizes (Figure S1-2).
Some of the most deeply sequenced [23, 24] or unusual [25] datasets were
analyzed unassembled. The endgame of the reference MAG accumulation process is
unclear. On the other hand, the drawbacks of short read MAGs are still
relevant even if MAGs are just treated as bags of genes [26, 27], or used for
composition inference [11, 28].

Metagenome assembly using accurate long reads, such as Pacbio HiFi reads, can
now recover haplotype-resolved near-complete closed MAGs and separate species
with sequence divergence as low as 1% [29, 30, 31]. This somewhat resolves the
first issue of short read MAGs, although the quality evaluation remains
reference-dependent. Binning re-emerged to be problemetic. Traditional binners
tend to collected contigs from multiple haplotypes, despite they are well
capable of creating high quality short read MAGs containing tens of hundreds
of short contigs. In HiFi assemblies, only when the binner pulls less than ten
contigs would a bin likely to pass quality checks, basically relying on the
recruitment of a few long contigs and nothing more. The second issue about
representation-completeness is an open question. However, with circular
contigs no longer relying on binning and full length 16S rRNA sequences
readily available from the HiFi reads, we think it is possible to approach it
from both k-mer specturm, leanring from evaluations in single sample large
genome assemblies, and 16S-based species-level operational taxonomic units
(OTUs). It is no less important about learning what is absent from the
assemblies than examining what has been recovered in order to properly compare
across samples and different microbiotas, and further, to improve the assembly
method.

This manuscript has three major components. We first introduce a graph
topology-based, reference-free postprocessing step for the hifiasm-meta [31]
assembler. The core is a circle-finding depth-first search. Although simple,
this could replace about half of regularly binned MAGs with pseudo circular
paths of identical or better quality, as well as adding in MAGs that are not
recoverable by regular binners. There is no external data needed, and the end
results were very close to the almost-best-possible MAG recovery (i.e. if
assisted by checkM). We then propose two evaluations as mentioned above and
apply them to the HiFi MAGs. We demonstrate that in the best cases, near-
complete MAGs alone could offer proper sample representations comparable to
that of a comprehensive reference catalog and delegate almost 90% abundant
OTUs (more than 30 16S copies, using 99% boundary).

## Results

### Topology-based cycle-finding on the contig graph as a binning method and
its merging with traditional binning

Established binners were developed for binning short read metagenome
assemblies. The process relies on tetranucleotide profiles and coverage
estimates [32, 33], and sometimes consults assembly graph for further
refinements [34]. Single-copy marker genes and other prior knowledge can also
be exploited [35, 36], though such binning complicates the final evaluations,
as assessing assemblies of empirical datasets heavily relies on checkM [4] or
similar methods that also utilize marker gene sets and phylogenetic placement.

In short read assemblies, the binners are capable of pooling hundreds of short
contigs to form very high quality MAGs. However, they do not perform as good
in HiFi assemblies. We found most near-complete or high-quality MAGs generated
by binners contain less than 10 contigs. The completeness of these bins mostly
come from only one or two long linear contig(s). While this does not defeat
the purpose of binning on HiFi assemblies, haplotypes that do not assemble
well would remain missing. Additionally, long contigs of hundreds or more
kilobases from closely related haplotypes are prone to be clustered together,
which results in complete but highly contaminated bins. The first issue was
especially hard to be improved. The second can be somewhat mitigated by
arbitrarily splitting bins based on intra-bin sequence similarity. We noticed
that the assembly graph of hifiasm-meta was simple enough to provide binning
hints when the assembly was visually [37] reasonably not bad (Figure 1). We
propose a simple topology-based assembly postprocessing in complement to
binning to achieve a reasonable complete representation of the library. The
method also has bonus contig ordering and circularization. We implemented this
approach in the hifiasm-meta assembler. The other two capable assemblers,
metaFlye and HiCanu, produces more tangled graphs or does not produce a graph,
respectively.

Figure 1: (a) Sankey plot showing flow of reads and contigs of sheep-gut-1b.
Left: reads to contig mappability (“minimap2 -cxmap-hifi”). Middle: contig
binning categories. Right: MAGs quality categories. Group heights normalized
by counts, except the left-most side which is normalized by base pairs. The
heights of “not aligned” and the “dismissed” have no meaning. (b) Bandage [37]
plot of sheep-gut-1b primary assembly: (top-left) primary assembly graph,
(top-right) circular path candidates to rescue marked in blue, (bottom-right)
high-quality-or-above rescued circles marked in green, and (bottom-left) best
possible high-quality-or-above MAGs marked in red. Note that the size of
circle-shaped subgraph that consists of numerous short contigs was exaggerated
by the graph layout method.

We use two assumptions: a bacterial chromosome is expected to be circular, and
its length to fall into the range of 1Mb-8Mb. We work on the primary assembly
graph. We first collect cycles of roughly 1Mb-10Mb long on the contig graph
via depth-first search (DFS). To balance redundancy and the risk of missing
promising cycles, the DFS always try to use a contig that has not been a node
in any of the collected cycles. The criteria are permissive and only
terminates the search if 100-200 nodes have been used. If a search extends too
long, we halve the stack and resume from halfway, in order to avoid missing
cycles due to the arbitrary ordering of contigs.

We then remove duplicated cycles using Mash distance [38]. Mash distance
correlates well with ANI in the range of 90-100% and we use bottom 1000 hashes
with k-mer size of 21 per recommendation in the Mash paper. If the two cycles
being compared differ in length for more than 1Mb, we do not drop any of them
regardless of Mash distance. We set the species boundary at 95% whole-genome
Average Nucleotide Identity (ANI), which was derived from DNA-DNA
hybridization values and later reconfirmed by large-scale studies [39, 40]. We
also tried 97% as a slightly relaxed threshold, and 99% which was the highest
possible ANI for most pairs of long circular contigs produced in HiFi
assemblies (hifiasm-meta, Hicanu or metaFlye). 97% threshold produced similar
completeness evaluation results to those of 95%. 99% threshold resulted in
more duplications, although also slightly improved the overall assembly
completeness. Some contigs might be used more than once. We see little such
duplication in the final result (i.e. checkM-passed MAGs). Note that the
comparison is restricted within the collection of rescued circles candidates.
A cycle that is very similar to a long circular contig would not be discarded
because of so.

We run MetaBAT2 on raw contigs independent of circle rescue and then merge
MetaBAT2 bins with rescued circles. We collect all rescued circles (set1) and
remember the contigs used by them (setb). Next, MetaBAT2 bins that are at
least 500kb, no more than 10Mb and do not only contain a single $>$1Mb
circular contig were examined. If a bin has more than 1Mb of its contents
coming from setb, or more than 10 of its contigs coming from setb, the bin is
rejected. Otherwise we accept it into set2 and update setb with its contigs.
We then collect linear contigs that are longer than 1Mb and not in setb into
set3. All circular contigs longer than 1Mb form set4. The binning outcome
(“merged MAGs” or “merged bins”) is the union of set1, 2, 3 and 4.

### MAG quality evaluations

We first compared MAG quality brackets from rescued circles, MetaBAT2 bins and
vamb bins (Table 1). There are 16 HiFi libraries available, all except one are
gut materials (Table S1, S4). We follow the minimum information criteria about
a single amplified genome (MISAG) convention [5] using checkM: “Near-complete”
means $\geq 90\%$ completeness and $<$5% contamination. “High-quality” means
$\geq 70\%$ completeness, $<$10% contamination but does not qualify for near-
complete. “Medium-qualify” means $\geq 50\%$ completeness, $\geq 50$ quality
score but does not qualify for the above two. Quality score of a MAG is
defined as “${\rm completeness}-5\times{\rm contamination}$”. In terms of the
yield in the near-complete quality bracket, the proposed topology-based rescue
performed occasionally comparable to vamb, but was consistently worse than
MetaBAT2 except for in env-digester-1. This was because binners were able to
pull together long contigs that are not connected on the assembly graph to
form a checkM-passing bin, while the rescue heuristic relies on graph
traversal. Below we demonstrate the rescuing and the merge with regular
MetaBAT2 bins using sheep-gut-1b, a diverse sheep fecal material sequenced
with extraordinary depth and diversity [41]; see Table 2 for results of all
samples available and Table S2 for binning details. We tried the heurstices on
rust-mdbg [42] assembly graphs. The rescued circles can recover up to
1/3$\sim$1/4 of HiFi MAGs in the corresponding samples, although with more
indels.

Table 1: CheckM evaluation of rescued circles, MetaBAT2 bins and vamb bins. Vamb requires at least 4096 contigs to run. sample | binning | near-complete | high-quality | near-complete | high-quality
---|---|---|---|---|---
| method | bins | bins | orphaned $>$1Mb | orphaned $>$1Mb
|  |  |  | linear contig | linear contig
sheep-gut-1a | rescue | 45 | 5 | 18 | 17
| metabat2 | 51 | 35 | 8 | 6
| vamb | 31 | 34 | 22 | 15
sheep-gut-1b | rescue | 99 | 15 | 85 | 59
| metabat2 | 128 | 103 | 20 | 21
| vamb | 64 | 71 | 79 | 62
chicken-gut-1 | rescue | 10 | 0 | 8 | 7
| metabat2 | 26 | 18 | 0 | 3
| vamb | 18 | 14 | 4 | 6
env-digester-1 | rescue | 12 | 3 | 3 | 7
| metabat2 | 9 | 18 | 1 | 1
| vamb | 4 | 14 | 4 | 7
human-gut-1 | rescue | 20 | 3 | 13 | 15
| metabat2 | 28 | 32 | 8 | 10
| vamb | 18 | 33 | 11 | 16
human-gut-2 | rescue | 21 | 2 | 9 | 24
| metabat2 | 33 | 49 | 6 | 14
| vamb | 23 | 41 | 8 | 21
human-gut-3 | rescue | 8 | 1 | 4 | 3
| metabat2 | 15 | 4 | 0 | 1
| vamb | 8 | 6 | 6 | 1
human-gut-4 | rescue | 7 | 4 | 2 | 7
| metabat2 | 24 | 30 | 0 | 1
| vamb | 18 | 28 | 2 | 6
human-gut-5 | rescue | 8 | 3 | 1 | 3
| metabat2 | 20 | 13 | 1 | 0
| vamb | N/A |  |  |
human-gut-6 | rescue | 5 | 3 | 2 | 2
| metabat2 | 6 | 10 | 1 | 2
| vamb | 8 | 10 | 2 | 2
human-gut-7 | rescue | 5 | 7 | 8 | 6
| metabat2 | 22 | 17 | 1 | 1
| vamb | 17 | 14 | 5 | 6
huamn-gut-8 | rescue | 3 | 3 | 1 | 4
| metabat2 | 13 | 12 | 0 | 0
| vamb | 8 | 8 | 1 | 5
human-gut-9 | rescue | 4 | 1 | 0 | 3
| metabat2 | 10 | 14 | 0 | 0
| vamb | 22 | 11 | 0 | 0
human-gut-10 | rescue | 12 | 3 | 4 | 6
| metabat2 | 23 | 21 | 0 | 1
| vamb | 19 | 20 | 2 | 3
Table 2: Merging the rescued circles and MetaBAT2 bins. “90-5”: near-complete. CheckM evaluation columns delimited by commas contain bin counts of in the order of near-complete, high-quality and medium-quality. Total MetaBAT2 bins did not double count the long circular contigs. sample | #circular | #rescue | #metabat2 | checkM | checkM | checkM | checkM | #orphan | #present/miss
---|---|---|---|---|---|---|---|---|---
| contig |  | accept/all | circular | rescue | metabat2 | merged | 90-5/ | 90-5
|  |  |  | contig |  |  |  | all | metabat2
sheep-gut-1a | 145 | 61 | 611/668 | 132,10,2 | 45,5,0 | 51,35,37 | 69,32,36 | 1/4 | 27/27
sheep-gut-1b | 279 | 161 | 1269/1431 | 249,22,2 | 99,15,2 | 130,104,111 | 188,102,96 | 4/25 | 37/41
chicken-gut-1 | 73 | 12 | 332/379 | 62,10,1 | 10,0,0 | 26,18,14 | 28,18,14 | 0/0 | 8/8
env-digester-1 | 21 | 41 | 292/316 | 18,2,0 | 12,3,1 | 9,19,11 | 20,19,12 | 0/2 | 1/1
human-gut-1 | 42 | 52 | 470/511 | 37,5,0 | 20,3,2 | 28,32,29 | 37,28,29 | 2/22 | 8/11
human-gut-2 | 46 | 39 | 545/581 | 37,7,1 | 21,2,1 | 34,49,40 | 44,44,37 | 2/11 | 9/11
human-gut-3 | 37 | 12 | 86/128 | 36,0,1 | 8,1,0 | 15,4,8 | 17,4,8 | 0/4 | 5/6
human-gut-4 | 18 | 37 | 234/267 | 16,2,0 | 7,4,3 | 24,30,40 | 17,28,39 | 0/9 | 6/14
human-gut-5 | 10 | 19 | 107/123 | 7,0,0 | 8,3,0 | 20,13,8 | 18,15,7 | 0/4 | 8/10
human-gut-6 | 7 | 15 | 76/93 | 4,1,2 | 5,3,1 | 6,10,11 | 8,10,12 | 0/1 | 2/3
human-gut-7 | 17 | 22 | 166/188 | 16,1,0 | 5,7,1 | 22,17,18 | 21,18,17 | 0/2 | 5/6
huamn-gut-8 | 4 | 15 | 112/124 | 3,0,0 | 3,3,0 | 13,12,15 | 10,11,15 | 0/0 | 3/6
human-gut-9 | 13 | 18 | 107/125 | 7,1,1 | 4,1,3 | 10,14,19 | 7,9,21 | 0/1 | 4/7
human-gut-10 | 23 | 30 | 180/206 | 19,3,0 | 12,3,4 | 23,21,24 | 26,21,26 | 0/6 | 5/9

For sheep-gut-1b (Figure 2), the circle-rescue heuristic first considered 3340
cycles. After deduplication, 161 were reported, 99 of which were near-
complete, 15 were high-quality and 1 was medium-quality. Additionally, 85
$>$1Mb linear contigs that were not used by medium-quality-and-above cycles
were near complete. Only 25 of them could align more than 60% of their length
to a long circular contig or a high-quality-or-above rescued circle (-c
-xasm20). The proportion of candidates failing the checkM quality check was
comparable to traditional binning. For example, MetaBAT2 reported 438 bins
that are at least 1Mb and less than 8Mb, 121 of them were near-complete, 80
were high-quality and 75 were medium-quality. We did not find obvious
difference between the rescued circles that failed checkM, and the rescued
circles in the better quality bracket, in terms of sequence homology to the
long circular contigs. The true positive rate from rescued circles could be
improved by discarding obviously wrong cycles (e.g. repeating similar
subsequences), but it is not useful in practice, as MAGs will always be
quality-controlled by checkM and other reference-based tools.

Figure 2: Clustering of near-complete MAGs of sheep-gut-1b. The middle-left
dense purple square represents Archaea. Mash distance more than 30% are shown
as 30% for a better visual.

When merging the rescued circles with the regular MetaBAT2 binning (Table 2),
most MetaBAT2 bins were accepted (1269/1431), which yielded 89 near-complete
bins and 87 high-quality bins. 25 $>$1Mb linear contigs were accepted due to
orphanage and all of them evaluated to be near-complete. Along with 279 $>$1Mb
circular contigs (249 near-complete and 22 high-quality), this summed to 462
near-complete MAGs. 41 near-complete MetaBAT2 bins were rejected due to
overlapping with rescued circles. We found 37 of them have counterparts in the
collection of the rescued circles with less than 1% mash distance (30 had less
than 0.1% mash distance match), effectively being replaced by them. Out of
these 37 pairs, 11 rescued circles matches were better than the rejected
MetaBAT2 bins in terms of better checkM completeness and/or contamination, 22
had identical checkM evaluations, 2 were slightly worse but comparable, and 2
were worse (rescued circle was wrong). Among the four MetaBAT2 bins without
representation, 3 contained a single long linear contig (wrongly used by
rescued circles), 1 contained two long linear contigs. Overall, we were able
to recover more MAGs than MetaBAT2 alone, and replace many MetaBAT2 bins with
pseudo-circles while having acceptable lose. We note that the lose would be
larger in worse assemblies where the rescued circle is prone to make more
errors, such as human-gut-4 shown in Table 2.

Since both circle-rescuing and the merging allow some contigs to be shared
between bins, there could be duplications. To check this, we calculate
pairwise mash distance between near-complete MAGs from the above. MetaBAT2 bin
with more than 1 contig was concatenated with 31 N-base paddings for
convenience. Mash sketch does hash k-mer with N-base instead of ignoring them,
but this has little influence on the distance estimation here. There were 441
MAGs, therefore 97020 unique pairs. 2/97020 had less than 1% mash distance.
One pair was between a linear long contig and a rescued circle’s path
belonging to the same not fully resolved subgraph (0.7%; both 91% complete and
0.6% contaminated). The other pair was between two circular contigs (0.7%;
both 100% complete and 1.1% contaminated). Overall, the MAG collection had
little redundancy.

The MAG merging procedure described above does not look at existing gene
annotations. It is possible to achieve better binning if we use checkM to
guide binning (e.g. DAS Tools [36]). We showcase a simple checkM-guided bin
merging to demonstrate the gap between merged MAGs and the almost-best-
possible outcome in the near-complete quality bracket (Table 3). We say
“almost” because we do not try to create new bins or swap bin contents, but
just to avoid shadowing valid MetaBAT2 bins with wrong rescued circles. First,
the circle-rescuing heuristics and MetaBAT2 binning were executed as described
above. We then ran checkM on the rescued circles, MetaBAT2 bins and all 1Mb
linear contigs separately. For each quality bracket, we first accept all
rescued circles that qualify, then accept MetaBAT2 bins that qualify and have
less than 100kb of its content used by rescued circles, followed by accepting
$>$1Mb linear contigs that qualify and have not been used by the above two
categories, and finally all $>$1Mb circular contigs that qualify.

Table 3: CheckM-assisted best-possible MAG recovery for each sample. First four numbers of each field give the count of MAGs from $>$1Mb circular contigs, rescued circles, metabat2 bins and $>$1Mb orphaned linear contigs of the corresponding quality bracket, respectively. The fifth number gives the sum. sample | #near-complete | #high-quality | #medium-quality
---|---|---|---
sheep-gut-1a | 132+45+26+6=209 | 10+5+33+8=56 | 2+0+37+19=58
sheep-gut-1b | 249+99+94+19=461 | 22+15+101+30=168 | 2+2+110+49=163
chicken-gut-1 | 62+10+18+0=90 | 10+0+18+7=35 | 1+0+14+13=28
env-digester-1 | 18+12+8+0=38 | 2+3+19+5=29 | 0+1+11+5=17
human-gut-1 | 37+20+21+7=85 | 5+3+31+15=54 | 0+2+29+30=61
human-gut-2 | 37+21+25+5=88 | 7+2+48+19=76 | 1+1+39+34=75
human-gut-3 | 36+8+10+0=54 | 0+1+4+2=7 | 1+0+8+6=15
human-gut-4 | 16+7+18+0=41 | 2+4+30+4=40 | 0+3+39+10=52
human-gut-5 | 7+8+15+0=30 | 0+3+13+3=19 | 0+0+8+17=25
human-gut-6 | 4+5+4+1=14 | 1+3+9+3=16 | 2+1+11+2=16
human-gut-7 | 16+5+19+1=41 | 1+7+14+6=28 | 0+1+18+3=22
huamn-gut-8 | 3+3+11+0=17 | 0+3+11+3=17 | 0+0+15+5=20
human-gut-9 | 7+4+6+0=17 | 1+1+14+1=17 | 1+3+19+4=27
human-gut-10 | 19+12+18+0=49 | 3+3+20+6=32 | 0+4+23+7=34

### Sample representation completeness assessed using k-mer spectrum

We borrow the k-mer spectrum plot presented such as in KAT [43] and mercury
[44]. In single genome assemblies, the main interest is to grasp assembly
redundancies, incompleteness, correction errors and phasing by looking at the
plot. For metagenome assembly, we focus on illustrating the first two aspects.
Metagenome libraries have extremely high counts of low multiplicity k-mers,
resulted from the combination of sequencing errors and reads from low
prevalence species. We defined ratios based on k-mer counts in read occurrence
ranges (rather than exact read occurrence values) to render the spectrum plot
visually easier to interpret (Figure 3). In the modified k-mer spectrum plot,
if haplotypes are not closely related, contain neglectable amount of repeats
and the assembly is perfect, we expect to the see the band representing
$1\times$ assembly multiplicity (saturated orange in such plots throughout
this manuscript; “$1\times$ band” from now on) to dominate the plot within
genome-wide read coverage range, and fall sharply outside of it. There will be
very few k-mers with read multiplicities higher than genome-wide coverage due
to shared sequences that are longer than the k-mers (see below). We also hope
to see few k-mers that exist in the reads but not in contigs, i.e. $0\times$
assembly multiplicity plotted saturated blue at the bottom of each plot, as
they imply unassembled contents. If haplotypes are somewhat related, we expect
to see several peaks formed by the $1\times$ band due to some k-mers being
shared by more than one haplotypes. As a real near-ideal empirical example, we
compared the human HiFi datasets with HRGM [19] (Figure S3), which is a
diverse yet non-redundant collection of 5414 species from human gut
microbiome. This plot marks what the plot would look like for a reasonably
good sample representation. If we can achieve a smaller $0\times$ band than
Figure S3 through de novo assembly, then the assembly offers no less
information than comparing the sample to a comprehensive reference set.

Figure 3: K-mer spectrum plots of all samples. In each subplot: let
$N^{(c)}_{x}$ be the number of k-mers occurring $\geq x$ times in reads and
exactly $c$ times in contigs. Then $N_{x}=\sum_{c}N^{(c)}_{x}$ is the number
of k-mers occurring $\geq x$ times in reads. The height of the blue area
intersecting at $x$ equals $N^{(0)}_{x}/N_{x}$; the height of the orange area
intersecting at $x$ equals $N^{(1)}_{x}/N_{x}$. Similarly we have green for
$c=2$, red for $c=3$, etc. Intuitively speaking, a large orange region
suggests a good assembly with most reads k-mers assembled into a single
contig; a large blue region suggests an incomplete assembly that misses many
high-abundance k-mers in reads; a large region of other colors suggests a
potentially duplicated assembly with k-mers assembled into multiple contigs.
We truncate at x if there are less than 1 million k-mers occurring x or more
times in reads.

Using the k-mer spectrum plots, we found that although hifiasm-meta could
generate more near-complete circular contigs than previous studies, the
circular contigs alone do not provide a complete view of the corresponding
libraries (Figure 3, left column). Unresolved subgraphs in the primary
assembly accounted for quite a few dominant haplotypes (Figure S4). All
samples did not receive good representation from just the circular contigs
except for sheep-gut-1a, but this could be greatly improved if the merged MAGs
are used instead (Figure 3, right column), except for human-gut-9 and env-
digester-1. In four libraries, merged MAG were close to the near-ideal
situation demonstrated above: chicken-gut-1, human-gut-3, sheep-gut-1a and
sheep-gut-1b. We extracted k-mers (n=31) of the $0\times$ band and aligned
them to MAGs in question with bwa aln [45] to see if they can be found by
allowing a few mismatches or indels. This had very insignificant impact to the
plot, however. We also examined the $2\times$ or higher bands in the high
read-multiplicity range by dumping all kmers from human-gut-10 that were from
$2\times$ to $15\times$ bands, and had read multiplicities higher than 800.
They were likely from ubiquitous genes (e.g. tRNAs) or horizontally
transferable sequences (see Methods).

### Cross validation with full length 16S rRNA compositional estimation

Complement to the the k-mer spectrum-based completeness evaluation is the 16S
rRNA-based methods. Most HiFi reads are a few folds longer than the full
length rRNA genes. This provides composition estimation for free. Predicting
the rRNA genes is well-studied [46, 47]. The performance was reasonable in
HiFi reads (e.g. for barrnap, around 13% 16S genes were marked as partial).
16S-based taxonomy annotation has been similarly extensively explored [48, 46,
49], but it suffers from reference scope bias. For example, the human-gut-9
library has 1.8M reads. 24445 16S sequences were identified (1.4%), 21954 of
them (90.0%) were assigned with confident genus level annotation. In contrast,
the env-digester-1 library has 1.0M reads. Similarly, 1.3% of them contain 16S
genes, but only 21.2% which would have confident genus level annotation. This
was not limited to 16S genes. The SRA minhash-based read taxonomy analysis
only identified 41.78% reads as of cellular organism origin. Non-human gut
samples were in the middleground: sheep-gut-1b has 1.2% reads containing 16S
genes, with 47.3% of them confidently annotated to genus level.

Based on these observations, although reference-based methods are widely used
in practice and could resolve to species level [12, 16], we use greedy
incremental clustering to define OTUs with boundary set to mismatch sequence
identity 99% [12] (see Methods). No assembly could recover all abundant OTUs,
but those evaluated to be better in k-mer spectrum approach missed less
(Figure 4; Figure S5 provides plots of all samples and other OTU boundaries).
Some MAGs had more than one OTU assignment due to possible erroneous OTU
clustering, 16S diversity within genome, or collapsing of very similar
haplotypes during assembly. Missing high-prevalence OTUs are usually partially
assembled, i.e. linear contigs or subgraphs that failed to form near-complete
MAGs (see Methods).

Figure 4: OTU recovery of merged bins, showing sheep-gut-1b and env-
digester-1, which evaluated poorly and good in k-mer spectrum plot,
respectively.

## Discussion

### Drawbacks of short read metagenome assemblies and the justification of the
cost of using accurate long reads.

With binning, the near-complete MAG yield from HiFi assemblies, usually
1$\sim$2 such MAGs per Gb per sample, is similar to or less than the yields of
past short read studies (Figure S1) and currently at a substantially higher
cost. It might seem that the only advantage of HiFi assembly is the contig
continuity, which is important but not so much for gene-focused approaches.
However, SR assemblies often miss the major species at MAG level, and what
have been missed is not easy to guess. While this could be overcame if using a
comprehensive catalog as reference and approach the sequencing library from
alignment-based or pangenome point of view, actual studies tend to analyze or
compare MAGs without accounting for possible missing species. To illustrate
the issue, we use Almeida et al.’s collection of human gut MAGs [18]. The
authors assembled 13133 human gut metagenomic datasets from 75 studies with a
unified pipeline (metaSPAdes and MetaBAT2, quality controlled by checkM),
ruling out impact of any pipeline difference. We collected k-mer spectrum
plots for 3490 assemblies (Table S3) that had at least 10 MAGs. We define a
good assembly to have more than half of its 0x band below the 50% watermark, a
nice assembly if using the 30% watermark, and a failed assembly otherwise.
Only 0.46% assemblies were nice. 2.1% were good. The completeness of the final
catalog might be at the mercy of assemblies of unrelated samples luckily
recovering each others’ most abundant species. In contrast, except for env-
digester-1 and human-gut-9, all other HiFi assemblies presented by the merged
MAGs were nice. Moreover, combined with the representation-completeness
evaluations from two different aspects, we hope that it will be possible to
transform the ideal sequencing depth from simply common practices (e.g.
$\leq$10 Gb for human gut samples) to be educated guesses. This would be
especially important when faced with novel environmental samples.

Other than incomplete sample representation, due to short contig length and
binning errors, SR MAGs are known to be less complete than they appear to be.
Meziti et al. noticed that on average, 23% population core genes and 50%
variable genes were missing from near 95% complete MAGs [50]. We found most SR
MAGs from short read studies to be smaller than their counterparts found in
the HiFi MAGs (Figure 5); mash distance $<$5%). 12 pairs between Almeida et
al. MAGs and HiFi contigs had checkM completeness $>$99% and contamination
$<$1%, yet the SR MAGs were smaller than the HiFi contigs for more than 100Kb.
When compared to HRGM, this trend held. However, some HiFi circular contigs
appeared smaller. HRGM is a non-redundant species-level collection. It is not
clear whether we actually had smaller genomes (or if HRGM favored larger
genomes when deduplicating), or the smaller size was caused by assembly
errors. Anyway, this is outside the scope of this paper and we will
investigate in future studies.

Figure 5: Most near-complete circular contigs and rescued cycles of human gut
samples are longer than assemblies in human gut short read MAG collection
(Almeida et al.) or reference catalog (HRGM, the UHGG extended). We pair our
MAGs with theirs if the pair’s mash distance is less than 5%. Reference MAGs
are allowed to be used more than one time. If there is a tie between reference
choices, we arbitrarily pick one. We have no requirements for MAG length when
pairing. The Jaccard index (which is the basis of mash distance) already
penalizes set size difference, i.e. comparing a set with its strict superset
gives imperfect score.

Haplotypes represented by the rescued circles do not appear to be particularly
harder to recover than those assembled into circular contigs in the general
sample pool (Figure 2). At species level (mash distance $<$5%), most of our
near-complete circular contigs (180/182; 163/180 if comparing to Almeida et
al.) or near-complete rescued circles (93/93; 90/93 if Almeida et al.) from
human samples had matches in HRGM. At strain level (mash distance $<$1%),
51/182 (28%) near-complete circular contigs and 16/93 (17%) near-complete
rescued circles were found in Almeida et al.’s (difference not significant
between the two: fisher’s exact, p-value=0.07).

### Undocumented or unseen MAGs.

For our near-complete MAGs in human samples (8 individuals plus 2 4-pooled
libraries), we noticed some species or genera were absent from public
repositories or MAG collections from large scale studies, using a cutoff at
mash dist 0.05 (roughly $\geq 95\%$ ANI). 47 (11%) unseen if compared to only
near-complete MAGs from Almeida et al, 17 (4%) if compared to all reported
MAGs from them. 217 (51%) were not in refseq. For $\geq$ 1Mb circular contigs,
that would be 33 (15%), 11 (5%) and 122 (56%), respectively. These ratios are
consistent with previous study. Jin et al. reported [51] 475 high quality
(checkM completeness $\geq 80\%$, contamination $\leq 5\%$, quality score
$\geq 60$) from pacbio-illumina hybrid assembly and binning. 4 MAGs were
circular. The MAGs were compared to the UHGG dataset [52] with ANI cutoff at
95% and reported 24 MAGs as novel. In our comparison settings, the number of
unseen MAGs are 43 (9%), 21 (4%) and 193 (41%), respectively. It is hard to
verify whether these observations are due to true compositional diversity
between individual samples [22] or representation-incomplete assembly of past
projects. We did not find significant coverage difference between our MAGs
that were unseen in previous works and the rest.

MAG catalogs other than the human gut microbiome’s are far from complete. We
compared the chicken-gut-1 and the sheep samples to their relevant MAG
catalogs, ICRGGC (chicken-gut-1) and PRJNA657473(ruminants), respectively.
ICRGGC [53] contained 12339 MAGs derived from 799 samples. PRJNA657473 [54]
contained 10371 MAGs derived from 370 samples (7 sample species, 10
gastrointestinal tract regions). Near-complete circular contigs and rescued
circles from chicken-gut-1 found 55/62 (contig) and 8/10 (rescued) species-
level matches, 7/62 (contig) and 0/10 (rescued) strain-level matches. sheep-
gut-1b found no matches for 249 contigs and 99 rescued circles in both
category, probably due to insufficient relevant reference.

## Conclusions

In this study, we implemented a topology-based binning heuristic on the contig
graph for hifiasm-meta and described two approaches, namely k-mer spectrum and
species-evel OTUs based on full length 16S rRNAs, to evaluate an metagenome
assembly’s completeness in terms of sample representation. We further showed
that de novo HiFi assemblies plus binning would have the potential to be both
genome-complete and representation-complete, bringing MAGs closer to their
original goal which is to delegate their biosample with minimum bias. Note
that inferring whether a sequencing library has sampled the microbiome of
interest exhaustively is hard and out of the scope of this work. There are
other drawbacks, too. For example, the binning heuristic and the traditional
binners still assume some reference bias due to being dependent on checkM as
the final genome quality validation.

Nonetheless, we anticipate that high quality metagenome assemblies and further
method improvements could transform previously inaccessible approaches, such
as analyzing horizontal gene transfers, de novo variant calling in unusual
samples and direct comparison between microbial communities.

## Methods

### Assembly, evaluation and simulation

We generated HiFi assemblies using hifiasm-meta r63 with default settings. See
Table S5 for versions of tools used. Short read assemblies were downloaded
from their releases. We used MetaBAT2 [32] with contig coverage Contig
coverage for MetaBAT2 and vamb binning was estimated with minimap2 alignment
and MetaBAT2’s jgi module. Coverage: we ran minimap2 [55] with “minimap2 -ak19
-w10 -I10G -g5k -r2k --lj-min-ratio 0.5 -A2 -B5 -O5,56 -E4,1 -z400,50 --sam-
hit-only contigs.fa reads.fa”. BAM file handling used SAMtools [56].

Coverage was estimated by “jgi_summa_rsize_bam_contig_depths \--outputDepth
depth.txt input.bam”. Binning: we ran MetaBAT2 using “metabat2 --seed 1 -i
contigs.fa -a depth.txt”, and vamb using “vamb -p 48 --outdir ./ --fasta
contigs.fa --jgi jgi_depth --minfasta 500000”. MetaBAT2’s random seed has
little influence. We separate circular contigs of 1Mb or longer into a
separate MAG if it is binned together with other contigs.

We used checkM module “lineage_wf” to evaluate MAG quality. Its outputs were
formatted by “checkm qa -o 2” before parsing. We did not try DAS tools’s
evaluation in this work, but it should give consistent but more generous
results.

We did HiFi read simulation using PBSIM2 [57]: “pbsim2 --depth INT”
“\--sample-fastq sample.fp”. We generated “sample.fq”, the empirical error
profile, by randomly sampling 100k reads from sheep-gut-1a. We used seqtk
(https://github.com/lh3/seqtk): “seqtk sample sheep-gut-1a.fq 100000 >
out.fq”.

We forked yak (https://github.com/lh3/yak) for k-mer spectrum plot. The fork
used k=31 rather than the k=27 default of KAT. This did not have significant
impact on the plots; k=21 generated overall similar outputs as well.

We ran rust-mdBG with “-k 21 -l 14 --density 0.003 -p asm”, then
“magic_simplify_meta asm” to generate the final assembly graph and the
sequences per developer’s recommendation. A freestanding implementation of the
circle-finding heuristics was used, and we used mash distance (cutoff: 90%) to
compare between the reported circular paths to merged MAGs of hifiasm-meta’s.
In sludge, humanO3 and sheepB, this reported 22, 22 and 120 rescued circles,
respectively. We did not try to do metaBAT2 and bin merging because checkM was
sensitive to indels.

### Examining high read- and assembly-multiplicity kmers

We use human-gut-10 as an example. We identified kmers with at least 2x and up
to 15x assembly multiplicity, and at least 800x read multiplicity, i.e. the
right-most part of k-mer spectrum plot above 1x band. Their location on the
contigs were collected. We merged overlapping intervals and dumped the
sequences. There were 5469 unique sequences (max length 18.6kb, N50 2.0kb). We
randomly select 20 from these (“seqtk subseq in.fa 20”) and did BLAST (blastn
web cgi, defaults) against nr/nt. All queries had full length BLAST hits with
low sequence divergence and frequently overlap with genes encoding DNA-related
enzymes, transposase and tRNA or rRNA (Table S6).

### 16S rRNA methods

We identified and annotated 16S rRNA genes from HiFi libraries with the following steps. First, HiFi reads that could aligned to SILVA reference were extracted, with base qualities stripped: “seqtk subseq hifi.fq <(minimap2 SILVA.fa hifi.fq | cut -f1 | uniq) | seqtk seq -A > SSUreads.fa” . We ran barrnap to identify rRNA genes: “barrnap --kingdom bac --outseq rRNA.fa SSUreads.fa”. INFERNAL cmsearch might identify a few more rRNAs than barrnap. We believe this would not have major influence on the conclusion based on previous observations. We then annotated rRNA genes with RDP classifier: “java -Xmx16g -jar RDPTools/classifier.jar classify -o RDP.tsv -h RDP.hier rRNA.fa” . We accept annotations of 16S rRNAs with genus scores of at least 0.9.

To define OTUs from HiFi reads, we first selected 16S genes not marked with
“partial” from the barrnap. We used greedy incremental clustering: we
initialize an empty collection $S$ to collect seed sequences. For each 16S
gene $q$, if it could align to any sequence $s$ in $S$ with alignment block
longer than 1000bp and at least 99% mismatch identity, it is assigned the same
OTU label as $s$. (If multiple seeds are available, the one with highest
identity will be chosen. If a tie, the seed is arbitrarily chosen from the
bests.) Otherwise, $q$ is added to $S$ and assigned a new OTU label. Alignment
is done with minimap2’s python binding, mappy, with “preset=map-hifi”.
Alignment block length is given by “mappy.Alignment.blen”. Mismatch identity
is calculated as “mappy.Alignment.mlen/mappy.Alignment.blen”. Assigning OTU
label for an unseen 16S copy is done similarly. If a sequence can not align to
any seed sequences, its OTU label is undefined.

There are two ways to assign OTU labels to MAGs: 1) collected reads belonging
to contigs of a MAG and their OTU labels, or 2)identify 16S copies from
contigs then assign labels. We did both and found them to be mostly
consistent. We only considered near-complete MAGs.

When evaluating MAGs as in Figure 4 and Figure S5, we drop OTUs with less than
10 16S copies to rule out artifacts from sequencing errors and to ignore
species with very low coverage. This was done after the OTU label assignment.
A MAG could have more than one OTU assignment. It was difficult to distinguish
wrong cases (i.e. suboptimal clustering result) from true cases, i.e. a genome
having multiple distinct 16S copies, therefore we simply accept all OTU labels
of a MAG. For example, if the read set yields 3 OTUs ($a$,$b$ and $c$), the
assembly has a single MAG from which we identify three 16S copies that are
labeled $a$, $a$, $b$. Then in the plot, both $a$ and $b$ would be colored as
“seen in MAG”.

## Supplementary Information

## Declarations

### Ethics approval and consent to participate

Not applicable.

### Consent for publication

Not applicable.

### Availability of data and materials

Hifiasm-meta is open source at https://github.com/xfengnefx/hifiasm-meta/. The
yak fork used is at https://github.com/xfengnefx/yam. Scripts and assemblies
are at https://github.com/xfengnefx/snpt_mtgncomp and
ftp://ftp.dfci.harvard.edu/pub/hli/hifiasm-meta/hifiasm-binning/asm/,
respectively. Hifi libraries [29, 41, 58, 59] are all publicly available and
can be found under their accession IDs (table S4). For convenience, in Table
S3 we supply urls to the short read libraries we used when evaluating short
read MAGs from Almeida et al (the accession IDs were originally available in
their supplementary tables).

### Competing interests

The authors declare that they have no competing interests.

### Funding

This work was funded by National Human Genome Research Institute (NHGRI)
R01HG010040 and U01HG010961, the Chan Zuckerberg Initiative and the Sloan
foundation.

### Authors’ contributions

HL and XF conceived the projected, carried out analysis and wrote the
manuscript. All authors read and approved the final manuscript.

### Acknowledgements

Not applicable.

## References

  * [1] Tully, B.J., Graham, E.D., Heidelberg, J.F.: The reconstruction of 2,631 draft metagenome-assembled genomes from the global oceans. Scientific data 5(1), 1–8 (2018)
  * [2] Howe, A.C., Jansson, J.K., Malfatti, S.A., Tringe, S.G., Tiedje, J.M., Brown, C.T.: Tackling soil diversity with the assembly of large, complex metagenomes. Proceedings of the National Academy of Sciences 111(13), 4904–4909 (2014)
  * [3] Kroeger, M.E., Delmont, T.O., Eren, A.M., Meyer, K.M., Guo, J., Khan, K., Rodrigues, J.L., Bohannan, B.J., Tringe, S.G., Borges, C.D., et al.: New biological insights into how deforestation in amazonia affects soil microbial communities using metagenomics and metagenome-assembled genomes. Frontiers in microbiology 9, 1635 (2018)
  * [4] Parks, D.H., Imelfort, M., Skennerton, C.T., Hugenholtz, P., Tyson, G.W.: CheckM: assessing the quality of microbial genomes recovered from isolates, single cells, and metagenomes. Genome Res. 25(7), 1043–1055 (2015)
  * [5] Bowers, R.M., Kyrpides, N.C., Stepanauskas, R., Harmon-Smith, M., Doud, D., Reddy, T., Schulz, F., Jarett, J., Rivers, A.R., Eloe-Fadrosh, E.A., et al.: Minimum information about a single amplified genome (misag) and a metagenome-assembled genome (mimag) of bacteria and archaea. Nature biotechnology 35(8), 725–731 (2017)
  * [6] Albertsen, M., Hugenholtz, P., Skarshewski, A., Nielsen, K.L., Tyson, G.W., Nielsen, P.H.: Genome sequences of rare, uncultured bacteria obtained by differential coverage binning of multiple metagenomes. Nature biotechnology 31(6), 533–538 (2013)
  * [7] Vicedomini, R., Quince, C., Darling, A.E., Chikhi, R.: Strainberry: automated strain separation in low-complexity metagenomes using long reads. Nature Communications 12(1), 1–14 (2021)
  * [8] Nayfach, S., Shi, Z.J., Seshadri, R., Pollard, K.S., Kyrpides, N.C.: New insights from uncultivated genomes of the global human gut microbiome. Nature 568(7753), 505–510 (2019)
  * [9] Luo, C., Tsementzi, D., Kyrpides, N.C., Konstantinidis, K.T.: Individual genome assembly from complex community short-read metagenomic datasets. The ISME journal 6(4), 898–901 (2012)
  * [10] Stackebrandt, E., Goebel, B.M.: Taxonomic note: a place for dna-dna reassociation and 16s rrna sequence analysis in the present species definition in bacteriology. International journal of systematic and evolutionary microbiology 44(4), 846–849 (1994)
  * [11] Brumfield, K.D., Huq, A., Colwell, R.R., Olds, J.L., Leddy, M.B.: Microbial resolution of whole genome shotgun and 16s amplicon metagenomic sequencing using publicly available neon data. PLoS One 15(2), 0228899 (2020)
  * [12] Edgar, R.C.: Updating the 97% identity threshold for 16s ribosomal rna otus. Bioinformatics 34(14), 2371–2375 (2018)
  * [13] Yuan, C., Lei, J., Cole, J., Sun, Y.: Reconstructing 16s rrna genes in metagenomic data. Bioinformatics 31(12), 35–43 (2015)
  * [14] Poretsky, R., Rodriguez-R, L.M., Luo, C., Tsementzi, D., Konstantinidis, K.T.: Strengths and limitations of 16s rrna gene amplicon sequencing in revealing temporal microbial community dynamics. PloS one 9(4), 93827 (2014)
  * [15] Kai, S., Matsuo, Y., Nakagawa, S., Kryukov, K., Matsukawa, S., Tanaka, H., Iwai, T., Imanishi, T., Hirota, K.: Rapid bacterial identification by direct pcr amplification of 16s rrna genes using the minion™ nanopore sequencer. FEBS open bio 9(3), 548–557 (2019)
  * [16] Johnson, J.S., Spakowicz, D.J., Hong, B.-Y., Petersen, L.M., Demkowicz, P., Chen, L., Leopold, S.R., Hanson, B.M., Agresta, H.O., Gerstein, M., et al.: Evaluation of 16s rrna gene sequencing for species and strain-level microbiome analysis. Nature communications 10(1), 1–11 (2019)
  * [17] Wang, W., Hu, H., Zijlstra, R.T., Zheng, J., Gänzle, M.G.: Metagenomic reconstructions of gut microbial metabolism in weanling pigs. Microbiome 7(1), 1–11 (2019)
  * [18] Almeida, A., Mitchell, A.L., Boland, M., Forster, S.C., Gloor, G.B., Tarkowska, A., Lawley, T.D., Finn, R.D.: A new genomic blueprint of the human gut microbiota. Nature 568(7753), 499–504 (2019)
  * [19] Kim, C.Y., Lee, M., Yang, S., Kim, K., Yong, D., Kim, H.R., Lee, I.: Human reference gut microbiome catalog including newly assembled genomes from under-represented asian metagenomes. Genome Medicine 13(1), 1–20 (2021)
  * [20] Pasolli, E., Asnicar, F., Manara, S., Zolfo, M., Karcher, N., Armanini, F., Beghini, F., Manghi, P., Tett, A., Ghensi, P., et al.: Extensive unexplored human microbiome diversity revealed by over 150,000 genomes from metagenomes spanning age, geography, and lifestyle. Cell 176(3), 649–662 (2019)
  * [21] Chen, C., Zhou, Y., Fu, H., Xiong, X., Fang, S., Jiang, H., Wu, J., Yang, H., Gao, J., Huang, L.: Expanded catalog of microbial genes and metagenome-assembled genomes from the pig gut microbiome. Nature communications 12(1), 1–13 (2021)
  * [22] Visconti, A., Le Roy, C.I., Rosa, F., Rossi, N., Martin, T.C., Mohney, R.P., Li, W., de Rinaldis, E., Bell, J.T., Venter, J.C., et al.: Interplay between the human gut microbiome and host metabolism. Nature communications 10(1), 1–10 (2019)
  * [23] Liu, P., Hu, S., He, Z., Feng, C., Dong, G., An, S., Liu, R., Xu, F., Chen, Y., Ying, X.: Towards strain-level complexity: Sequencing depth required for comprehensive single-nucleotide polymorphism analysis of the human gut microbiome. Frontiers in Microbiology 13 (2022)
  * [24] Hillmann, B., Al-Ghalith, G.A., Shields-Cutler, R.R., Zhu, Q., Gohl, D.M., Beckman, K.B., Knight, R., Knights, D.: Evaluating the information content of shallow shotgun metagenomics. Msystems 3(6), 00069–18 (2018)
  * [25] Institute, D.J.G. https://www.ncbi.nlm.nih.gov/bioproject/?term=Wetland+microbial+communities+from+Twitchell+Island (2011)
  * [26] Kanehisa, M., Sato, Y., Morishima, K.: Blastkoala and ghostkoala: Kegg tools for functional characterization of genome and metagenome sequences. Journal of molecular biology 428(4), 726–731 (2016)
  * [27] Frioux, C., Singh, D., Korcsmaros, T., Hildebrand, F.: From bag-of-genes to bag-of-genomes: metabolic modelling of communities in the era of metagenome-assembled genomes. Computational and Structural Biotechnology Journal 18, 1722–1734 (2020)
  * [28] Laudadio, I., Fulci, V., Palone, F., Stronati, L., Cucchiara, S., Carissimi, C.: Quantitative assessment of shotgun metagenomics and 16s rdna amplicon sequencing in the study of human gut microbiome. OMICS: A Journal of Integrative Biology 22(4), 248–254 (2018)
  * [29] Kolmogorov, M., Bickhart, D.M., Behsaz, B., Gurevich, A., Rayko, M., Shin, S.B., Kuhn, K., Yuan, J., Polevikov, E., Smith, T.P.L., Pevzner, P.A.: metaFlye: scalable long-read metagenome assembly using repeat graphs. Nat. Methods 17(11), 1103–1110 (2020)
  * [30] Nurk, S., Walenz, B.P., Rhie, A., Vollger, M.R., Logsdon, G.A., Grothe, R., Miga, K.H., Eichler, E.E., Phillippy, A.M., Koren, S.: HiCanu: accurate assembly of segmental duplications, satellites, and allelic variants from high-fidelity long reads. Genome Res. 30(9), 1291–1305 (2020)
  * [31] Feng, X., Cheng, H., Portik, D., Li, H.: Metagenome assembly of high-fidelity long reads with hifiasm-meta. Nature Methods, 1–4 (2022)
  * [32] Kang, D.D., Li, F., Kirton, E., Thomas, A., Egan, R., An, H., Wang, Z.: MetaBAT 2: an adaptive binning algorithm for robust and efficient genome reconstruction from metagenome assemblies. PeerJ 7, 7359 (2019)
  * [33] Nissen, J.N., Johansen, J., Allesøe, R.L., Sønderby, C.K., Armenteros, J.J.A., Grønbech, C.H., Jensen, L.J., Nielsen, H.B., Petersen, T.N., Winther, O., et al.: Improved metagenome binning and assembly using deep variational autoencoders. Nature biotechnology 39(5), 555–560 (2021)
  * [34] Mallawaarachchi, V., Wickramarachchi, A., Lin, Y.: Graphbin: refined binning of metagenomic contigs using assembly graphs. Bioinformatics 36(11), 3307–3313 (2020)
  * [35] Wu, Y.-W., Tang, Y.-H., Tringe, S.G., Simmons, B.A., Singer, S.W.: Maxbin: an automated binning method to recover individual genomes from metagenomes using an expectation-maximization algorithm. Microbiome 2(1), 1–18 (2014)
  * [36] Sieber, C.M., Probst, A.J., Sharrar, A., Thomas, B.C., Hess, M., Tringe, S.G., Banfield, J.F.: Recovery of genomes from metagenomes via a dereplication, aggregation and scoring strategy. Nature microbiology 3(7), 836–843 (2018)
  * [37] Wick, R.R., Schultz, M.B., Zobel, J., Holt, K.E.: Bandage: interactive visualization of de novo genome assemblies. Bioinformatics 31(20), 3350–3352 (2015)
  * [38] Ondov, B.D., Treangen, T.J., Melsted, P., Mallonee, A.B., Bergman, N.H., Koren, S., Phillippy, A.M.: Mash: fast genome and metagenome distance estimation using minhash. Genome biology 17(1), 1–14 (2016)
  * [39] Jain, C., Rodriguez-R, L.M., Phillippy, A.M., Konstantinidis, K.T., Aluru, S.: High throughput ani analysis of 90k prokaryotic genomes reveals clear species boundaries. Nature communications 9(1), 1–8 (2018)
  * [40] Parks, D.H., Chuvochina, M., Chaumeil, P.-A., Rinke, C., Mussig, A.J., Hugenholtz, P.: A complete domain-to-species taxonomy for bacteria and archaea. Nature biotechnology 38(9), 1079–1086 (2020)
  * [41] Bickhart, D.M., Kolmogorov, M., Tseng, E., Portik, D.M., Korobeynikov, A., Tolstoganov, I., Uritskiy, G., Liachko, I., Sullivan, S.T., Shin, S.B., et al.: Generating lineage-resolved, complete metagenome-assembled genomes from complex microbial communities. Nature biotechnology, 1–9 (2022)
  * [42] Ekim, B., Berger, B., Chikhi, R.: Minimizer-space de bruijn graphs: Whole-genome assembly of long reads in minutes on a personal computer. Cell systems 12(10), 958–968 (2021)
  * [43] Mapleson, D., Garcia Accinelli, G., Kettleborough, G., Wright, J., Clavijo, B.J.: Kat: a k-mer analysis toolkit to quality control ngs datasets and genome assemblies. Bioinformatics 33(4), 574–576 (2017)
  * [44] Rhie, A., Walenz, B.P., Koren, S., Phillippy, A.M.: Merqury: reference-free quality, completeness, and phasing assessment for genome assemblies. Genome biology 21(1), 1–27 (2020)
  * [45] Li, H., Durbin, R.: Fast and accurate short read alignment with burrows–wheeler transform. bioinformatics 25(14), 1754–1760 (2009)
  * [46] Seemann, T., Booth, T.: Barrnap: BAsic Rapid Ribosomal RNA Predictor (2018)
  * [47] Nawrocki, E.P., Eddy, S.R.: Infernal 1.1: 100-fold faster rna homology searches. Bioinformatics 29(22), 2933–2935 (2013)
  * [48] Wang, Q., Garrity, G.M., Tiedje, J.M., Cole, J.R.: Naive bayesian classifier for rapid assignment of rrna sequences into the new bacterial taxonomy. Applied and environmental microbiology 73(16), 5261–5267 (2007)
  * [49] Edgar, R.C.: Accuracy of taxonomy prediction for 16s rrna and fungal its sequences. PeerJ 6, 4652 (2018)
  * [50] Meziti, A., Rodriguez-R, L.M., Hatt, J.K., Peña-Gonzalez, A., Levy, K., Konstantinidis, K.T.: The reliability of metagenome-assembled genomes (mags) in representing natural populations: Insights from comparing mags against isolate genomes derived from the same fecal sample. Applied and environmental microbiology 87(6), 02593–20 (2021)
  * [51] Jin, H., You, L., Zhao, F., Li, S., Ma, T., Kwok, L.-Y., Xu, H., Sun, Z.: Hybrid, ultra-deep metagenomic sequencing enables genomic and functional characterization of low-abundance species in the human gut microbiome. Gut microbes 14(1), 2021790 (2022)
  * [52] Almeida, A., Nayfach, S., Boland, M., Strozzi, F., Beracochea, M., Shi, Z.J., Pollard, K.S., Sakharova, E., Parks, D.H., Hugenholtz, P., et al.: A unified catalog of 204,938 reference genomes from the human gut microbiome. Nature biotechnology 39(1), 105–114 (2021)
  * [53] Feng, Y., Wang, Y., Zhu, B., Gao, G.F., Guo, Y., Hu, Y.: Metagenome-assembled genomes and gene catalog from the chicken gut microbiome aid in deciphering antibiotic resistomes. Communications biology 4(1), 1–9 (2021)
  * [54] Xie, F., Jin, W., Si, H., Yuan, Y., Tao, Y., Liu, J., Wang, X., Yang, C., Li, Q., Yan, X., et al.: An integrated gene catalog and over 10,000 metagenome-assembled genomes from the gastrointestinal microbiome of ruminants. Microbiome 9(1), 1–20 (2021)
  * [55] Li, H.: Minimap2: pairwise alignment for nucleotide sequences. Bioinformatics 34(18), 3094–3100 (2018)
  * [56] Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N., Marth, G., Abecasis, G., Durbin, R.: The sequence alignment/map format and samtools. Bioinformatics 25(16), 2078–2079 (2009)
  * [57] Ono, Y., Asai, K., Hamada, M.: Pbsim2: a simulator for long-read sequencers with a novel generative model of quality scores. Bioinformatics 37(5), 589–595 (2021)
  * [58] Sereika, M., Kirkegaard, R.H., Karst, S.M., Michaelsen, T.Y., Sørensen, E.A., Wollenberg, R.D., Albertsen, M.: Oxford nanopore r10. 4 long-read sequencing enables the generation of near-finished bacterial genomes from pure cultures and metagenomes without short-read or reference polishing. Nature methods 19(7), 823–826 (2022)
  * [59] Kim, C.Y., Ma, J., Lee, I.: Hifi metagenomic sequencing enables assembly of accurate and complete genomes from human gut microbiota. bioRxiv (2022)

## Figures

## Supplementary Figures

Supplementary Figure 1: Almeida et al. unified assembly and MAG quality
evaluations for 13,133 human gut metagenomic datasets from 75 short-read-based
studies. Most libraries are smaller than 10Gbp (99.2%), and per Gbp MAG
recovery visually negatively correlated with library size. For an “average”
library, we expect 1-2 near-complete MAGs. HiFi assemblies provide closed
contigs and much lower intra-MAG heterozygosity (0.06% near-complete circular
contigs have ¿25% checkM heterogenous score, versus 23.8% for SR) at similar
per Gbp yields.  Supplementary Figure 2: A few runs from short read metagenome
assembly studies achieved reasonable sample representation. However, they were
fragmented and larger library size did not imply a better outcome.
Supplementary Figure 3: K-mer spectrum plots using HRGM assemblies as the
MAGs, and HiFi reads as the library. Intend to show how good would a
relatively complete library represent a new sample.  Supplementary Figure 4:
Coverage of sheep-gut-1b merged MAGs by category.

Supplementary Figure 5: OTU recovery in samples, with OTU boundary set at 99%,
95% and 90%. The plots for human-gut-1 and human-gut-2 were omitted, because
the each of these libraries are a pool of four different samples.

## Tables

## Additional Files

### Table S1. Supplementary tables and data release used a different sample
naming convention. This file provides the name mapping.

### Table S2. Binning information and evaluation of bins.

A tab-delimited table.

### Table S3. The 3490 samples from Almeida et al. used in this manuscript.

List of files (sra ftq).

### Table S4. Sample and HiFi library information.

### Table S5. List of tools their versions used in this manuscript.

### Table S6. BLAST result summary of 20 high read- and contig- multiplicity
k-mers intervals.