Model building of protein complexes from intermediate-resolution cryo-EM maps with deep learning-guided automatic assembly.

Applying main-chain probability maps to other methods

One important strategy in EMBuild is the use of the main-chain probability maps. To investigate whether such main-chain probability map can help model building of other approaches, we also applied the main-chain probability maps to phenix.dock_in_map and gmfit on the test set of 47 single-particle EM maps. Specifically, for each input protein chain, only main-chain atoms (N, C, and Cα) are kept during fitting into the main-chain probability map, and the rest atoms are ignored until the assembling is finished. The comparisons for model building on the raw EM density map and main-chain probability map are shown in Supplementary Fig. 3 . The detailed results are listed in Supplementary Data 2 . It can be seen from the figure that the performance of phenix.dock_in_map was significantly improved on the main-chain probability maps. The average TM-score, RMSD, CC_box, and CC_mask values achieved by phenix.dock_in_map on the main-chain probability map are 0.807, 4.03 Å, 0.6265, and 0.5686, compared to 0.746, 4.52 Å, 0.5975, and 0.5323 on the original EM density map. The improvement of phenix.dock_in_map on the main-chain probability maps demonstrated the importance of the main-chain probability maps in model building. However, gmfit did not benefit from main-chain probability maps. This is understandable because gmfit converts the input map into a Gaussian mixture model in the fitting procedure, which may miss detailed structural information in the main-chain probability maps.

Methods

Data collection

We have collected a non-redundant dataset of cryo-EM maps from the EMDB. The primary maps, which are normally sharpened maps, are used in this study.All the single-particle EM entries at 4–8 Å resolution that have associated PDB models were downloaded from the EMDB 17 and PDB 18 . Specifically, the PDB structure for each entry was directly downloaded from the PDB at https://www.rcsb.org/ , and its corresponding density volume within 4.0 Å of the PDB structure was segmented out from the whole map. The EM map and its corresponding PDB structure that have any of the following features were removed: (i) without side-chain atoms, (ii) including unknown residues (UNK), (iii) including missing chain or non-protein chain, (iv) having nonorthogonal map axis, (v) corresponding to multiple PDB or EMDB entries, and (vi) having severe misfits between the PDB model and EM map. In order to ensure the feature (vi), we calculated the cross-correlation between the deposited map and the map simulated from the PDB model at the same resolution using the UCSF chimera 30 . Any map and its associated PDB model that have a cross-correlation of less than 0.65 were excluded. Afterwards, all the remaining maps were manually checked. The remaining cases were retained as the initial dataset. To remove redundancy, the initial dataset of cases was clustered using greedy algorithm. Two models are considered to be similar if any chain in the first model has >30% sequence identity with any chain in the second model. The one with the largest number of similar cases is chosen as the representative of the corresponding cluster, and then the rest cases in the cluster are removed. This procedure is repeated until all the cases are clustered. The resulted non-redundant training set consists of the representatives of each cluster. A total of 262 pairs of EM maps and associated PDB structures were retained. In order to train a deep learning model for predicting main-chain probability map, 209 maps were randomly selected as the final training set (Supplementary Data 1 ). To build an independent test set with a sufficient number of valuable cases, all the cases in the initial dataset that have >30% sequence identity with any case in the training set are excluded. Then, any pair of EM map and its associated PDB structure that meet the following criteria are removed: (i) having only one chain and (ii) having >30% gap in the structure according to the gene sequence. For saving time in evaluations, we have also removed those cases with more than 10 chains, though our method can work with any number of chains. The remaining cases are clustered using a sequence identity cutoff of 70% by the similar greedy algorithm described above. The final test set contains 47 single-particle cryo-EM maps with resolutions ranging from 4.0 to 8.0 Å (Supplementary Data 2 ). In the modeling and evaluating processes, short chains that have less than 20 residues are ignored. To evaluate the performance of EMBuild on the maps by subtomogram averaging of cryo-ET data, we further constructed another independent test set of subtomogram averaging maps. All the subtomogram averaging maps with resolutions within 10.0 Å that have associated PDB models are filtered and clustered using the same way as that for the test set of single-particle maps, except for allowing more than 10 chains. The final test set contains 16 subtomogram averaging maps with resolutions ranging from 3.7 to 9.3 Å (Supplementary Data 3 ).

Main-chain probability prediction by deep learning

We used a deep learning model to predict the main-chain probability map from the input raw EM density map through a nested U-net (UNet++) 51 . The network architecture consists of three encoder sub-blocks and three decoder sub-blocks with dense skip connections, where 3D convolution layers with a kernel size of 3 × 3 × 3 are applied (Fig. 1a ). The 3D maxpooling layer with stride of 2 is adopted for down-sampling, and the trilinear interpolation layer with zoom factor of 2 is adopted for up-sampling. During the prediction, the raw EM map is rescaled to have a grid interval of 1.0 Å by applying cubic interpolation and cut into overlapping chunks of size 40 × 40 × 40 with slide strides of 10 voxels, which are input into the deep learning model. Then, the predicted main-chain probability chunks with the same size and grid interval are re-assembled into the final main-chain probability map by averaging overlapping parts. A set of pairs of experimental EM density and main-chain probability maps were used to train the deep learning model. For the EM maps, the grid interval is rescaled to 1.0 Å by applying cubic interpolation. The density values are clipped to be equal or greater than 0.0 and are normalized to the range 0.0–1.0 by the 99.999-percentile density value of the map. For a given EM density map, its main-chain probability map is generated from its associated PDB structure with 1.0 Å grid interval, where the main-chain probability p for a grid point x is calculated as follows 1 documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$p({{{{{bf{x}}}}}})=,max {{e}^{-lambda {Vert {{{{{bf{x}}}}}}-{{{{{bf{a}}}}}}Vert }^{2}},forall {{{{{bf{a}}}}}}in A}$$end{document} p ( x ) = max { e − λ ∥ x − a ∥ 2 , ∀ a ∈ A } where A stands for the set of position vectors of all main-chain atoms (N, C, and Cα). The value of λ is defined as 2 documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$lambda ={(pi /(2.4+0.8R))}^{2}$$end{document} λ = ( π / ( 2.4 + 0.8 R ) ) 2 where R is the map resolution 43 . For training, the EM density maps and their corresponding main-chain probability maps are cut into pairs of overlapping boxes of size 60 × 60 × 60 with slide strides of 30 voxels, where non-positive boxes are excluded for effective training. During the training, 20% of the maps are randomly selected from the training set as the validation set. The training data are augmented through random 90° rotations and randomly cropping a 40 × 40 × 40 chunk from the input 60 × 60 × 60 box. The network is implemented in Python with Pytorch1.8.1 + cuda11.1. For each model, the network is trained for at most 300 epochs with 160 boxes employed in one batch. The Adam optimizer is adopted to minimize the loss of the prediction, where the total loss is a sum of two different loss functions. One is the smooth L1 loss, which calculates the numerical difference in the probability values between the predicted chunk X and target chunk Y as follows. 3 documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$${{{{{rm{Smoot}}}}}}{{{{{rm{h}}}}}}{{{{{rm{L}}}}}}1{{{{{rm{Loss}}}}}}(X,Y)=mathop{sum }limits_{i=1}^{N}mathop{sum }limits_{j=1}^{N}mathop{sum }limits_{k=1}^{N}frac{{L}_{i,j,k}}{{N}^{3}}$$end{document} Smoot h L 1 Loss ( X , Y ) = ∑ i = 1 N ∑ j = 1 N ∑ k = 1 N L i , j , k N 3 where N is the chunk size, and documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$${L}_{i,j,k}$$end{document} L i , j , k is the Smooth L1 distance between X and Y at position ( i , j , k ). The other is the structural similarity (SSIM) loss which compares the contrast and structure similarity between a predicted chunk X with its target chunk Y according to the following formula 4 documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$${{{{{rm{SSIMLoss}}}}}}left(X,Yright)=1-frac{2{sigma }_{{XY}}+varepsilon }{{sigma }_{X}^{2}+{sigma }_{Y}^{2}+varepsilon }$$end{document} SSIMLoss X , Y = 1 − 2 σ X Y + ε σ X 2 + σ Y 2 + ε where σ X and σ Y are the standard deviations for the predicted chunk X and target chunk Y , σ XY is the covariance between X and Y , and ε is set to be a small constant ( documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$varepsilon ={10}^{-6}$$end{document} ε = 10 − 6 in this study) to prevent dividing by zero. The learning rate is initially set to 10 −3 and will be reduced to 1/2 of its current value if the average loss on the training set does not decrease for 4 epochs. The training process continues until the learning rate reaches a minimum value of 10 −5 . The network model with the least validation loss is selected. As seen from the learning curves, the training and validating losses converge well for our deep learning model (Supplementary Fig. 6 ). Fitting protein chains into the main-chain probability map The structure models of individual chains can be predicted from their sequences by a protein structure prediction program. In this study, AlphaFold2 was selected for this purpose, given its excellent performance in protein structure prediction 52 . To mimic real situations, the corresponding PDB structure and those newer structures are excluded from the templates by setting the “max template data” to the day before the released date of the corresponding PDB structure 52 , and full-length gene sequences are used as the input. For consistency, the predicted model of each chain is cropped to have identical starting and ending residues with the sequence in the PDB structure. It should be noted that during the implementation of EMBuild there is no structural gap and/or uninterpreted region in our predicted chain models, while they may possibly present in the PDB structure. For each chain model predicted by AlphaFold2, we use SWORD to assign its structural domains 64 . SWORD will generate multiple assignments for one chain. The assignment with the most domains with no less than 30 residues and a κ value of no less than 3.0 is chosen as the final domain assignment. For computational efficiency, instead of directly using the predicted main-chain probability map, a reduced representation of main-chain probabilities is adopted by EMBuild. Namely, a mean shift algorithm is applied to generate representative points of main-chain probabilities. Specifically, starting from the positions of positive grid points on the main-chain probability map, the seed points documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$${{{{{{bf{z}}}}}}}_{i}^{t}(i=1,...,{N}^{{prime} }{{{{{rm{;}}}}}},t=0,1,...)$$end{document} z i t ( i = 1 , . . . , N ′ ; t = 0 , 1 , . . . ) is iteratively shifted to local maxima of probabilities as 5 documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$${{{{{{bf{z}}}}}}}_{i}^{t+1}=frac{{sum }_{n=1}^{N}K({{{{{{bf{z}}}}}}}_{i}^{t}-{{{{{{bf{x}}}}}}}_{n})p({{{{{{bf{x}}}}}}}_{n}){{{{{{bf{x}}}}}}}_{n}}{{sum }_{{n}^{{prime} }=1}^{N}K({{{{{{bf{z}}}}}}}_{i}^{t}-{{{{{{bf{x}}}}}}}_{{n}^{{prime} }})p({{{{{{bf{x}}}}}}}_{{n}^{{prime} }})}$$end{document} z i t + 1 = ∑ n = 1 N K ( z i t − x n ) p ( x n ) x n ∑ n ′ = 1 N K ( z i t − x n ′ ) p ( x n ′ ) where documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$${{{{{{bf{x}}}}}}}_{n}(n=1,...,N)$$end{document} x n ( n = 1 , . . . , N ) are the position vectors of grid points, documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$K({{{{{{bf{z}}}}}}}_{i}^{t}-{{{{{{bf{x}}}}}}}_{n})$$end{document} K ( z i t − x n ) is the Gaussian kernel function, and documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$p({{{{{{bf{x}}}}}}}_{n})$$end{document} p ( x n ) is the main-chain probability of grid point x n . The Gaussian kernel function is described as documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$K({{{{{{bf{z}}}}}}}_{i}^{t}-{{{{{{bf{x}}}}}}}_{n})={e}^{-lambda {Vert {{{{{{bf{z}}}}}}}_{i}^{t}-{{{{{{bf{x}}}}}}}_{n}Vert }^{2}}$$end{document} K ( z i t − x n ) = e − λ ∥ z i t − x n ∥ 2 . The main-chain probability of a shifted seed point documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$P({{{{{{bf{z}}}}}}}_{i}^{t})$$end{document} P ( z i t ) is computed as documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$P({{{{{{bf{z}}}}}}}_{i}^{t})=frac{1}{N}{sum }_{n=1}^{N}K({{{{{{bf{z}}}}}}}_{i}^{t}-{{{{{{bf{x}}}}}}}_{n})p({{{{{{bf{x}}}}}}}_{n})$$end{document} P ( z i t ) = 1 N ∑ n = 1 N K ( z i t − x n ) p ( x n ) . After the mean shift procedure is converged, the seed points that are closer than a threshold distance are clustered and the one with the highest probability value is chosen as the representative of each cluster. The resulted points documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$${{{{{{bf{z}}}}}}}_{i}(i=1,...,L)in Z$$end{document} z i ( i = 1 , . . . , L ) ∈ Z are referred to as the main-chain points. We adopt a fast Fourier transform (FFT)-based matching strategy to globally fitting the protein model for each chain to the main-chain probability map (main-chain points) 65 . To perform an exhaustive FFT-based search, both the protein model and main-chain points are first mapped onto a three-dimensional (3D) grid of shape M × M × M with 1.5 Å grid interval. The main-chain probabilities are assigned to the grid of the main-chain points (say grid A ) and the grid of the main-chain atoms in the protein model (say grid B ), according to Eq. ( 1 ). With the above main-chain probability mapped on grids, the match score S for a superimposition between the protein model and main-chain probability map can be generally expressed by the following formula 6 documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$S(i,j,k)=-theta mathop{sum }limits_{l=1}^{M}mathop{sum }limits_{m=1}^{M}mathop{sum }limits_{n=1}^{M}({A}_{l,m,n}times {B}_{l+i,m+j,n+k})$$end{document} S ( i , j , k ) = − θ ∑ l = 1 M ∑ m = 1 M ∑ n = 1 M ( A l , m , n × B l + i , m + j , n + k ) where i , j , and k are the numbers of grid points by which the protein model is shifted with respect to the main-chain points in three translational dimensions, and θ is the resolution-dependent factor defined as documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$theta ={(lambda /pi )}^{1.5}$$end{document} θ = ( λ / π ) 1.5 . The match scores for all the M 3 translations can be computed through one round of FFT-based calculation. The rotational search is conducted by exploring a large set of rotation angles. That is, for each rotation of the protein structure, an FFT-based translational search is carried out. For EMBuild, an angle interval of 15° is used to evenly discretize the Euler space, which results in a total of 4392 evenly distributed orientations. The fitting results of the exhaustive search are further optimized through a SIMPLEX method. The match score s ′ of main-chain atom documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$${{{{{{bf{y}}}}}}}_{q}(q=1,...,Q)$$end{document} y q ( q = 1 , . . . , Q ) in protein pose Y can be calculated as 7 documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$s^{prime} ({{{{{{bf{y}}}}}}}_{q};Z)=-max {theta P({{{{{bf{z}}}}}}){e}^{-lambda Vert {{{{{{bf{y}}}}}}}_{q}-{{{{{{bf{z}}}}}}}{Vert}^{2}},forall {{{{{bf{z}}}}}}in Z}$$end{document} s ′ ( y q ; Z ) = − max { θ P ( z ) e − λ ∥ y q − z ∥ 2 , ∀ z ∈ Z } where documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$${{{{{bf{z}}}}}}in Z$$end{document} z ∈ Z is the position of the main-chain point and documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$P({{{{{bf{z}}}}}})$$end{document} P ( z ) is its main-chain probability. The match score for pose Y is the summation of scores of included atoms as documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$${S}^{{prime} }({Y;Z})={sum }_{q}^{Q}{s}^{{prime} }({{{{{{bf{y}}}}}}}_{q}{;Z})$$end{document} S ′ ( Y ; Z ) = ∑ q Q s ′ ( y q ; Z ) . Finally, the fitting results are ranked by the match scores and those top-scored poses are retained. After the protein model of a chain is rigidly fitted to the main-chain probability map, we employed a semi-flexible domain refinement strategy to further improve the fitness between the protein models the map, as illustrated in Fig. 2a . For each protein model, n short structure domains are assigned by SWORD 64 . A simple graph is built based on the domain assignment, where two domains with connecting residues are connected by an undirected edge. The domain refinement is applied on each of the M top-scored poses from rigid fitting. Starting from a selected domain as the seed domain, the positions of all domains are optimized one after another. Specifically, starting from the seed domain, a SIMPLEX optimization is conducted to find a locally best match of the current domain, and then the optimization is carried out to the neighboring domain of current domain. This procedure is repeated until all the domains are locally optimized. The order of the refinement is determined by breadth-first search (BFS) on the domain graph. Taking each domain as the seed domain, the domain refinement will output n different models. The results of domain refinement for one protein chain are M × n refined models plus M rigidly fitted models, which are ranked by their match scores. Assembling chains into protein complex After fitting and refinement of individual chains, the last step of EMBuild is to assemble individual chains into a reliable complex structure. To prevent clash between different chains, we defined a clash score c of one atom a of chain A with respect to all atoms b of another chain B as follows 8 documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$c({{{{{bf{a}}}}}};B)=,max {{e}^{-lambda {({{max }}{Vert {{{{{bf{a}}}}}}-{{{{{bf{b}}}}}}Vert -{d}_{{{{{{rm{clash}}}}}}},0.0})}^{2}},forall {{{{{bf{b}}}}}}in B}$$end{document} c ( a ; B ) = max { e − λ ( max { ∥ a − b ∥ − d clash , 0.0 } ) 2 , ∀ b ∈ B } where documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$${d}_{{{{{{rm{clash}}}}}}}$$end{document} d clash stands for the cutoff distance within which the clash score is set to 1.0. The clash score C of chain A with respect to chain B is the average clash score of all n A atoms with documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$C({A;},B)=frac{1}{{n}_{A}}mathop{sum}limits_{{{{{{bf{a}}}}}}{in}{{{{{boldsymbol{in }}}}}}A}c({{{{{bf{a}}}}}};,B)$$end{document} C ( A ; B ) = 1 n A ∑ a ∈ ∈ A c ( a ; B ) . According to the match scores of individual chains and the clash scores between different chains, the problem of assembling chains into a complex will become a Maximum Clique Problem. First, an undirected graph is built by EMBuild, of which the vertices are the match scores of individual chains, and the edges connecting two different chains are the corresponding clash scores. The different fitting poses of one chain are not connected. Then, the edges that have a clash score exceeding a given threshold C thr are removed from the graph. After the graph is built, the Bron-Kerbosch algorithm is used to find the best combination of chains from the graph that has the highest total match score. However, though only a single cycle of Bron–Kerbosch algorithm, some chains may not be assembled into the protein complex. Therefore, we adopt an iterative strategy in EMBuild, as illustrated in Fig. 2b . Namely, after each round of the Bron-Kerbosch algorithm, the clash scores of main-chain points are calculated according to the current assembled complex structure D , which are used to update the probability documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$Pleft({{{{{{bf{z}}}}}}}_{i}right)(i=1,ldots ,L)$$end{document} P z i ( i = 1 , … , L ) of main-chain points, according to the following formula 9 documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$${P}^{{prime} }left({{{{{{bf{z}}}}}}}_{i}{{{{{rm{;}}}}}},Dright)=Pleft({{{{{{bf{z}}}}}}}_{i}right)times (1.0-c({{{{{{bf{z}}}}}}}_{i}{{{{{rm{;}}}}}},D))$$end{document} P ′ z i ; D = P z i × ( 1.0 − c ( z i ; D ) ) By updating the probabilities of main-chain points, the regions with fitted structures are removed from further assembling. The remaining chains are iteratively assembled to the complex structure by fitting and refining in accordance with the updated main-chain points. Finally, the resulted complex model is refined in the EM density map using phenix.real_space_refine 47 .

Evaluation metrics

Three types of metrics are used to evaluate the quality of the protein complex model built by EMBuild. The first type of metrics is one that measure the closeness between the built model and the PDB structure. In this respect, we adopt the TM-scores and RMSDs between the built complex model and the corresponding PDB structure calculated by MM-align 54 , 55 . The second type of metrics is those fit-to-map metrics of the built models 56 , which measure the consistency between the built model and the EM density map. In this regard, we report CC_box and CC_mask values calculated by phenix.map_model_cc and the map-model FSC05 calculated by phenix.mtriage 57 . It is noted that phenix.mtriage will fail to give a valid FSC05 value if the built model does not conform to the map, where the maximum value of the map-model FSC is only around or even below 0.5 over the entire resolution range. Besides the above metrics, we also reported coordinates-only metrics including the Ramachandran scores and MolProbity score calculated by MolProbity 59 .

Quality assessment of built models

To assess the quality of the model built by EMBuild, we propose a metric of the main-chain match score, which measures the fitness between the built model and the main-chain probability map. The match scores for individual main-chain atoms are calculated according to Eq. ( 7 ). Starting from the match scores of main-chain atoms, it is easy to calculate the average score of three main-chain atoms as the match score for each residue. By analogy, we further calculate the match score of each continuous secondary structure fragment as the average match score of its containing residues. The secondary structure is assigned by STRIDE 66 . The final match score of a fragment is defined as a combination of its initial score and the match score of its domain as 10 documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$${S}_{{{{{{rm{fragment}}}}}}}^{{prime} }=0.7times {S}_{{{{{{rm{fragment}}}}}}}+0.3times {S}_{{{{{{rm{domain}}}}}}}$$end{document} S fragment ′ = 0.7 × S fragment + 0.3 × S domain where documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$${S}_{{{{{{rm{domain}}}}}}}$$end{document} S domain is the average match score of all residues in the domain that contains the fragment. Similar to the fragment match score, the final match score of a residue is defined as a combination of its initial score and the match score of its fragment as 11 documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$${S}_{{{{{{rm{residue}}}}}}}^{{prime} }=0.7times {S}_{{{{{{rm{residue}}}}}}}+0.3times {S}_{{{{{{rm{fragment}}}}}}}$$end{document} S residue ′ = 0.7 × S residue + 0.3 × S fragment The resulted residue scores are further smoothed along the chain with a sliding window of weights 1:2:4:8:16:8:4:2:1 centered at each residue. To measure the quality of modeled fragments with respect to the PDB structure, we define an alignment score, which is described as 12 documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$${{{{{rm{alignment}}}}}}-{{{{{rm{score}}}}}}=frac{1}{L}mathop{sum }limits_{i=1}^{L}frac{1}{1+{d}_{i}^{2}/{d}_{0}^{2}}$$end{document} alignment − score = 1 L ∑ i = 1 L 1 1 + d i 2 / d 0 2 where L is the length of the fragment, d i is the distance of the i th pair of the aligned residues between the built model and the fragment in the PDB structure, and d 0 is a scale factor. Although our alignment score takes a similar expression to TM-score 55 , two significant differences should be noticed. One is that no superposition is applied before calculating the alignment score. The other is that d 0 is set to a fixed distance of 3.0 Å.

Comparison with related methods

EMBuild is compared with phenix.dock_in_map, DEMO-EM, and gmfit on the test sets of cryo-EM maps. For each test case, the protein models of individual chains are predicted from sequences by AlphaFold2, which are used as the input for different methods. To be general, the symmetry information of test cases is ignored during the evaluation. phenix.dock_in_map uses both the secondary structure matching and convolution-based shape searches to find a part of a map that is similar to a protein model, which can be used to place any number of copies of any number of unique molecules 42 . gmfit is a program for fitting subunits into a density map using GMM (Gaussian Mixture Model) 31 , 32 . To convert the density maps and protein chains into a GMM used by gmfit, the number of Gaussian functions for a density map is set to 20 multiplied by the number of chains in its associated protein complex. The number of Gaussian functions for a query protein chain is set to 20. For gmfit, the number of randomly generated initial configurations is set to 100000, the number of configurations for search is set to 20000, and the number of configurations for refinement is set to 4000. DEMO-EM is a hierarchical method to assemble multi-domain protein structures from cryo-EM density maps 53 . We have tested DEMO-EM on the test set of single-particle maps in two ways. One is inputting protein sequences, where structure predictions are carried out by DEMO-EM itself. The other is providing structures predicted by AlphaFold2 to DEMO-EM. It should be noted that 22 cases with a total sequence length >2000 are not accepted by the DEMO-EM server. Thus, the rest 25 cases were submitted to the DEMO-EM online server, of which EMD-22216 failed using AlphaFold2 structures as the input. For every model building method, the resulted complex model is refined in the EM density map using phenix.real_space_refine 47 . Reporting summary Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Comparison with related methods

EMBuild is compared with phenix.dock_in_map, DEMO-EM, and gmfit on the test sets of cryo-EM maps. For each test case, the protein models of individual chains are predicted from sequences by AlphaFold2, which are used as the input for different methods. To be general, the symmetry information of test cases is ignored during the evaluation. phenix.dock_in_map uses both the secondary structure matching and convolution-based shape searches to find a part of a map that is similar to a protein model, which can be used to place any number of copies of any number of unique molecules 42 . gmfit is a program for fitting subunits into a density map using GMM (Gaussian Mixture Model) 31 , 32 . To convert the density maps and protein chains into a GMM used by gmfit, the number of Gaussian functions for a density map is set to 20 multiplied by the number of chains in its associated protein complex. The number of Gaussian functions for a query protein chain is set to 20. For gmfit, the number of randomly generated initial configurations is set to 100000, the number of configurations for search is set to 20000, and the number of configurations for refinement is set to 4000. DEMO-EM is a hierarchical method to assemble multi-domain protein structures from cryo-EM density maps 53 . We have tested DEMO-EM on the test set of single-particle maps in two ways. One is inputting protein sequences, where structure predictions are carried out by DEMO-EM itself. The other is providing structures predicted by AlphaFold2 to DEMO-EM. It should be noted that 22 cases with a total sequence length >2000 are not accepted by the DEMO-EM server. Thus, the rest 25 cases were submitted to the DEMO-EM online server, of which EMD-22216 failed using AlphaFold2 structures as the input. For every model building method, the resulted complex model is refined in the EM density map using phenix.real_space_refine 47 .

Supplementary information Supplementary Information Peer Review File Description of Additional Supplementary Files Supplementary Data 1 Supplementary Data 2 Supplementary Data 3 Supplementary Data 4 Supplementary Data 5 Supplementary Data 6 Supplementary Data 7 Supplementary Data 8 Reporting Summary Soure Data