/* command line argument parsing and document of US-align main program */ #include "MMalign.h" #include "SOIalign.h" #include "flexalign.h" using namespace std; void print_version() { cout << "\n" " ********************************************************************\n" " * US-align (Version 20260819) *\n" " * Universal Structure Alignment of Proteins and Nucleic Acids *\n" " * Reference: C Zhang, L Freddolino, Y Zhang. (2026) Nat Protoc *\n" " * C Zhang, M Shine, AM Pyle, Y Zhang. (2022) Nat Methods*\n" " * C Zhang, AM Pyle (2022) iScience. *\n" " * Please email comments and suggestions to zhang@zhanggroup.org *\n" " ********************************************************************" << endl; } void print_extra_help() { cout << "Additional options:\n" " -v Print the version of US-align\n" "\n" " -a TM-score normalized by the average length of two structures\n" " T or F, (default F). -a does not change the final alignment.\n" "\n" " -fast Fast but slightly inaccurate alignment\n" "\n" " -dir Perform all-against-all alignment among the list of PDB\n" " chains listed by 'chain_list' under 'chain_folder'.\n" " $ USalign -dir chain_folder/ chain_list\n" "\n" //"-dirpair Perform batch alignment for each pair of chains listed by\n" //" 'chain_pair_list' under 'chain_folder'. Each line consist of\n" //" two chains, separated by tab or space.\n" //" $ USalign -dirpair chain_folder/ chain_pair_list\n" //"\n" " -dir1 Use chain2 to search a list of PDB chains listed by 'chain1_list'\n" " under 'chain1_folder'.\n" " $ USalign -dir1 chain1_folder/ chain1_list chain2\n" "\n" " -dir2 Use chain1 to search a list of PDB chains listed by 'chain2_list'\n" " under 'chain2_folder'\n" " $ USalign chain1 -dir2 chain2_folder/ chain2_list\n" "\n" " -suffix (Only when -dir1 and/or -dir2 are set, default is empty)\n" " add file name suffix to files listed by chain1_list or chain2_list\n" "\n" " -atom 4-character atom name used to represent a residue.\n" " Default is \" C3'\" for RNA/DNA and \" CA \" for proteins\n" " (note the spaces before and after CA).\n" "\n" " -split Whether to split PDB file into multiple chains\n" " 0: treat the whole structure as one single chain\n" " (default for -TMscore 2)\n" " 1: treat each MODEL as a separate chain\n" " 2: (default for other cases) treat each chain as a separate chain\n" "\n" " -outfmt Output format\n" " 0: (default) full output\n" " 1: fasta format compact output\n" " 2: tabular format very compact output\n" " -1: full output, but without version or citation information\n" "\n" " -TMcut -1: (default) do not consider TMcut\n" " Values in [0.5,1): Do not proceed with TM-align for this\n" " structure pair if TM-score is unlikely to reach TMcut.\n" " TMcut is normalized as set by -a option:\n" " -2: normalized by longer structure length\n" " -1: normalized by shorter structure length\n" " 0: (default, same as F) normalized by second structure\n" " 1: same as T, normalized by average structure length\n" "\n" " -mirror Whether to align the mirror image of input structure\n" " 0: (default) do not align mirrored structure\n" " 1: align mirror of Structure_1 to origin Structure_2,\n" " which usually requires the '-het 1' option:\n" " $ USalign 4glu.pdb 3p9w.pdb -mirror 1 -het 1\n" "\n" " -het Whether to align residues marked as 'HETATM' in addition to 'ATOM '\n" " 0: (default) only align 'ATOM ' residues\n" " 1: align both 'ATOM ' and 'HETATM' residues\n" " 2: align both 'ATOM ' and MSE residues\n" "\n" " -full Whether to show full pairwise alignment of individual chains for\n" " -mm 2 or 4. T or F, (default F)\n" //"\n" //" -closeK Number of closest atoms used for sequence order independent\n" //" initial alignment. default: 5\n" "\n" " -hinge Maximum number of hinge allowed in flexible alignment.\n" " Only functional with '-mm 7'. default: 9\n" "\n" " -afp Enable AFP-enhancement mechanism. Only functional with '-mm 7'.\n" "\n" " -TMpass Early stopping threshold for AFP-enhancement mechanism.\n" " Only functional with '-mm 7 -afp'. default: 0.85\n" "\n" " -se Do not perform superposition. Useful for extracting alignment from\n" " superposed structure pairs\n" "\n" " -infmt1 Input format for structure_1\n" " -infmt2 Input format for structure_2\n" " -1: (default) automatically detect PDB or PDBx/mmCIF format\n" " 0: PDB format\n" " 1: SPICKER format\n" //" 2: xyz format\n" " 3: PDBx/mmCIF format\n" "\n" "-chainmap (only useful for -mm 1) use the final chain mapping 'chainmap.txt'\n" " specified by user. 'chainmap.txt' is a tab-seperated text with two\n" " columns, one for each complex\n" "\n" "-chain1 Chains to parse in structure_1\n" "-chain2 Chains to parse in structure_2. Use _ for a chain without chain ID.\n" " Multiple chains can be separated by commas, e.g.,\n" " USalign -chain1 C,D,E,F 5jdo.pdb -chain2 A,B,C,D 3wtg.pdb -ter 0\n" "\n" "-model1 Models to parse in structure_1\n" "-model2 Models to parse in structure_2.\n" " Multiple models can be separated by commas, e.g.,\n" " USalign -model1 1,2 1a03.pdb -model2 3,4 1a0n.pdb -ter 0\n" "\n" "Advanced usage 1 (generate an image for a pair of superposed structures):\n" " USalign 1cpc.pdb 1mba.pdb -o sup\n" " pymol -c -d @sup_all_atm.pml -g sup_all_atm.png\n" "\n" "Advanced usage 2 (a quick search of query.pdb against I-TASSER PDB library):\n" " wget https://zhanggroup.org/library/PDB.tar.bz2\n" " tar -xjvf PDB.tar.bz2\n" " USalign query.pdb -dir2 PDB/ PDB/list -suffix .pdb -outfmt 2 -fast\n" << endl; } void print_help(bool h_opt = false) { print_version(); cout << "\n" "Usage: USalign PDB1.pdb PDB2.pdb [Options]\n" "\n" "Options:\n" " -mol Type of molecule(s) to align.\n" " auto: (default) align both protein and nucleic acids.\n" " prot: only align proteins in a structure.\n" " RNA : only align RNA and DNA in a structure.\n" "\n" " -mm Multimeric alignment option:\n" " 0: (default) alignment of two monomeric structures\n" " 1: alignment of two multi-chain oligomeric structures\n" " 2: alignment of individual chains to an oligomeric structure\n" " $ USalign -dir1 monomers/ list oligomer.pdb -ter 0 -mm 2\n" " 3: alignment of circularly permuted structure\n" " 4: MSTA, i.e., alignment of multiple monomeric chains into a\n" " consensus alignment\n" " $ USalign -dir chains/ list -suffix .pdb -mm 4\n" " 5: fully non-sequential (fNS) alignment\n" " 6: semi-non-sequential (sNS) alignment\n" " 7: flexible structure alignment\n" " To use -mm 1 or -mm 2, '-ter' option must be 0 or 1.\n" "\n" " -ter Number of chains to align.\n" " 0: align all chains from all models (recommended for aligning\n" " biological assemblies, i.e. biounits)\n" " 1: align all chains of the first model (recommended for aligning\n" " asymmetric units, default for -mm 1,2 and -TMscore 2,6,7)\n" " 2: (default for other cases) only align the first chain\n" " 3: only align the first chain, or the first segment of the\n" " first chain as marked by the 'TER' string in PDB file\n" "\n" " -TMscore Whether to perform TM-score superposition without structure-based\n" " alignment. The same as -byresi.\n" " 0: (default) sequence independent structure alignment\n" " 1: superpose two structures by assuming that a pair of residues\n" " with the same residue index are equivalent between the two\n" " structures\n" " 2: superpose two complex structures, assuming that a pair of\n" " residues with the same residue index and the same chain ID\n" " are equivalent between the two structures\n" //" 3: (similar to TMscore '-c' option; used with -ter 0 or 1)\n" //" align by residue index and order of chain\n" //" 4: sequence dependent alignment: perform Needleman-Wunsch\n" //" global sequence alignment, followed by TM-score superposition\n" " 5: sequence dependent alignment: perform glocal sequence\n" " alignment followed by TM-score superposition.\n" " -byresi 5 is the same as -seq\n" " 6: superpose two complex structures by first deriving optimal\n" " chain mapping, followed by TM-score superposition for residues\n" " with the same residue ID\n" " 7: sequence dependent alignment of two complex structures:\n" " perform global sequence alignment of each chain pair, derive\n" " optimal chain mapping, and then superpose two complex\n" " structures by TM-score\n" "\n" " -I Use the final alignment specified by FASTA file 'align.txt'\n" "\n" " -i Use alignment specified by 'align.txt' as an initial alignment\n" "\n" " -m Output rotation matrix for superposition, e.g., '-m matrix.txt'\n" " prints the matrix to 'matrix.txt'; '-m -' prints to stdout.\n" "\n" " -d TM-score scaled by an assigned d0, e.g., '-d 3.5' reports MaxSub\n" " score, where d0 is 3.5 Angstrom. -d does not change final alignment.\n" "\n" " -u TM-score normalized by an assigned length. It should be >= length\n" " of protein to avoid TM-score >1. -u does not change final alignment.\n" "\n" " -o Output superposed structure1 to sup.* for PyMOL viewing.\n" " $ USalign structure1.pdb structure2.pdb -o sup\n" " $ pymol -d @sup.pml # C-alpha trace aligned region\n" " $ pymol -d @sup_all.pml # C-alpha trace whole chain\n" " $ pymol -d @sup_atm.pml # full-atom aligned region\n" " $ pymol -d @sup_all_atm.pml # full-atom whole chain\n" " $ pymol -d @sup_all_atm_lig.pml # full-atom with all molecules\n" "\n" " -rasmol Output superposed structure1 to sup.* for RasMol viewing.\n" " $ USalign structure1.pdb structure2.pdb -rasmol sup\n" " $ rasmol -script sup # C-alpha trace aligned region\n" " $ rasmol -script sup_all # C-alpha trace whole chain\n" " $ rasmol -script sup_atm # full-atom aligned region\n" " $ rasmol -script sup_all_atm # full-atom whole chain\n" " $ rasmol -script sup_all_atm_lig # full-atom with all molecules\n" "\n" "-chimerax Output superposed structure1 to sup.* for ChimeraX viewing.\n" " $ USalign structure1.pdb structure2.pdb -chimerax sup\n" " $ chimerax --script sup.cxc # C-alpha trace aligned region\n" " $ chimerax --script sup_all.cxc # C-alpha trace whole chain\n" " $ chimerax --script sup_atm.cxc # full-atom aligned region\n" " $ chimerax --script sup_all_atm.cxc # full-atom whole chain\n" " $ chimerax --script sup_all_atm_lig.cxc # full-atom with all molecules\n" "\n" " -do Output distance of aligned residue pairs\n" "\n" //" -h Print the full help message, including additional options\n" //"\n" "Example usages ('gunzip' program is needed to read .gz compressed files):\n" " USalign 101m.cif.gz 1mba.pdb # pairwise monomeric protein alignment\n" " USalign 1qf6.cif 5yyn.pdb.gz -mol RNA # pairwise monomeric RNA alignment\n" " USalign model.pdb native.pdb -TMscore 1 # calculate TM-score between two conformations of a monomer\n" " USalign 4v4a.cif 4v49.cif -mm 1 -ter 1 # oligomeric alignment for asymmetic units\n" " USalign 3ksc.pdb1 4lej.pdb1 -mm 1 -ter 0 # oligomeric alignment for biological units\n" " USalign 1ajk.pdb.gz 2ayh.pdb.gz -mm 3 # circular permutation alignment\n" << endl; // if (h_opt) print_extra_help(); exit(EXIT_SUCCESS); } /* TMalign, RNAalign, CPalign, TMscore */ int TMalign(string &xname, string &yname, const string &fname_super, const string &fname_lign, const string &fname_matrix, vector &sequence, const double Lnorm_ass, const double d0_scale, const bool m_opt, const int i_opt, const int o_opt, const int a_opt, const bool u_opt, const bool d_opt, const double TMcut, const int infmt1_opt, const int infmt2_opt, const int ter_opt, const int split_opt, const int outfmt_opt, const bool fast_opt, const int cp_opt, const int mirror_opt, const int het_opt, const string &atom_opt, const bool autojustify, const string &mol_opt, const string &dir_opt, const string &dirpair_opt, const string &dir1_opt, const string &dir2_opt, const vector &chain2parse1, const vector &chain2parse2, const vector &model2parse1, const vector &model2parse2, const int byresi_opt, const vector &chain1_list, const vector &chain2_list, const bool se_opt, const bool do_opt) { /* declare previously global variables */ vector >PDB_lines1; // text of chain1 vector >PDB_lines2; // text of chain2 vector mol_vec1; // molecule type of chain1, RNA if >0 vector mol_vec2; // molecule type of chain2, RNA if >0 vector chainID_list1; // list of chainID1 vector chainID_list2; // list of chainID2 int i, j; // file index int chain_i, chain_j; // chain index int r; // residue index int xlen, ylen; // chain length int xchainnum, ychainnum; // number of chains in a PDB file char *seqx, *seqy; // for the protein sequence char *secx, *secy; // for the secondary structure double **xa, **ya; // for input vectors xa[0...xlen-1][0..2] and // ya[0...ylen-1][0..2], in general, // ya is regarded as native structure // --> superpose xa onto ya vector resi_vec1; // residue index for chain1 vector resi_vec2; // residue index for chain2 int read_resi = byresi_opt; // whether to read residue index if (byresi_opt == 0 && o_opt) read_resi = 2; /* loop over file names */ for (i = 0; i < chain1_list.size(); i++) { /* parse chain 1 */ xname = chain1_list[i]; xchainnum = get_PDB_lines(xname, PDB_lines1, chainID_list1, mol_vec1, ter_opt, infmt1_opt, atom_opt, autojustify, split_opt, het_opt, chain2parse1, model2parse1); if (!xchainnum) { cerr << "Warning! Cannot parse file: " << xname << ". Chain number 0." << endl; continue; } for (chain_i = 0; chain_i < xchainnum; chain_i++) { xlen = PDB_lines1[chain_i].size(); if (mol_opt == "RNA") mol_vec1[chain_i] = 1; else if (mol_opt == "protein") mol_vec1[chain_i] = -1; if (!xlen) { cerr << "Warning! Cannot parse file: " << xname << ". Chain length 0." << endl; continue; } else if (xlen < 3) { cerr << "Sequence is too short <3!: " << xname << endl; continue; } NewArray(&xa, xlen, 3); seqx = new char[xlen + 1]; secx = new char[xlen + 1]; xlen = read_PDB(PDB_lines1[chain_i], xa, seqx, resi_vec1, read_resi); if (mirror_opt) for (r = 0; r < xlen; r++) xa[r][2] = -xa[r][2]; if (mol_vec1[chain_i] > 0) make_sec(seqx, xa, xlen, secx, atom_opt); else make_sec(xa, xlen, secx); // secondary structure assignment for (j = (dir_opt.size() > 0) * (i + 1); j < chain2_list.size(); j++) { if (dirpair_opt.size() && j != i) continue; /* parse chain 2 */ if (PDB_lines2.size() == 0) { yname = chain2_list[j]; ychainnum = get_PDB_lines(yname, PDB_lines2, chainID_list2, mol_vec2, ter_opt, infmt2_opt, atom_opt, autojustify, split_opt, het_opt, chain2parse2, model2parse2); if (!ychainnum) { cerr << "Warning! Cannot parse file: " << yname << ". Chain number 0." << endl; continue; } } for (chain_j = 0; chain_j < ychainnum; chain_j++) { ylen = PDB_lines2[chain_j].size(); if (mol_opt == "RNA") mol_vec2[chain_j] = 1; else if (mol_opt == "protein") mol_vec2[chain_j] = -1; if (!ylen) { cerr << "Warning! Cannot parse file: " << yname << ". Chain length 0." << endl; continue; } else if (ylen < 3) { cerr << "Sequence is too short <3!: " << yname << endl; continue; } NewArray(&ya, ylen, 3); seqy = new char[ylen + 1]; secy = new char[ylen + 1]; ylen = read_PDB(PDB_lines2[chain_j], ya, seqy, resi_vec2, read_resi); if (mol_vec2[chain_j] > 0) make_sec(seqy, ya, ylen, secy, atom_opt); else make_sec(ya, ylen, secy); if (byresi_opt) extract_aln_from_resi(sequence, seqx, seqy, resi_vec1, resi_vec2, byresi_opt); /* declare variable specific to this pair of TMalign */ double t0[3], u0[3][3]; double TM1, TM2; double TM3, TM4, TM5; // for a_opt, u_opt, d_opt double d0_0, TM_0; double d0A, d0B, d0u, d0a; double d0_out = 5.0; string seqM, seqxA, seqyA; // for output alignment double rmsd0 = 0.0; int L_ali; // Aligned length in standard_TMscore double Liden = 0; double TM_ali, rmsd_ali; // TMscore and rmsd in standard_TMscore int n_ali = 0; int n_ali8 = 0; bool force_fast_opt = (getmin(xlen, ylen) > 1500) ? true : fast_opt; vector do_vec; /* entry function for structure alignment */ if (cp_opt) CPalign_main( xa, ya, seqx, seqy, secx, secy, t0, u0, TM1, TM2, TM3, TM4, TM5, d0_0, TM_0, d0A, d0B, d0u, d0a, d0_out, seqM, seqxA, seqyA, do_vec, rmsd0, L_ali, Liden, TM_ali, rmsd_ali, n_ali, n_ali8, xlen, ylen, sequence, Lnorm_ass, d0_scale, i_opt, a_opt, u_opt, d_opt, force_fast_opt, mol_vec1[chain_i] + mol_vec2[chain_j], TMcut); else if (se_opt) { int *invmap = new int[ylen + 1]; u0[0][0] = u0[1][1] = u0[2][2] = 1; u0[0][1] = u0[0][2] = u0[1][0] = u0[1][2] = u0[2][0] = u0[2][1] = t0[0] = t0[1] = t0[2] = 0; se_main(xa, ya, seqx, seqy, TM1, TM2, TM3, TM4, TM5, d0_0, TM_0, d0A, d0B, d0u, d0a, d0_out, seqM, seqxA, seqyA, do_vec, rmsd0, L_ali, Liden, TM_ali, rmsd_ali, n_ali, n_ali8, xlen, ylen, sequence, Lnorm_ass, d0_scale, i_opt, a_opt, u_opt, d_opt, mol_vec1[chain_i] + mol_vec2[chain_j], outfmt_opt, invmap); if (outfmt_opt >= 2) { Liden = L_ali = 0; int r1, r2; for (r2 = 0; r2 < ylen; r2++) { r1 = invmap[r2]; if (r1 < 0) continue; L_ali += 1; Liden += (seqx[r1] == seqy[r2]); } } delete[] invmap; } else TMalign_main( xa, ya, seqx, seqy, secx, secy, t0, u0, TM1, TM2, TM3, TM4, TM5, d0_0, TM_0, d0A, d0B, d0u, d0a, d0_out, seqM, seqxA, seqyA, do_vec, rmsd0, L_ali, Liden, TM_ali, rmsd_ali, n_ali, n_ali8, xlen, ylen, sequence, Lnorm_ass, d0_scale, i_opt, a_opt, u_opt, d_opt, force_fast_opt, mol_vec1[chain_i] + mol_vec2[chain_j], TMcut, 0); /* print result */ if (outfmt_opt == 0) print_version(); int left_num = 0; int right_num = 0; int left_aln_num = 0; int right_aln_num = 0; bool after_cp = false; if (cp_opt) after_cp = output_cp( xname.substr(dir1_opt.size() + dir_opt.size()), yname.substr(dir2_opt.size() + dir_opt.size()), seqxA, seqyA, outfmt_opt, left_num, right_num, left_aln_num, right_aln_num); output_results( xname.substr(dir1_opt.size() + dir_opt.size() + dirpair_opt.size()), yname.substr(dir2_opt.size() + dir_opt.size() + dirpair_opt.size()), chainID_list1[chain_i], chainID_list2[chain_j], xlen, ylen, t0, u0, TM1, TM2, TM3, TM4, TM5, rmsd0, d0_out, seqM.c_str(), seqxA.c_str(), seqyA.c_str(), Liden, n_ali8, L_ali, TM_ali, rmsd_ali, TM_0, d0_0, d0A, d0B, Lnorm_ass, d0_scale, d0a, d0u, (m_opt ? fname_matrix : "").c_str(), outfmt_opt, ter_opt, false, split_opt, o_opt, fname_super, i_opt, a_opt, u_opt, d_opt, mirror_opt, resi_vec1, resi_vec2); if (do_opt || (cp_opt && outfmt_opt <= 0)) { cout << "###############\t###############\t#########" << endl; cout << "#Aligned atom 1\tAligned atom 2 \tDistance#" << endl; size_t r1 = right_num; size_t r2 = 0; size_t r; int postcp = 0; for (r = 0; r < seqxA.size(); r++) { r1 += seqxA[r] != '-'; r2 += seqyA[r] != '-'; if (seqxA[r] == '*') { cout << "###### Circular\tPermutation ###\t#########\n"; r1 = 0; postcp = 1; } else if (seqxA[r] != '-' && seqyA[r] != '-') { cout << PDB_lines1[chain_i][r1 - 1].substr(12, 15) << '\t' << PDB_lines2[chain_j][r2 - 1].substr(12, 15) << '\t' << setw(9) << setiosflags(ios::fixed) << setprecision(3) << do_vec[r - postcp] << '\n'; } } cout << "###############\t###############\t#########" << endl; } /* Done! Free memory */ seqM.clear(); seqxA.clear(); seqyA.clear(); DeleteArray(&ya, ylen); delete[] seqy; delete[] secy; resi_vec2.clear(); do_vec.clear(); } // chain_j if (chain2_list.size() > 1) { yname.clear(); for (chain_j = 0; chain_j < ychainnum; chain_j++) PDB_lines2[chain_j].clear(); PDB_lines2.clear(); chainID_list2.clear(); mol_vec2.clear(); } } // j PDB_lines1[chain_i].clear(); DeleteArray(&xa, xlen); delete[] seqx; delete[] secx; resi_vec1.clear(); } // chain_i xname.clear(); PDB_lines1.clear(); chainID_list1.clear(); mol_vec1.clear(); } // i if (chain2_list.size() == 1) { yname.clear(); for (chain_j = 0; chain_j < ychainnum; chain_j++) PDB_lines2[chain_j].clear(); PDB_lines2.clear(); resi_vec2.clear(); chainID_list2.clear(); mol_vec2.clear(); } return 0; } /* MMalign if more than two chains. TMalign if only one chain */ int MMalign(const string &xname, const string &yname, const string &fname_super, const string &fname_lign, const string &fname_matrix, vector &sequence, const double d0_scale, const bool m_opt, const int o_opt, const int a_opt, const bool d_opt, const bool full_opt, const double TMcut, const int infmt1_opt, const int infmt2_opt, const int ter_opt, const int split_opt, const int outfmt_opt, bool fast_opt, const int mirror_opt, const int het_opt, const string &atom_opt, const bool autojustify, const string &mol_opt, const string &dir1_opt, const string &dir2_opt, const vector &chain2parse1, const vector &chain2parse2, const vector &model2parse1, const vector &model2parse2, const vector &chain1_list, const vector &chain2_list, const int byresi_opt, const string &chainmapfile, const bool se_opt) { /* declare previously global variables */ vector > > xa_vec;// structure of complex1 vector > > ya_vec;// structure of complex2 vector > seqx_vec; // sequence of complex1 vector > seqy_vec; // sequence of complex2 vector > secx_vec; // secondary structure of complex1 vector > secy_vec; // secondary structure of complex2 vector mol_vec1; // molecule type of complex1, RNA if >0 vector mol_vec2; // molecule type of complex2, RNA if >0 vector chainID_list1; // list of chainID1 vector chainID_list2; // list of chainID2 vector xlen_vec; // length of complex1 vector ylen_vec; // length of complex2 int i, j; // chain index int xlen, ylen; // chain length double **xa, **ya; // structure of single chain char *seqx, *seqy; // for the protein sequence char *secx, *secy; // for the secondary structure int xlen_aa, ylen_aa; // total length of protein int xlen_na, ylen_na; // total length of RNA/DNA vector resi_vec1; // residue index for chain1 vector resi_vec2; // residue index for chain2 /* parse complex */ parse_chain_list(chain1_list, xa_vec, seqx_vec, secx_vec, mol_vec1, xlen_vec, chainID_list1, ter_opt, split_opt, mol_opt, infmt1_opt, atom_opt, autojustify, mirror_opt, het_opt, xlen_aa, xlen_na, o_opt, resi_vec1, chain2parse1, model2parse1); if (xa_vec.size() == 0) PrintErrorAndQuit("ERROR! 0 chain in complex 1"); parse_chain_list(chain2_list, ya_vec, seqy_vec, secy_vec, mol_vec2, ylen_vec, chainID_list2, ter_opt, split_opt, mol_opt, infmt2_opt, atom_opt, autojustify, 0, het_opt, ylen_aa, ylen_na, o_opt, resi_vec2, chain2parse2, model2parse2); if (ya_vec.size() == 0) PrintErrorAndQuit("ERROR! 0 chain in complex 2"); int len_aa = getmin(xlen_aa, ylen_aa); int len_na = getmin(xlen_na, ylen_na); if (a_opt) { len_aa = (xlen_aa + ylen_aa) / 2; len_na = (xlen_na + ylen_na) / 2; } int i_opt = 0; if (byresi_opt) i_opt = 3; map chainmap; if (chainmapfile.size()) { string line; int chainidx1, chainidx2; vector line_vec; ifstream fin; bool fromStdin = (chainmapfile == "-"); if (!fromStdin) fin.open(chainmapfile.c_str()); while (fromStdin ? cin.good() : fin.good()) { if (fromStdin) getline(cin, line); else getline(fin, line); if (line.size() == 0 || line[0] == '#') continue; split(line, line_vec, '\t'); if (line_vec.size() == 2) { chainidx1 = -1; chainidx2 = -1; for (i = 0; i < chainID_list1.size(); i++) { if (line_vec[0] == chainID_list1[i] || ":" + line_vec[0] == chainID_list1[i] || ":1," + line_vec[0] == chainID_list1[i]) { chainidx1 = i; break; } } for (i = 0; i < chainID_list2.size(); i++) { if (line_vec[1] == chainID_list2[i] || ":" + line_vec[1] == chainID_list2[i] || ":1," + line_vec[1] == chainID_list2[i]) { chainidx2 = i; break; } } if (chainidx1 >= 0 && chainidx2 >= 0) { if (chainmap.count(chainidx1)) cerr << "ERROR! " << line_vec[0] << " already mapped" << endl; chainmap[chainidx1] = chainidx2; } else cerr << "ERROR! Cannot map " << line << endl; } else cerr << "ERROR! Cannot map " << line << endl; for (i = 0; i < line_vec.size(); i++) line_vec[i].clear(); line_vec.clear(); } if (!fromStdin) fin.close(); if (chainmap.size() == 0) cerr << "ERROR! cannot map any chain pair from " << chainmapfile << endl; } /* perform monomer alignment if there is only one chain */ if (xa_vec.size() == 1 && ya_vec.size() == 1) { xlen = xlen_vec[0]; ylen = ylen_vec[0]; seqx = new char[xlen + 1]; seqy = new char[ylen + 1]; secx = new char[xlen + 1]; secy = new char[ylen + 1]; NewArray(&xa, xlen, 3); NewArray(&ya, ylen, 3); copy_chain_data(xa_vec[0], seqx_vec[0], secx_vec[0], xlen, xa, seqx, secx); copy_chain_data(ya_vec[0], seqy_vec[0], secy_vec[0], ylen, ya, seqy, secy); /* declare variable specific to this pair of TMalign */ double t0[3], u0[3][3]; double TM1, TM2; double TM3, TM4, TM5; // for a_opt, u_opt, d_opt double d0_0, TM_0; double d0A, d0B, d0u, d0a; double d0_out = 5.0; string seqM, seqxA, seqyA; // for output alignment double rmsd0 = 0.0; int L_ali; // Aligned length in standard_TMscore double Liden = 0; double TM_ali, rmsd_ali; // TMscore and rmsd in standard_TMscore int n_ali = 0; int n_ali8 = 0; vector do_vec; if (byresi_opt) extract_aln_from_resi(sequence, seqx, seqy, resi_vec1, resi_vec2, byresi_opt); /* entry function for structure alignment */ if (se_opt) { int *invmap = new int[ylen + 1]; u0[0][0] = u0[1][1] = u0[2][2] = 1; u0[0][1] = u0[0][2] = u0[1][0] = u0[1][2] = u0[2][0] = u0[2][1] = t0[0] = t0[1] = t0[2] = 0; se_main(xa, ya, seqx, seqy, TM1, TM2, TM3, TM4, TM5, d0_0, TM_0, d0A, d0B, d0u, d0a, d0_out, seqM, seqxA, seqyA, do_vec, rmsd0, L_ali, Liden, TM_ali, rmsd_ali, n_ali, n_ali8, xlen, ylen, sequence, 0, d0_scale, i_opt, a_opt, false, d_opt, mol_vec1[0] + mol_vec2[0], outfmt_opt, invmap); if (outfmt_opt >= 2) { Liden = L_ali = 0; int r1, r2; for (r2 = 0; r2 < ylen; r2++) { r1 = invmap[r2]; if (r1 < 0) continue; L_ali += 1; Liden += (seqx[r1] == seqy[r2]); } } delete[] invmap; } else TMalign_main(xa, ya, seqx, seqy, secx, secy, t0, u0, TM1, TM2, TM3, TM4, TM5, d0_0, TM_0, d0A, d0B, d0u, d0a, d0_out, seqM, seqxA, seqyA, do_vec, rmsd0, L_ali, Liden, TM_ali, rmsd_ali, n_ali, n_ali8, xlen, ylen, sequence, 0, d0_scale, i_opt, a_opt, false, d_opt, fast_opt, mol_vec1[0] + mol_vec2[0], TMcut, 0); /* print result */ if (outfmt_opt == 0) print_version(); output_results( xname.substr(dir1_opt.size()), yname.substr(dir2_opt.size()), chainID_list1[0], chainID_list2[0], xlen, ylen, t0, u0, TM1, TM2, TM3, TM4, TM5, rmsd0, d0_out, seqM.c_str(), seqxA.c_str(), seqyA.c_str(), Liden, n_ali8, L_ali, TM_ali, rmsd_ali, TM_0, d0_0, d0A, d0B, 0, d0_scale, d0a, d0u, (m_opt ? fname_matrix : "").c_str(), outfmt_opt, ter_opt, true, split_opt, o_opt, fname_super, 0, a_opt, false, d_opt, mirror_opt, resi_vec1, resi_vec2); /* clean up */ seqM.clear(); seqxA.clear(); seqyA.clear(); delete[] seqx; delete[] seqy; delete[] secx; delete[] secy; DeleteArray(&xa, xlen); DeleteArray(&ya, ylen); do_vec.clear(); vector > >().swap(xa_vec); // structure of complex1 vector > >().swap(ya_vec); // structure of complex2 vector >().swap(seqx_vec); // sequence of complex1 vector >().swap(seqy_vec); // sequence of complex2 vector >().swap(secx_vec); // secondary structure of complex1 vector >().swap(secy_vec); // secondary structure of complex2 mol_vec1.clear(); // molecule type of complex1, RNA if >0 mol_vec2.clear(); // molecule type of complex2, RNA if >0 chainID_list1.clear(); // list of chainID1 chainID_list2.clear(); // list of chainID2 xlen_vec.clear(); // length of complex1 ylen_vec.clear(); // length of complex2 return 0; } /* declare TM-score tables */ int chain1_num = xa_vec.size(); int chain2_num = ya_vec.size(); int chain_num = MAX(chain1_num, chain2_num); vector tmp_str_vec(chain2_num, ""); double **TMave_mat; double **ut_mat; // rotation matrices for all-against-all alignment int ui,uj,ut_idx; NewArray(&TMave_mat,chain1_num,chain2_num); NewArray(&ut_mat,chain1_num*chain2_num,4*3); vector >seqxA_mat(chain1_num,tmp_str_vec); vector > seqM_mat(chain1_num,tmp_str_vec); vector >seqyA_mat(chain1_num,tmp_str_vec); double maxTMmono = -1; int maxTMmono_i, maxTMmono_j; /* get all-against-all alignment */ if (len_aa + len_na > 500) fast_opt = true; for (i = 0; i < chain1_num; i++) { xlen = xlen_vec[i]; if (xlen < 3) { for (j=0;j do_vec; int Lnorm_tmp = len_aa; if (mol_vec1[i] + mol_vec2[j] > 0) Lnorm_tmp = len_na; if (byresi_opt) { int total_aln = extract_aln_from_resi(sequence, seqx, seqy, resi_vec1, resi_vec2, xlen_vec, ylen_vec, i, j, byresi_opt); seqxA_mat[i][j] = sequence[0]; seqyA_mat[i][j] = sequence[1]; if (total_aln > xlen + ylen - 3) { for (ui=0;ui<3;ui++) for (uj=0;uj<3;uj++) ut_mat[ut_idx][ui*3+uj]=(ui==uj)?1:0; for (uj=0;uj<3;uj++) ut_mat[ut_idx][9+uj]=0; TMave_mat[i][j]=0; seqM.clear(); seqxA.clear(); seqyA.clear(); delete[] seqy; delete[] secy; DeleteArray(&ya, ylen); continue; } } /* entry function for structure alignment */ if (se_opt) { int *invmap = new int[ylen + 1]; u0[0][0] = u0[1][1] = u0[2][2] = 1; u0[0][1] = u0[0][2] = u0[1][0] = u0[1][2] = u0[2][0] = u0[2][1] = t0[0] = t0[1] = t0[2] = 0; se_main(xa, ya, seqx, seqy, TM1, TM2, TM3, TM4, TM5, d0_0, TM_0, d0A, d0B, d0u, d0a, d0_out, seqM, seqxA, seqyA, do_vec, rmsd0, L_ali, Liden, TM_ali, rmsd_ali, n_ali, n_ali8, xlen, ylen, sequence, Lnorm_tmp, d0_scale, i_opt, false, true, false, mol_vec1[i] + mol_vec2[j], outfmt_opt, invmap); if (outfmt_opt >= 2) { Liden = L_ali = 0; int r1, r2; for (r2 = 0; r2 < ylen; r2++) { r1 = invmap[r2]; if (r1 < 0) continue; L_ali += 1; Liden += (seqx[r1] == seqy[r2]); } } delete[] invmap; } else TMalign_main(xa, ya, seqx, seqy, secx, secy, t0, u0, TM1, TM2, TM3, TM4, TM5, d0_0, TM_0, d0A, d0B, d0u, d0a, d0_out, seqM, seqxA, seqyA, do_vec, rmsd0, L_ali, Liden, TM_ali, rmsd_ali, n_ali, n_ali8, xlen, ylen, sequence, Lnorm_tmp, d0_scale, i_opt, false, true, false, fast_opt, mol_vec1[i] + mol_vec2[j], TMcut, 0); /* store result */ for (ui=0;ui<3;ui++) for (uj=0;uj<3;uj++) ut_mat[ut_idx][ui*3+uj]=u0[ui][uj]; for (uj=0;uj<3;uj++) ut_mat[ut_idx][9+uj]=t0[uj]; seqxA_mat[i][j]=seqxA; seqyA_mat[i][j]=seqyA; TMave_mat[i][j]=TM4*Lnorm_tmp; if (TMave_mat[i][j]>maxTMmono) { maxTMmono = TMave_mat[i][j]; maxTMmono_i = i; maxTMmono_j = j; } /* clean up */ seqM.clear(); seqxA.clear(); seqyA.clear(); delete[] seqy; delete[] secy; DeleteArray(&ya, ylen); do_vec.clear(); } delete[] seqx; delete[] secx; DeleteArray(&xa, xlen); } /* calculate initial chain-chain assignment */ int *assign1_list; // value is index of assigned chain2 int *assign2_list; // value is index of assigned chain1 assign1_list = new int[chain1_num]; assign2_list = new int[chain2_num]; double total_score = enhanced_greedy_search(TMave_mat, assign1_list, assign2_list, chain1_num, chain2_num); if (total_score <= 0) PrintErrorAndQuit("ERROR! No assignable chain"); /* refine alignment for large oligomers */ int aln_chain_num = count_assign_pair(assign1_list, chain1_num); bool is_oligomer = (aln_chain_num >= 3); if (aln_chain_num == 2 && chainmap.size() == 0 && !se_opt) // dimer alignment { int na_chain_num1, na_chain_num2, aa_chain_num1, aa_chain_num2; count_na_aa_chain_num(na_chain_num1, aa_chain_num1, mol_vec1); count_na_aa_chain_num(na_chain_num2, aa_chain_num2, mol_vec2); /* align protein-RNA hybrid dimer to another hybrid dimer */ if (na_chain_num1 == 1 && na_chain_num2 == 1 && aa_chain_num1 == 1 && aa_chain_num2 == 1) is_oligomer = false; /* align pure protein dimer or pure RNA dimer */ else if ((getmin(na_chain_num1, na_chain_num2) == 0 && aa_chain_num1 == 2 && aa_chain_num2 == 2) || (getmin(aa_chain_num1, aa_chain_num2) == 0 && na_chain_num1 == 2 && na_chain_num2 == 2)) { adjust_dimer_assignment(xa_vec, ya_vec, xlen_vec, ylen_vec, mol_vec1, mol_vec2, assign1_list, assign2_list, seqxA_mat, seqyA_mat); is_oligomer = false; // cannot refiner further } else is_oligomer = true; /* align oligomers to dimer */ } if ((aln_chain_num >= 3 || is_oligomer) && chainmap.size() == 0 && !se_opt) // oligomer alignment { /* extract centroid coordinates */ double **xcentroids; double **ycentroids; NewArray(&xcentroids, chain1_num, 3); NewArray(&ycentroids, chain2_num, 3); double d0MM = getmin( calculate_centroids(xa_vec, chain1_num, xcentroids), calculate_centroids(ya_vec, chain2_num, ycentroids)); /* refine enhanced greedy search with centroid superposition */ // double het_deg=check_heterooligomer(TMave_mat, chain1_num, chain2_num); homo_refined_greedy_search(TMave_mat, assign1_list, assign2_list, chain1_num, chain2_num, xcentroids, ycentroids, d0MM, len_aa + len_na, ut_mat); if (chain1_num <= chain2_num) { hetero_refined_greedy_search(TMave_mat, assign1_list, assign2_list, chain1_num, chain2_num, xcentroids, ycentroids, d0MM, len_aa + len_na); } else { hetero_refined_greedy_search(TMave_mat, assign2_list, assign1_list, chain2_num, chain1_num, ycentroids, xcentroids, d0MM, len_aa + len_na); } /* clean up */ DeleteArray(&xcentroids, chain1_num); DeleteArray(&ycentroids, chain2_num); } /* store initial assignment */ int init_pair_num = count_assign_pair(assign1_list, chain1_num); int *assign1_init, *assign2_init; assign1_init = new int[chain1_num]; assign2_init = new int[chain2_num]; double **TMave_init; NewArray(&TMave_init, chain1_num, chain2_num); vector > seqxA_init(chain1_num, tmp_str_vec); vector > seqyA_init(chain1_num, tmp_str_vec); vector sequence_init; copy_chain_assign_data(chain1_num, chain2_num, sequence_init, seqxA_mat, seqyA_mat, assign1_list, assign2_list, TMave_mat, seqxA_init, seqyA_init, assign1_init, assign2_init, TMave_init); /* perform iterative alignment */ double max_total_score = 0; // ignore old total_score because previous // score was from monomeric chain superpositions int max_iter = 5 - (int)((len_aa + len_na) / 200); if (max_iter < 2) max_iter = 2; // if (byresi_opt==0) if (!se_opt) MMalign_iter(max_total_score, max_iter, xa_vec, ya_vec, seqx_vec, seqy_vec, secx_vec, secy_vec, mol_vec1, mol_vec2, xlen_vec, ylen_vec, xa, ya, seqx, seqy, secx, secy, len_aa, len_na, chain1_num, chain2_num, TMave_mat, seqxA_mat, seqyA_mat, assign1_list, assign2_list, sequence, d0_scale, fast_opt, chainmap, byresi_opt); if (byresi_opt && aln_chain_num >= 4 && is_oligomer && chainmap.size() == 0 && !se_opt) // oligomer alignment { MMalign_final(xname.substr(dir1_opt.size()), yname.substr(dir2_opt.size()), chainID_list1, chainID_list2, fname_super, fname_lign, fname_matrix, xa_vec, ya_vec, seqx_vec, seqy_vec, secx_vec, secy_vec, mol_vec1, mol_vec2, xlen_vec, ylen_vec, xa, ya, seqx, seqy, secx, secy, len_aa, len_na, chain1_num, chain2_num, TMave_mat, seqxA_mat, seqM_mat, seqyA_mat, assign1_list, assign2_list, sequence, d0_scale, 1, 0, 5, ter_opt, split_opt, 0, 0, true, true, mirror_opt, resi_vec1, resi_vec2); /* extract centroid coordinates */ double **xcentroids; double **ycentroids; NewArray(&xcentroids, chain1_num, 3); NewArray(&ycentroids, chain2_num, 3); double d0MM = getmin( calculate_centroids(xa_vec, chain1_num, xcentroids), calculate_centroids(ya_vec, chain2_num, ycentroids)); /* refine enhanced greedy search with centroid superposition */ // double het_deg=check_heterooligomer(TMave_mat, chain1_num, chain2_num); homo_refined_greedy_search(TMave_mat, assign1_list, assign2_list, chain1_num, chain2_num, xcentroids, ycentroids, d0MM, len_aa + len_na, ut_mat); hetero_refined_greedy_search(TMave_mat, assign1_list, assign2_list, chain1_num, chain2_num, xcentroids, ycentroids, d0MM, len_aa + len_na); /* clean up */ DeleteArray(&xcentroids, chain1_num); DeleteArray(&ycentroids, chain2_num); } /* sometime MMalign_iter is even worse than monomer alignment */ if (byresi_opt == 0 && max_total_score < maxTMmono) { copy_chain_assign_data(chain1_num, chain2_num, sequence, seqxA_init, seqyA_init, assign1_init, assign2_init, TMave_init, seqxA_mat, seqyA_mat, assign1_list, assign2_list, TMave_mat); for (i = 0; i < chain1_num; i++) { if (i != maxTMmono_i) assign1_list[i] = -1; else assign1_list[i] = maxTMmono_j; } for (j = 0; j < chain2_num; j++) { if (j != maxTMmono_j) assign2_list[j] = -1; else assign2_list[j] = maxTMmono_i; } sequence[0] = seqxA_mat[maxTMmono_i][maxTMmono_j]; sequence[1] = seqyA_mat[maxTMmono_i][maxTMmono_j]; max_total_score = maxTMmono; MMalign_iter(max_total_score, max_iter, xa_vec, ya_vec, seqx_vec, seqy_vec, secx_vec, secy_vec, mol_vec1, mol_vec2, xlen_vec, ylen_vec, xa, ya, seqx, seqy, secx, secy, len_aa, len_na, chain1_num, chain2_num, TMave_mat, seqxA_mat, seqyA_mat, assign1_list, assign2_list, sequence, d0_scale, fast_opt, chainmap); } /* perform cross chain alignment * in some cases, this leads to dramatic improvement, esp for homodimer */ int iter_pair_num = count_assign_pair(assign1_list, chain1_num); if (iter_pair_num >= init_pair_num) copy_chain_assign_data( chain1_num, chain2_num, sequence_init, seqxA_mat, seqyA_mat, assign1_list, assign2_list, TMave_mat, seqxA_init, seqyA_init, assign1_init, assign2_init, TMave_init); double max_total_score_cross = max_total_score; if (byresi_opt == 0 && len_aa + len_na < 10000) { MMalign_dimer(max_total_score_cross, xa_vec, ya_vec, seqx_vec, seqy_vec, secx_vec, secy_vec, mol_vec1, mol_vec2, xlen_vec, ylen_vec, xa, ya, seqx, seqy, secx, secy, len_aa, len_na, chain1_num, chain2_num, TMave_init, seqxA_init, seqyA_init, assign1_init, assign2_init, sequence_init, d0_scale, fast_opt); if (max_total_score_cross > max_total_score) { max_total_score = max_total_score_cross; copy_chain_assign_data(chain1_num, chain2_num, sequence, seqxA_init, seqyA_init, assign1_init, assign2_init, TMave_init, seqxA_mat, seqyA_mat, assign1_list, assign2_list, TMave_mat); } } /* final alignment */ if (outfmt_opt == 0) print_version(); if (se_opt) MMalign_se_final(xname.substr(dir1_opt.size()), yname.substr(dir2_opt.size()), chainID_list1, chainID_list2, fname_super, fname_lign, fname_matrix, xa_vec, ya_vec, seqx_vec, seqy_vec, secx_vec, secy_vec, mol_vec1, mol_vec2, xlen_vec, ylen_vec, xa, ya, seqx, seqy, secx, secy, len_aa, len_na, chain1_num, chain2_num, TMave_mat, seqxA_mat, seqM_mat, seqyA_mat, assign1_list, assign2_list, sequence, d0_scale, m_opt, o_opt, outfmt_opt, ter_opt, split_opt, a_opt, d_opt, fast_opt, full_opt, mirror_opt, resi_vec1, resi_vec2); else MMalign_final(xname.substr(dir1_opt.size()), yname.substr(dir2_opt.size()), chainID_list1, chainID_list2, fname_super, fname_lign, fname_matrix, xa_vec, ya_vec, seqx_vec, seqy_vec, secx_vec, secy_vec, mol_vec1, mol_vec2, xlen_vec, ylen_vec, xa, ya, seqx, seqy, secx, secy, len_aa, len_na, chain1_num, chain2_num, TMave_mat, seqxA_mat, seqM_mat, seqyA_mat, assign1_list, assign2_list, sequence, d0_scale, m_opt, o_opt, outfmt_opt, ter_opt, split_opt, a_opt, d_opt, fast_opt, full_opt, mirror_opt, resi_vec1, resi_vec2); /* clean up everything */ delete [] assign1_list; delete [] assign2_list; DeleteArray(&TMave_mat,chain1_num); DeleteArray(&ut_mat, chain1_num*chain2_num); vector >().swap(seqxA_mat); vector >().swap(seqM_mat); vector >().swap(seqyA_mat); vector().swap(tmp_str_vec); delete[] assign1_init; delete[] assign2_init; DeleteArray(&TMave_init, chain1_num); vector >().swap(seqxA_init); vector >().swap(seqyA_init); vector > >().swap(xa_vec); // structure of complex1 vector > >().swap(ya_vec); // structure of complex2 vector >().swap(seqx_vec); // sequence of complex1 vector >().swap(seqy_vec); // sequence of complex2 vector >().swap(secx_vec); // secondary structure of complex1 vector >().swap(secy_vec); // secondary structure of complex2 mol_vec1.clear(); // molecule type of complex1, RNA if >0 mol_vec2.clear(); // molecule type of complex2, RNA if >0 vector().swap(chainID_list1); // list of chainID1 vector().swap(chainID_list2); // list of chainID2 xlen_vec.clear(); // length of complex1 ylen_vec.clear(); // length of complex2 vector().swap(resi_vec1); // residue index for chain1 vector().swap(resi_vec2); // residue index for chain2 map().swap(chainmap); return 1; } /* alignment individual chains to a complex. */ int MMdock(const string &xname, const string &yname, const string &fname_super, const string &fname_matrix, vector &sequence, const double Lnorm_ass, const double d0_scale, const bool m_opt, const int o_opt, const int a_opt, const bool u_opt, const bool d_opt, const double TMcut, const int infmt1_opt, const int infmt2_opt, const int ter_opt, const int split_opt, const int outfmt_opt, bool fast_opt, const int mirror_opt, const int het_opt, const string &atom_opt, const bool autojustify, const string &mol_opt, const string &dir1_opt, const string &dir2_opt, const vector &chain2parse1, const vector &chain2parse2, const vector &model2parse1, const vector &model2parse2, const vector &chain1_list, const vector &chain2_list, const bool do_opt) { /* declare previously global variables */ vector > > xa_vec; // structure of complex1 vector > > ya_vec; // structure of complex2 vector > seqx_vec; // sequence of complex1 vector > seqy_vec; // sequence of complex2 vector > secx_vec; // secondary structure of complex1 vector > secy_vec; // secondary structure of complex2 vector mol_vec1; // molecule type of complex1, RNA if >0 vector mol_vec2; // molecule type of complex2, RNA if >0 vector chainID_list1; // list of chainID1 vector chainID_list2; // list of chainID2 vector xlen_vec; // length of complex1 vector ylen_vec; // length of complex2 int i, j; // chain index int xlen, ylen; // chain length double **xa, **ya; // structure of single chain char *seqx, *seqy; // for the protein sequence char *secx, *secy; // for the secondary structure int xlen_aa, ylen_aa; // total length of protein int xlen_na, ylen_na; // total length of RNA/DNA vector resi_vec1; // residue index for chain1 vector resi_vec2; // residue index for chain2 /* parse complex */ parse_chain_list(chain1_list, xa_vec, seqx_vec, secx_vec, mol_vec1, xlen_vec, chainID_list1, ter_opt, split_opt, mol_opt, infmt1_opt, atom_opt, autojustify, mirror_opt, het_opt, xlen_aa, xlen_na, o_opt, resi_vec1, chain2parse1, model2parse1); if (xa_vec.size() == 0) PrintErrorAndQuit("ERROR! 0 individual chain"); parse_chain_list(chain2_list, ya_vec, seqy_vec, secy_vec, mol_vec2, ylen_vec, chainID_list2, ter_opt, split_opt, mol_opt, infmt2_opt, atom_opt, autojustify, 0, het_opt, ylen_aa, ylen_na, o_opt, resi_vec2, chain2parse2, model2parse2); if (xa_vec.size() > ya_vec.size()) PrintErrorAndQuit( "ERROR! more individual chains to align than number of chains in complex template"); int len_aa = getmin(xlen_aa, ylen_aa); int len_na = getmin(xlen_na, ylen_na); if (a_opt) { len_aa = (xlen_aa + ylen_aa) / 2; len_na = (xlen_na + ylen_na) / 2; } /* perform monomer alignment if there is only one chain */ if (xa_vec.size() == 1 && ya_vec.size() == 1) { xlen = xlen_vec[0]; ylen = ylen_vec[0]; seqx = new char[xlen + 1]; seqy = new char[ylen + 1]; secx = new char[xlen + 1]; secy = new char[ylen + 1]; NewArray(&xa, xlen, 3); NewArray(&ya, ylen, 3); copy_chain_data(xa_vec[0], seqx_vec[0], secx_vec[0], xlen, xa, seqx, secx); copy_chain_data(ya_vec[0], seqy_vec[0], secy_vec[0], ylen, ya, seqy, secy); /* declare variable specific to this pair of TMalign */ double t0[3], u0[3][3]; double TM1, TM2; double TM3, TM4, TM5; // for a_opt, u_opt, d_opt double d0_0, TM_0; double d0A, d0B, d0u, d0a; double d0_out = 5.0; string seqM, seqxA, seqyA; // for output alignment double rmsd0 = 0.0; int L_ali; // Aligned length in standard_TMscore double Liden = 0; double TM_ali, rmsd_ali; // TMscore and rmsd in standard_TMscore int n_ali = 0; int n_ali8 = 0; vector do_vec; /* entry function for structure alignment */ TMalign_main(xa, ya, seqx, seqy, secx, secy, t0, u0, TM1, TM2, TM3, TM4, TM5, d0_0, TM_0, d0A, d0B, d0u, d0a, d0_out, seqM, seqxA, seqyA, do_vec, rmsd0, L_ali, Liden, TM_ali, rmsd_ali, n_ali, n_ali8, xlen, ylen, sequence, Lnorm_ass, d0_scale, 0, a_opt, u_opt, d_opt, fast_opt, mol_vec1[0] + mol_vec2[0], TMcut, 0); /* print result */ output_results( xname.substr(dir1_opt.size()), yname.substr(dir2_opt.size()), chainID_list1[0], chainID_list2[0], xlen, ylen, t0, u0, TM1, TM2, TM3, TM4, TM5, rmsd0, d0_out, seqM.c_str(), seqxA.c_str(), seqyA.c_str(), Liden, n_ali8, L_ali, TM_ali, rmsd_ali, TM_0, d0_0, d0A, d0B, Lnorm_ass, d0_scale, d0a, d0u, (m_opt ? fname_matrix : "").c_str(), (outfmt_opt == 2 ? outfmt_opt : 3), ter_opt, true, split_opt, o_opt, fname_super, 0, a_opt, false, d_opt, mirror_opt, resi_vec1, resi_vec2); if (outfmt_opt == 2) printf("%s%s\t%s%s\t%.4f\n", xname.substr(dir1_opt.size()).c_str(), chainID_list1[0].c_str(), yname.substr(dir2_opt.size()).c_str(), chainID_list2[0].c_str(), sqrt((TM1 * TM1 + TM2 * TM2) / 2)); /* clean up */ seqM.clear(); seqxA.clear(); seqyA.clear(); delete[] seqx; delete[] seqy; delete[] secx; delete[] secy; DeleteArray(&xa, xlen); DeleteArray(&ya, ylen); do_vec.clear(); vector > >().swap(xa_vec); // structure of complex1 vector > >().swap(ya_vec); // structure of complex2 vector >().swap(seqx_vec); // sequence of complex1 vector >().swap(seqy_vec); // sequence of complex2 vector >().swap(secx_vec); // secondary structure of complex1 vector >().swap(secy_vec); // secondary structure of complex2 mol_vec1.clear(); // molecule type of complex1, RNA if >0 mol_vec2.clear(); // molecule type of complex2, RNA if >0 chainID_list1.clear(); // list of chainID1 chainID_list2.clear(); // list of chainID2 xlen_vec.clear(); // length of complex1 ylen_vec.clear(); // length of complex2 return 0; } /* declare TM-score tables */ int chain1_num = xa_vec.size(); int chain2_num = ya_vec.size(); vector tmp_str_vec(chain2_num, ""); double **TMave_mat; NewArray(&TMave_mat, chain1_num, chain2_num); vector > seqxA_mat(chain1_num, tmp_str_vec); vector > seqM_mat(chain1_num, tmp_str_vec); vector > seqyA_mat(chain1_num, tmp_str_vec); /* trimComplex */ vector > > ya_trim_vec; // structure of complex2 vector > seqy_trim_vec; // sequence of complex2 vector > secy_trim_vec; // secondary structure of complex2 vector ylen_trim_vec; // length of complex2 int Lchain_aa_max1 = 0; int Lchain_na_max1 = 0; for (i = 0; i < chain1_num; i++) { xlen = xlen_vec[i]; if (mol_vec1[i] > 0 && xlen > Lchain_na_max1) Lchain_na_max1 = xlen; else if (mol_vec1[i] <= 0 && xlen > Lchain_aa_max1) Lchain_aa_max1 = xlen; } int trim_chain_count = trimComplex(ya_trim_vec, seqy_trim_vec, secy_trim_vec, ylen_trim_vec, ya_vec, seqy_vec, secy_vec, ylen_vec, mol_vec2, Lchain_aa_max1, Lchain_na_max1); int ylen_trim; // chain length double **ya_trim; // structure of single chain char *seqy_trim; // for the protein sequence char *secy_trim; // for the secondary structure double **xt; /* get all-against-all alignment */ if (len_aa + len_na > 500) fast_opt = true; for (i = 0; i < chain1_num; i++) { xlen = xlen_vec[i]; if (xlen < 3) { for (j = 0; j < chain2_num; j++) TMave_mat[i][j] = -1; continue; } seqx = new char[xlen + 1]; secx = new char[xlen + 1]; NewArray(&xa, xlen, 3); copy_chain_data(xa_vec[i], seqx_vec[i], secx_vec[i], xlen, xa, seqx, secx); for (j = 0; j < chain2_num; j++) { if (mol_vec1[i] * mol_vec2[j] < 0) // no protein-RNA alignment { TMave_mat[i][j] = -1; continue; } ylen = ylen_vec[j]; if (ylen < 3) { TMave_mat[i][j] = -1; continue; } seqy = new char[ylen + 1]; secy = new char[ylen + 1]; NewArray(&ya, ylen, 3); copy_chain_data(ya_vec[j], seqy_vec[j], secy_vec[j], ylen, ya, seqy, secy); /* declare variable specific to this pair of TMalign */ double t0[3], u0[3][3]; double TM1, TM2; double TM3, TM4, TM5; // for a_opt, u_opt, d_opt double d0_0, TM_0; double d0A, d0B, d0u, d0a; double d0_out = 5.0; string seqM, seqxA, seqyA; // for output alignment double rmsd0 = 0.0; int L_ali; // Aligned length in standard_TMscore double Liden = 0; double TM_ali, rmsd_ali; // TMscore and rmsd in standard_TMscore int n_ali = 0; int n_ali8 = 0; vector do_vec; int Lnorm_tmp = len_aa; if (mol_vec1[i] + mol_vec2[j] > 0) Lnorm_tmp = len_na; /* entry function for structure alignment */ if (trim_chain_count && ylen_trim_vec[j] < ylen) { ylen_trim = ylen_trim_vec[j]; seqy_trim = new char[ylen_trim + 1]; secy_trim = new char[ylen_trim + 1]; NewArray(&ya_trim, ylen_trim, 3); copy_chain_data(ya_trim_vec[j], seqy_trim_vec[j], secy_trim_vec[j], ylen_trim, ya_trim, seqy_trim, secy_trim); TMalign_main(xa, ya_trim, seqx, seqy_trim, secx, secy_trim, t0, u0, TM1, TM2, TM3, TM4, TM5, d0_0, TM_0, d0A, d0B, d0u, d0a, d0_out, seqM, seqxA, seqyA, do_vec, rmsd0, L_ali, Liden, TM_ali, rmsd_ali, n_ali, n_ali8, xlen, ylen_trim, sequence, Lnorm_tmp, d0_scale, 0, false, true, false, fast_opt, mol_vec1[i] + mol_vec2[j], TMcut, 0); seqxA.clear(); seqyA.clear(); delete[] seqy_trim; delete[] secy_trim; DeleteArray(&ya_trim, ylen_trim); NewArray(&xt, xlen, 3); do_rotation(xa, xt, xlen, t0, u0); int *invmap = new int[ylen + 1]; se_main(xt, ya, seqx, seqy, TM1, TM2, TM3, TM4, TM5, d0_0, TM_0, d0A, d0B, d0u, d0a, d0_out, seqM, seqxA, seqyA, do_vec, rmsd0, L_ali, Liden, TM_ali, rmsd_ali, n_ali, n_ali8, xlen, ylen, sequence, Lnorm_tmp, d0_scale, 0, false, 2, false, mol_vec1[i] + mol_vec2[j], 1, invmap); delete[] invmap; if (sequence.size() < 2) sequence.push_back(""); if (sequence.size() < 2) sequence.push_back(""); sequence[0] = seqxA; sequence[1] = seqyA; TMalign_main(xt, ya, seqx, seqy, secx, secy, t0, u0, TM1, TM2, TM3, TM4, TM5, d0_0, TM_0, d0A, d0B, d0u, d0a, d0_out, seqM, seqxA, seqyA, do_vec, rmsd0, L_ali, Liden, TM_ali, rmsd_ali, n_ali, n_ali8, xlen, ylen, sequence, Lnorm_tmp, d0_scale, 2, false, true, false, fast_opt, mol_vec1[i] + mol_vec2[j], TMcut, 0); DeleteArray(&xt, xlen); } else { TMalign_main(xa, ya, seqx, seqy, secx, secy, t0, u0, TM1, TM2, TM3, TM4, TM5, d0_0, TM_0, d0A, d0B, d0u, d0a, d0_out, seqM, seqxA, seqyA, do_vec, rmsd0, L_ali, Liden, TM_ali, rmsd_ali, n_ali, n_ali8, xlen, ylen, sequence, Lnorm_tmp, d0_scale, 0, false, true, false, fast_opt, mol_vec1[i] + mol_vec2[j], TMcut, 0); } /* store result */ seqxA_mat[i][j] = seqxA; seqyA_mat[i][j] = seqyA; TMave_mat[i][j] = TM4 * Lnorm_tmp; /* clean up */ seqM.clear(); seqxA.clear(); seqyA.clear(); do_vec.clear(); delete[] seqy; delete[] secy; DeleteArray(&ya, ylen); } delete[] seqx; delete[] secx; DeleteArray(&xa, xlen); } vector > >().swap(ya_trim_vec); vector >().swap(seqy_trim_vec); vector >().swap(secy_trim_vec); vector().swap(ylen_trim_vec); /* calculate initial chain-chain assignment */ int *assign1_list; // value is index of assigned chain2 int *assign2_list; // value is index of assigned chain1 assign1_list = new int[chain1_num]; assign2_list = new int[chain2_num]; enhanced_greedy_search(TMave_mat, assign1_list, assign2_list, chain1_num, chain2_num); /* final alignment */ if (outfmt_opt == 0) print_version(); double **ut_mat; // rotation matrices for all-against-all alignment NewArray(&ut_mat, chain1_num, 4 * 3); int ui, uj; vector xname_vec; vector yname_vec; vector TM_vec; for (i = 0; i < chain1_num; i++) { j = assign1_list[i]; xname_vec.push_back(xname + chainID_list1[i]); if (j < 0) { cerr << "Warning! " << chainID_list1[i] << " cannot be alighed" << endl; for (ui = 0; ui < 3; ui++) { for (uj = 0; uj < 4; uj++) ut_mat[i][ui * 3 + uj] = 0; ut_mat[i][ui * 3 + ui] = 1; } yname_vec.push_back(yname); continue; } yname_vec.push_back(yname + chainID_list2[j]); xlen = xlen_vec[i]; seqx = new char[xlen + 1]; secx = new char[xlen + 1]; NewArray(&xa, xlen, 3); copy_chain_data(xa_vec[i], seqx_vec[i], secx_vec[i], xlen, xa, seqx, secx); ylen = ylen_vec[j]; seqy = new char[ylen + 1]; secy = new char[ylen + 1]; NewArray(&ya, ylen, 3); copy_chain_data(ya_vec[j], seqy_vec[j], secy_vec[j], ylen, ya, seqy, secy); /* declare variable specific to this pair of TMalign */ double t0[3], u0[3][3]; double TM1, TM2; double TM3, TM4, TM5; // for a_opt, u_opt, d_opt double d0_0, TM_0; double d0A, d0B, d0u, d0a; double d0_out = 5.0; string seqM, seqxA, seqyA; // for output alignment double rmsd0 = 0.0; int L_ali; // Aligned length in standard_TMscore double Liden = 0; double TM_ali, rmsd_ali; // TMscore and rmsd in standard_TMscore int n_ali = 0; int n_ali8 = 0; vector do_vec; int c; for (c = 0; c < sequence.size(); c++) sequence[c].clear(); sequence.clear(); sequence.push_back(seqxA_mat[i][j]); sequence.push_back(seqyA_mat[i][j]); /* entry function for structure alignment */ TMalign_main(xa, ya, seqx, seqy, secx, secy, t0, u0, TM1, TM2, TM3, TM4, TM5, d0_0, TM_0, d0A, d0B, d0u, d0a, d0_out, seqM, seqxA, seqyA, do_vec, rmsd0, L_ali, Liden, TM_ali, rmsd_ali, n_ali, n_ali8, xlen, ylen, sequence, Lnorm_ass, d0_scale, 3, a_opt, u_opt, d_opt, fast_opt, mol_vec1[i] + mol_vec2[j], -1, 0); for (ui = 0; ui < 3; ui++) for (uj = 0; uj < 3; uj++) ut_mat[i][ui * 3 + uj] = u0[ui][uj]; for (uj = 0; uj < 3; uj++) ut_mat[i][9 + uj] = t0[uj]; TM_vec.push_back(TM1); TM_vec.push_back(TM2); if (outfmt_opt < 2) output_results( xname.c_str(), yname.c_str(), chainID_list1[i], chainID_list2[j], xlen, ylen, t0, u0, TM1, TM2, TM3, TM4, TM5, rmsd0, d0_out, seqM.c_str(), seqxA.c_str(), seqyA.c_str(), Liden, n_ali8, L_ali, TM_ali, rmsd_ali, TM_0, d0_0, d0A, d0B, Lnorm_ass, d0_scale, d0a, d0u, "", outfmt_opt, ter_opt, false, split_opt, false, "", // o_opt, fname_super+chainID_list1[i], false, a_opt, u_opt, d_opt, mirror_opt, resi_vec1, resi_vec2); /* clean up */ seqM.clear(); seqxA.clear(); seqyA.clear(); delete[] seqy; delete[] secy; DeleteArray(&ya, ylen); delete[] seqx; delete[] secx; DeleteArray(&xa, xlen); do_vec.clear(); } if (outfmt_opt == 2) { double TM = 0; for (i = 0; i < TM_vec.size(); i++) TM += TM_vec[i] * TM_vec[i]; TM = sqrt(TM / TM_vec.size()); string query_name = xname; string template_name = yname; for (i = 0; i < chain1_num; i++) { j = assign1_list[i]; if (j < 0) continue; query_name += chainID_list1[i]; template_name += chainID_list2[j]; } printf("%s\t%s\t%.4f\n", query_name.c_str(), template_name.c_str(), TM); query_name.clear(); template_name.clear(); } if (m_opt) output_dock_rotation_matrix(fname_matrix.c_str(), xname_vec, yname_vec, ut_mat, assign1_list); if (o_opt) output_dock(chain1_list, ter_opt, split_opt, infmt1_opt, atom_opt, mirror_opt, ut_mat, fname_super); /* clean up everything */ vector().swap(TM_vec); vector().swap(xname_vec); vector().swap(yname_vec); delete[] assign1_list; delete[] assign2_list; DeleteArray(&TMave_mat, chain1_num); DeleteArray(&ut_mat, chain1_num); vector >().swap(seqxA_mat); vector >().swap(seqM_mat); vector >().swap(seqyA_mat); vector().swap(tmp_str_vec); vector > >().swap(xa_vec); // structure of complex1 vector > >().swap(ya_vec); // structure of complex2 vector >().swap(seqx_vec); // sequence of complex1 vector >().swap(seqy_vec); // sequence of complex2 vector >().swap(secx_vec); // secondary structure of complex1 vector >().swap(secy_vec); // secondary structure of complex2 mol_vec1.clear(); // molecule type of complex1, RNA if >0 mol_vec2.clear(); // molecule type of complex2, RNA if >0 vector().swap(chainID_list1); // list of chainID1 vector().swap(chainID_list2); // list of chainID2 xlen_vec.clear(); // length of complex1 ylen_vec.clear(); // length of complex2 return 1; } int mTMalign(string &xname, string &yname, const string &fname_super, const string &fname_matrix, vector &sequence, double Lnorm_ass, const double d0_scale, const bool m_opt, const int i_opt, const int o_opt, const int a_opt, bool u_opt, const bool d_opt, const bool full_opt, const double TMcut, const int infmt_opt, const int ter_opt, const int split_opt, const int outfmt_opt, bool fast_opt, const int het_opt, const string &atom_opt, const bool autojustify, const string &mol_opt, const string &dir_opt, const int byresi_opt, const vector &chain_list, const vector &chain2parse, const vector &model2parse, const bool se_opt) { /* declare previously global variables */ vector > > a_vec; // atomic structure vector > > ua_vec; // unchanged atomic structure vector > seq_vec; // sequence of complex vector > sec_vec; // secondary structure of complex vector mol_vec; // molecule type of complex1, RNA if >0 vector chainID_list; // list of chainID vector len_vec; // length of complex int i, j; // chain index int xlen, ylen; // chain length double **xa, **ya; // structure of single chain char *seqx, *seqy; // for the protein sequence char *secx, *secy; // for the secondary structure int len_aa, len_na; // total length of protein and RNA/DNA vector resi_vec; // residue index for chain /* parse chain list */ parse_chain_list(chain_list, a_vec, seq_vec, sec_vec, mol_vec, len_vec, chainID_list, ter_opt, split_opt, mol_opt, infmt_opt, atom_opt, autojustify, false, het_opt, len_aa, len_na, o_opt, resi_vec, chain2parse, model2parse); int chain_num = a_vec.size(); if (chain_num <= 1) PrintErrorAndQuit("ERROR! <2 chains for multiple alignment"); if (m_opt || o_opt) for (i = 0; i < chain_num; i++) ua_vec.push_back(a_vec[i]); int mol_type = 0; int total_len = 0; xlen = 0; for (i = 0; i < chain_num; i++) { if (len_vec[i] > xlen) xlen = len_vec[i]; total_len += len_vec[i]; mol_type += mol_vec[i]; } if (!u_opt) Lnorm_ass = total_len / chain_num; u_opt = true; total_len -= xlen; if (total_len > 750) fast_opt = true; /* get all-against-all alignment */ double **TMave_mat; NewArray(&TMave_mat, chain_num, chain_num); vector tmp_str_vec(chain_num, ""); vector > seqxA_mat(chain_num, tmp_str_vec); vector > seqyA_mat(chain_num, tmp_str_vec); for (i = 0; i < chain_num; i++) for (j = 0; j < chain_num; j++) TMave_mat[i][j] = 0; for (i = 0; i < chain_num; i++) { xlen = len_vec[i]; if (xlen < 3) continue; seqx = new char[xlen + 1]; secx = new char[xlen + 1]; NewArray(&xa, xlen, 3); copy_chain_data(a_vec[i], seq_vec[i], sec_vec[i], xlen, xa, seqx, secx); seqxA_mat[i][i] = seqyA_mat[i][i] = (string)(seqx); for (j = i + 1; j < chain_num; j++) { ylen = len_vec[j]; if (ylen < 3) continue; seqy = new char[ylen + 1]; secy = new char[ylen + 1]; NewArray(&ya, ylen, 3); copy_chain_data(a_vec[j], seq_vec[j], sec_vec[j], ylen, ya, seqy, secy); /* declare variable specific to this pair of TMalign */ double t0[3], u0[3][3]; double TM1, TM2; double TM3, TM4, TM5; // for a_opt, u_opt, d_opt double d0_0, TM_0; double d0A, d0B, d0u, d0a; double d0_out = 5.0; string seqM, seqxA, seqyA; // for output alignment double rmsd0 = 0.0; int L_ali; // Aligned length in standard_TMscore double Liden = 0; double TM_ali, rmsd_ali; // TMscore and rmsd in standard_TMscore int n_ali = 0; int n_ali8 = 0; vector do_vec; /* entry function for structure alignment */ if (se_opt) { int *invmap = new int[ylen + 1]; u0[0][0] = u0[1][1] = u0[2][2] = 1; u0[0][1] = u0[0][2] = u0[1][0] = u0[1][2] = u0[2][0] = u0[2][1] = t0[0] = t0[1] = t0[2] = 0; se_main(xa, ya, seqx, seqy, TM1, TM2, TM3, TM4, TM5, d0_0, TM_0, d0A, d0B, d0u, d0a, d0_out, seqM, seqxA, seqyA, do_vec, rmsd0, L_ali, Liden, TM_ali, rmsd_ali, n_ali, n_ali8, xlen, ylen, sequence, Lnorm_ass, d0_scale, 0, false, u_opt, false, mol_type, outfmt_opt, invmap); if (outfmt_opt >= 2) { Liden = L_ali = 0; int r1, r2; for (r2 = 0; r2 < ylen; r2++) { r1 = invmap[r2]; if (r1 < 0) continue; L_ali += 1; Liden += (seqx[r1] == seqy[r2]); } } delete[] invmap; } else TMalign_main(xa, ya, seqx, seqy, secx, secy, t0, u0, TM1, TM2, TM3, TM4, TM5, d0_0, TM_0, d0A, d0B, d0u, d0a, d0_out, seqM, seqxA, seqyA, do_vec, rmsd0, L_ali, Liden, TM_ali, rmsd_ali, n_ali, n_ali8, xlen, ylen, sequence, Lnorm_ass, d0_scale, 0, false, u_opt, false, fast_opt, mol_type, TMcut, 0); /* store result */ TMave_mat[i][j] = TMave_mat[j][i] = TM4; seqxA_mat[i][j] = seqyA_mat[j][i] = seqxA; seqyA_mat[i][j] = seqxA_mat[j][i] = seqyA; // cout< xname_vec; for (i = 0; i < chain_num; i++) xname_vec.push_back( chain_list[i].substr(dir_opt.size()) + chainID_list[i]); vector yname_vec; double *TMave_list; TMave_list = new double[chain_num]; int *assign_list; assign_list = new int[chain_num]; vector msa(ylen, ""); // row is position along msa; column is sequence int compare_num; double TM1_total, TM2_total; double TM3_total, TM4_total, TM5_total; // for a_opt, u_opt, d_opt double d0_0_total, TM_0_total; double d0A_total, d0B_total, d0u_total, d0a_total; double d0_out_total; double rmsd0_total; int L_ali_total; // Aligned length in standard_TMscore double Liden_total; double TM_ali_total, rmsd_ali_total; // TMscore and rmsd in standard_TMscore int n_ali_total; int n_ali8_total; int xlen_total, ylen_total; double TM4_total_max = 0; int max_iter = 5 - (int)(total_len / 200); if (max_iter < 2) max_iter = 2; int iter = 0; vector TM_vec(chain_num, 0); vector d0_vec(chain_num, 0); vector seqID_vec(chain_num, 0); vector > TM_mat(chain_num, TM_vec); vector > d0_mat(chain_num, d0_vec); vector > seqID_mat(chain_num, seqID_vec); for (iter = 0; iter < max_iter; iter++) { /* select representative */ for (j = 0; j < chain_num; j++) TMave_list[j] = 0; for (i = 0; i < chain_num; i++) { for (j = 0; j < chain_num; j++) { // cout<<'\t'< do_vec; /* entry function for structure alignment */ if (se_opt) { int *invmap = new int[ylen + 1]; u0[0][0] = u0[1][1] = u0[2][2] = 1; u0[0][1] = u0[0][2] = u0[1][0] = u0[1][2] = u0[2][0] = u0[2][1] = t0[0] = t0[1] = t0[2] = 0; se_main(xa, ya, seqx, seqy, TM1, TM2, TM3, TM4, TM5, d0_0, TM_0, d0A, d0B, d0u, d0a, d0_out, seqM, seqxA, seqyA, do_vec, rmsd0, L_ali, Liden, TM_ali, rmsd_ali, n_ali, n_ali8, xlen, ylen, sequence, Lnorm_ass, d0_scale, 2, a_opt, u_opt, d_opt, mol_type, outfmt_opt, invmap); if (outfmt_opt >= 2) { Liden = L_ali = 0; int r1, r2; for (r2 = 0; r2 < ylen; r2++) { r1 = invmap[r2]; if (r1 < 0) continue; L_ali += 1; Liden += (seqx[r1] == seqy[r2]); } } delete[] invmap; } else TMalign_main(xa, ya, seqx, seqy, secx, secy, t0, u0, TM1, TM2, TM3, TM4, TM5, d0_0, TM_0, d0A, d0B, d0u, d0a, d0_out, seqM, seqxA, seqyA, do_vec, rmsd0, L_ali, Liden, TM_ali, rmsd_ali, n_ali, n_ali8, xlen, ylen, sequence, Lnorm_ass, d0_scale, 2, a_opt, u_opt, d_opt, fast_opt, mol_type, -1, 0); if (outfmt_opt < 0) output_results( xname_vec[i].c_str(), xname_vec[j].c_str(), "", "", xlen, ylen, t0, u0, TM1, TM2, TM3, TM4, TM5, rmsd0, d0_out, seqM.c_str(), seqxA.c_str(), seqyA.c_str(), Liden, n_ali8, L_ali, TM_ali, rmsd_ali, TM_0, d0_0, d0A, d0B, Lnorm_ass, d0_scale, d0a, d0u, "", 2, // outfmt_opt, ter_opt, false, split_opt, false, "", // o_opt, fname_super+chainID_list1[i], false, a_opt, u_opt, d_opt, false, resi_vec, resi_vec); NewArray(&xt, xlen, 3); do_rotation(xa, xt, xlen, t0, u0); for (r = 0; r < xlen; r++) { a_vec[i][r][0] = xt[r][0]; a_vec[i][r][1] = xt[r][1]; a_vec[i][r][2] = xt[r][2]; } DeleteArray(&xt, xlen); /* clean up */ seqM.clear(); seqxA.clear(); seqyA.clear(); sequence[0].clear(); sequence[1].clear(); delete[] seqx; delete[] secx; DeleteArray(&xa, xlen); delete[] seqy; delete[] secy; DeleteArray(&ya, ylen); do_vec.clear(); } ylen = len_vec[repr_idx]; seqy = new char[ylen + 1]; secy = new char[ylen + 1]; NewArray(&ya, ylen, 3); copy_chain_data(a_vec[repr_idx], seq_vec[repr_idx], sec_vec[repr_idx], ylen, ya, seqy, secy); /* recover alignment */ int ylen_ext = ylen; // chain length double **ya_ext; // structure of single chain char *seqy_ext; // for the protein sequence char *secy_ext; // for the secondary structure for (r = 0; r < msa.size(); r++) msa[r].clear(); msa.clear(); msa.assign(ylen, ""); // row is position along msa; column is sequence vector msa_ext; // row is position along msa; column is sequence for (r = 0; r < ylen; r++) msa[r] = seqy[r]; // for (r=0;r do_vec; se_main(xa, ya, seqx, seqy, TM1, TM2, TM3, TM4, TM5, d0_0, TM_0, d0A, d0B, d0u, d0a, d0_out, seqM, seqxA, seqyA, do_vec, rmsd0, L_ali, Liden, TM_ali, rmsd_ali, n_ali, n_ali8, xlen, ylen, sequence, Lnorm_ass, d0_scale, 0, a_opt, u_opt, d_opt, mol_type, 1, invmap); int rx = 0, ry = 0; ylen_ext = seqxA.size(); NewArray(&ya_ext, ylen_ext, 3); // structure of single chain seqy_ext = new char[ylen_ext + 1]; // for the protein sequence secy_ext = new char[ylen_ext + 1]; // for the secondary structure string tmp_gap = ""; for (r = 0; r < msa[0].size(); r++) tmp_gap += '-'; for (r = msa_ext.size(); r < ylen_ext; r++) msa_ext.push_back(""); // cout<<"x:"<().swap(msa_ext); vector >().swap(TM_pair_vec); for (i = 0; i < chain_num; i++) { tm_idx = assign_list[i]; if (tm_idx < 0) continue; seqyA_mat[i][i] = ""; for (r = 0; r < ylen; r++) seqyA_mat[i][i] += msa[r][tm_idx]; seqxA_mat[i][i] = seqyA_mat[i][i]; // cout< do_vec; se_main(xa, ya, seqx, seqy, TM1, TM2, TM3, TM4, TM5, d0_0, TM_0, d0A, d0B, d0u, d0a, d0_out, seqM, seqxA, seqyA, do_vec, rmsd0, L_ali, Liden, TM_ali, rmsd_ali, n_ali, n_ali8, xlen, ylen, sequence, Lnorm_ass, d0_scale, true, a_opt, u_opt, d_opt, mol_type, 1, invmap); if (xlen <= ylen) { xlen_total += xlen; ylen_total += ylen; TM1_total += TM1; TM2_total += TM2; d0A_total += d0A; d0B_total += d0B; } else { xlen_total += ylen; ylen_total += xlen; TM1_total += TM2; TM2_total += TM1; d0A_total += d0B; d0B_total += d0A; } TM_mat[i][j] = TM2; TM_mat[j][i] = TM1; d0_mat[i][j] = d0B; d0_mat[j][i] = d0A; seqID_mat[i][j] = 1. * Liden / xlen; seqID_mat[j][i] = 1. * Liden / ylen; TM3_total += TM3; TM4_total += TM4; TM5_total += TM5; d0_0_total += d0_0; TM_0_total += TM_0; d0u_total += d0u; d0_out_total += d0_out; rmsd0_total += rmsd0; L_ali_total += L_ali; // Aligned length in standard_TMscore Liden_total += Liden; TM_ali_total += TM_ali; rmsd_ali_total += rmsd_ali; // TMscore and rmsd in standard_TMscore n_ali_total += n_ali; n_ali8_total += n_ali8; /* clean up */ delete[] invmap; seqM.clear(); seqxA.clear(); seqyA.clear(); delete[] seqy; delete[] secy; DeleteArray(&ya, ylen); do_vec.clear(); } delete[] seqx; delete[] secx; DeleteArray(&xa, xlen); } if (TM4_total <= TM4_total_max) break; TM4_total_max = TM4_total; } for (i = 0; i < chain_num; i++) { for (j = 0; j < chain_num; j++) { if (i == j) continue; TM_vec[i] += TM_mat[i][j]; d0_vec[i] += d0_mat[i][j]; seqID_vec[i] += seqID_mat[i][j]; } TM_vec[i] /= (chain_num - 1); d0_vec[i] /= (chain_num - 1); seqID_vec[i] /= (chain_num - 1); } xlen_total /= compare_num; ylen_total /= compare_num; TM1_total /= compare_num; TM2_total /= compare_num; d0A_total /= compare_num; d0B_total /= compare_num; TM3_total /= compare_num; TM4_total /= compare_num; TM5_total /= compare_num; d0_0_total /= compare_num; TM_0_total /= compare_num; d0u_total /= compare_num; d0_out_total /= compare_num; rmsd0_total /= compare_num; L_ali_total /= compare_num; Liden_total /= compare_num; TM_ali_total /= compare_num; rmsd_ali_total /= compare_num; n_ali_total /= compare_num; n_ali8_total /= compare_num; xname = "shorter"; yname = "longer"; string seqM = ""; string seqxA = ""; string seqyA = ""; double t0[3]; double u0[3][3]; stringstream buf; for (i = 0; i < chain_num; i++) { if (assign_list[i] < 0) continue; buf << ">" << xname_vec[i] << "\tL=" << len_vec[i] << "\td0=" << setiosflags(ios::fixed) << setprecision(2) << d0_vec[i] << "\tseqID=" << setiosflags(ios::fixed) << setprecision(3) << seqID_vec[i] << "\tTM-score=" << setiosflags(ios::fixed) << setprecision(5) << TM_vec[i]; if (i == repr_idx) buf << "\t*"; buf << '\n' << seqxA_mat[i][i] << endl; } seqM = buf.str(); seqM = seqM.substr(0, seqM.size() - 1); buf.str(string()); // MergeAlign(seqxA_mat,seqyA_mat,repr_idx,xname_vec,chain_num,seqM); if (outfmt_opt == 0) print_version(); output_mTMalign_results(xname, yname, "", "", xlen_total, ylen_total, t0, u0, TM1_total, TM2_total, TM3_total, TM4_total, TM5_total, rmsd0_total, d0_out_total, seqM.c_str(), seqxA.c_str(), seqyA.c_str(), Liden_total, n_ali8_total, L_ali_total, TM_ali_total, rmsd_ali_total, TM_0_total, d0_0_total, d0A_total, d0B_total, Lnorm_ass, d0_scale, d0a_total, d0u_total, "", outfmt_opt, ter_opt, 0, split_opt, false, "", false, a_opt, u_opt, d_opt, false, resi_vec, resi_vec); if (m_opt || o_opt) { double **ut_mat; // rotation matrices for all-against-all alignment int ui, uj; double t[3], u[3][3]; double rmsd; NewArray(&ut_mat, chain_num, 4 * 3); for (i = 0; i < chain_num; i++) { xlen = ylen = a_vec[i].size(); NewArray(&xa, xlen, 3); NewArray(&ya, xlen, 3); for (r = 0; r < xlen; r++) { xa[r][0] = ua_vec[i][r][0]; xa[r][1] = ua_vec[i][r][1]; xa[r][2] = ua_vec[i][r][2]; ya[r][0] = a_vec[i][r][0]; ya[r][1] = a_vec[i][r][1]; ya[r][2] = a_vec[i][r][2]; } Kabsch(xa, ya, xlen, 1, &rmsd, t, u); for (ui = 0; ui < 3; ui++) for (uj = 0; uj < 3; uj++) ut_mat[i][ui * 3 + uj] = u[ui][uj]; for (uj = 0; uj < 3; uj++) ut_mat[i][9 + uj] = t[uj]; DeleteArray(&xa, xlen); DeleteArray(&ya, xlen); } vector > >().swap(ua_vec); if (m_opt) { assign_list[repr_idx] = -1; output_dock_rotation_matrix(fname_matrix.c_str(), xname_vec, yname_vec, ut_mat, assign_list); } // if (o_opt) output_dock(chain_list, ter_opt, split_opt, // infmt_opt, atom_opt, false, ut_mat, fname_super); if (o_opt) output_mTMalign_pymol(chain_list, infmt_opt, ut_mat, fname_super, o_opt); DeleteArray(&ut_mat, chain_num); } /* clean up */ vector().swap(msa); vector().swap(tmp_str_vec); vector >().swap(seqxA_mat); vector >().swap(seqyA_mat); vector().swap(xname_vec); vector().swap(yname_vec); delete[] TMave_list; DeleteArray(&TMave_mat, chain_num); vector > >().swap(a_vec); // structure of complex vector >().swap(seq_vec); // sequence of complex vector >().swap(sec_vec); // secondary structure of complex vector().swap(mol_vec); // molecule type of complex1, RNA if >0 vector().swap(chainID_list); // list of chainID vector().swap(len_vec); // length of complex vector().swap(TM_vec); vector().swap(d0_vec); vector().swap(seqID_vec); vector >().swap(TM_mat); vector >().swap(d0_mat); vector >().swap(seqID_mat); return 1; } /* sequence order independent alignment */ int SOIalign(string &xname, string &yname, const string &fname_super, const string &fname_lign, const string &fname_matrix, vector &sequence, const double Lnorm_ass, const double d0_scale, const bool m_opt, const int i_opt, const int o_opt, const int a_opt, const bool u_opt, const bool d_opt, const double TMcut, const int infmt1_opt, const int infmt2_opt, const int ter_opt, const int split_opt, const int outfmt_opt, const bool fast_opt, const int cp_opt, const int mirror_opt, const int het_opt, const string &atom_opt, const bool autojustify, const string &mol_opt, const string &dir_opt, const string &dirpair_opt, const string &dir1_opt, const string &dir2_opt, const vector &chain2parse1, const vector &chain2parse2, const vector &model2parse1, const vector &model2parse2, const vector &chain1_list, const vector &chain2_list, const bool se_opt, const int closeK_opt, const int mm_opt) { /* declare previously global variables */ vector > PDB_lines1; // text of chain1 vector > PDB_lines2; // text of chain2 vector mol_vec1; // molecule type of chain1, RNA if >0 vector mol_vec2; // molecule type of chain2, RNA if >0 vector chainID_list1; // list of chainID1 vector chainID_list2; // list of chainID2 int i, j; // file index int chain_i, chain_j; // chain index int r; // residue index int xlen, ylen; // chain length int xchainnum, ychainnum; // number of chains in a PDB file char *seqx, *seqy; // for the protein sequence char *secx, *secy; // for the secondary structure int **secx_bond; // boundary of secondary structure int **secy_bond; // boundary of secondary structure double **xa, **ya; // for input vectors xa[0...xlen-1][0..2] and // ya[0...ylen-1][0..2], in general, // ya is regarded as native structure // --> superpose xa onto ya double **xk, **yk; // k closest residues vector resi_vec1; // residue index for chain1 vector resi_vec2; // residue index for chain2 int read_resi = 0; // whether to read residue index if (o_opt) read_resi = 2; /* loop over file names */ for (i = 0; i < chain1_list.size(); i++) { /* parse chain 1 */ xname = chain1_list[i]; xchainnum = get_PDB_lines(xname, PDB_lines1, chainID_list1, mol_vec1, ter_opt, infmt1_opt, atom_opt, autojustify, split_opt, het_opt, chain2parse1, model2parse1); if (!xchainnum) { cerr << "Warning! Cannot parse file: " << xname << ". Chain number 0." << endl; continue; } for (chain_i = 0; chain_i < xchainnum; chain_i++) { xlen = PDB_lines1[chain_i].size(); if (mol_opt == "RNA") mol_vec1[chain_i] = 1; else if (mol_opt == "protein") mol_vec1[chain_i] = -1; if (!xlen) { cerr << "Warning! Cannot parse file: " << xname << ". Chain length 0." << endl; continue; } else if (xlen < 3) { cerr << "Sequence is too short <3!: " << xname << endl; continue; } NewArray(&xa, xlen, 3); if (closeK_opt >= 3) NewArray(&xk, xlen * closeK_opt, 3); seqx = new char[xlen + 1]; secx = new char[xlen + 1]; xlen = read_PDB(PDB_lines1[chain_i], xa, seqx, resi_vec1, read_resi); if (mirror_opt) for (r = 0; r < xlen; r++) xa[r][2] = -xa[r][2]; if (mol_vec1[chain_i] > 0) make_sec(seqx, xa, xlen, secx, atom_opt); else make_sec(xa, xlen, secx); // secondary structure assignment if (closeK_opt >= 3) getCloseK(xa, xlen, closeK_opt, xk); if (mm_opt == 6) { NewArray(&secx_bond, xlen, 2); assign_sec_bond(secx_bond, secx, xlen); } for (j = (dir_opt.size() > 0) * (i + 1); j < chain2_list.size(); j++) { if (dirpair_opt.size() && i != j) continue; /* parse chain 2 */ if (PDB_lines2.size() == 0) { yname = chain2_list[j]; ychainnum = get_PDB_lines(yname, PDB_lines2, chainID_list2, mol_vec2, ter_opt, infmt2_opt, atom_opt, autojustify, split_opt, het_opt, chain2parse2, model2parse2); if (!ychainnum) { cerr << "Warning! Cannot parse file: " << yname << ". Chain number 0." << endl; continue; } } for (chain_j = 0; chain_j < ychainnum; chain_j++) { ylen = PDB_lines2[chain_j].size(); if (mol_opt == "RNA") mol_vec2[chain_j] = 1; else if (mol_opt == "protein") mol_vec2[chain_j] = -1; if (!ylen) { cerr << "Warning! Cannot parse file: " << yname << ". Chain length 0." << endl; continue; } else if (ylen < 3) { cerr << "Sequence is too short <3!: " << yname << endl; continue; } NewArray(&ya, ylen, 3); if (closeK_opt >= 3) NewArray(&yk, ylen * closeK_opt, 3); seqy = new char[ylen + 1]; secy = new char[ylen + 1]; ylen = read_PDB(PDB_lines2[chain_j], ya, seqy, resi_vec2, read_resi); if (mol_vec2[chain_j] > 0) make_sec(seqy, ya, ylen, secy, atom_opt); else make_sec(ya, ylen, secy); if (closeK_opt >= 3) getCloseK(ya, ylen, closeK_opt, yk); if (mm_opt == 6) { NewArray(&secy_bond, ylen, 2); assign_sec_bond(secy_bond, secy, ylen); } /* declare variable specific to this pair of TMalign */ double t0[3], u0[3][3]; double TM1, TM2; double TM3, TM4, TM5; // for a_opt, u_opt, d_opt double d0_0, TM_0; double d0A, d0B, d0u, d0a; double d0_out = 5.0; string seqM, seqxA, seqyA; // for output alignment double rmsd0 = 0.0; int L_ali; // Aligned length in standard_TMscore double Liden = 0; double TM_ali, rmsd_ali; // TMscore and rmsd in standard_TMscore int n_ali = 0; int n_ali8 = 0; bool force_fast_opt = (getmin(xlen, ylen) > 1500) ? true : fast_opt; int *invmap = new int[ylen + 1]; double *dist_list = new double[ylen + 1]; /* entry function for structure alignment */ if (se_opt) { u0[0][0] = u0[1][1] = u0[2][2] = 1; u0[0][1] = u0[0][2] = u0[1][0] = u0[1][2] = u0[2][0] = u0[2][1] = t0[0] = t0[1] = t0[2] = 0; soi_se_main(xa, ya, seqx, seqy, TM1, TM2, TM3, TM4, TM5, d0_0, TM_0, d0A, d0B, d0u, d0a, d0_out, seqM, seqxA, seqyA, rmsd0, L_ali, Liden, TM_ali, rmsd_ali, n_ali, n_ali8, xlen, ylen, Lnorm_ass, d0_scale, i_opt, a_opt, u_opt, d_opt, mol_vec1[chain_i] + mol_vec2[chain_j], outfmt_opt, invmap, dist_list, secx_bond, secy_bond, mm_opt); if (outfmt_opt >= 2) { Liden = L_ali = 0; int r1, r2; for (r2 = 0; r2 < ylen; r2++) { r1 = invmap[r2]; if (r1 < 0) continue; L_ali += 1; Liden += (seqx[r1] == seqy[r2]); } } } else SOIalign_main(xa, ya, xk, yk, closeK_opt, seqx, seqy, secx, secy, t0, u0, TM1, TM2, TM3, TM4, TM5, d0_0, TM_0, d0A, d0B, d0u, d0a, d0_out, seqM, seqxA, seqyA, invmap, rmsd0, L_ali, Liden, TM_ali, rmsd_ali, n_ali, n_ali8, xlen, ylen, sequence, Lnorm_ass, d0_scale, i_opt, a_opt, u_opt, d_opt, force_fast_opt, mol_vec1[chain_i] + mol_vec2[chain_j], dist_list, secx_bond, secy_bond, mm_opt); /* print result */ if (outfmt_opt == 0) print_version(); output_results( xname.substr(dir1_opt.size() + dir_opt.size() + dirpair_opt.size()), yname.substr(dir2_opt.size() + dir_opt.size() + dirpair_opt.size()), chainID_list1[chain_i], chainID_list2[chain_j], xlen, ylen, t0, u0, TM1, TM2, TM3, TM4, TM5, rmsd0, d0_out, seqM.c_str(), seqxA.c_str(), seqyA.c_str(), Liden, n_ali8, L_ali, TM_ali, rmsd_ali, TM_0, d0_0, d0A, d0B, Lnorm_ass, d0_scale, d0a, d0u, (m_opt ? fname_matrix : "").c_str(), outfmt_opt, ter_opt, false, split_opt, o_opt, fname_super, i_opt, a_opt, u_opt, d_opt, mirror_opt, resi_vec1, resi_vec2); if (outfmt_opt <= 0) { cout << "###############\t###############\t#########" << endl; cout << "#Aligned atom 1\tAligned atom 2 \tDistance#" << endl; int r1, r2; for (r2 = 0; r2 < ylen; r2++) { r1 = invmap[r2]; if (r1 < 0) continue; cout << PDB_lines1[chain_i][r1].substr(12, 15) << '\t' << PDB_lines2[chain_j][r2].substr(12, 15) << '\t' << setw(9) << setiosflags(ios::fixed) << setprecision(3) << dist_list[r2] << '\n'; } cout << "###############\t###############\t#########" << endl; } /* Done! Free memory */ delete[] invmap; delete[] dist_list; seqM.clear(); seqxA.clear(); seqyA.clear(); DeleteArray(&ya, ylen); if (closeK_opt >= 3) DeleteArray(&yk, ylen * closeK_opt); delete[] seqy; delete[] secy; resi_vec2.clear(); if (mm_opt == 6) DeleteArray(&secy_bond, ylen); } // chain_j if (chain2_list.size() > 1) { yname.clear(); for (chain_j = 0; chain_j < ychainnum; chain_j++) PDB_lines2[chain_j].clear(); PDB_lines2.clear(); chainID_list2.clear(); mol_vec2.clear(); } } // j PDB_lines1[chain_i].clear(); DeleteArray(&xa, xlen); if (closeK_opt >= 3) DeleteArray(&xk, xlen * closeK_opt); delete[] seqx; delete[] secx; resi_vec1.clear(); if (mm_opt == 6) DeleteArray(&secx_bond, xlen); } // chain_i xname.clear(); PDB_lines1.clear(); chainID_list1.clear(); mol_vec1.clear(); } // i if (chain2_list.size() == 1) { yname.clear(); for (chain_j = 0; chain_j < ychainnum; chain_j++) PDB_lines2[chain_j].clear(); PDB_lines2.clear(); resi_vec2.clear(); chainID_list2.clear(); mol_vec2.clear(); } return 0; } // Unified engine replacing flexalign_greedy and flexalign_usbcat int flexalign_unified(string &xname, string &yname, const string &fname_super, const string &fname_lign, const string &fname_matrix, vector &sequence, const double Lnorm_ass, const double d0_scale, const bool m_opt, const int i_opt, const int o_opt, const int a_opt, const bool u_opt, const bool d_opt, const double TMcut, const int infmt1_opt, const int infmt2_opt, const int ter_opt, const int split_opt, const int outfmt_opt, const bool fast_opt, const int mirror_opt, const int het_opt, const string &atom_opt, const bool autojustify, const string &mol_opt, const string &dir_opt, const string &dirpair_opt, const string &dir1_opt, const string &dir2_opt, const vector &chain2parse1, const vector &chain2parse2, const vector &model2parse1, const vector &model2parse2, const int byresi_opt, const vector &chain1_list, const vector &chain2_list, const int hinge_opt, const int ss_opt, FlexAlignMode mode = FLEX_STANDARD, bool hinge_set = false, double TMpass = 0.85) { vector > PDB_lines1; vector > PDB_lines2; vector mol_vec1; vector mol_vec2; vector chainID_list1; vector chainID_list2; int i, j, chain_i, chain_j, r, xlen, ylen, xchainnum, ychainnum; char *seqx, *seqy, *secx, *secy; double **xa, **ya; vector resi_vec1; vector resi_vec2; int read_resi = byresi_opt; if (byresi_opt == 0 && o_opt) read_resi = 2; for (i = 0; i < chain1_list.size(); i++) { xname = chain1_list[i]; xchainnum = get_PDB_lines(xname, PDB_lines1, chainID_list1, mol_vec1, ter_opt, infmt1_opt, atom_opt, autojustify, split_opt, het_opt, chain2parse1, model2parse1); if (!xchainnum) { cerr << "Warning! Cannot parse file: " << xname << ". Chain number 0." << endl; continue; } for (chain_i = 0; chain_i < xchainnum; chain_i++) { xlen = PDB_lines1[chain_i].size(); if (mol_opt == "RNA") mol_vec1[chain_i] = 1; else if (mol_opt == "protein") mol_vec1[chain_i] = -1; if (xlen < 3) continue; NewArray(&xa, xlen, 3); seqx = new char[xlen + 1]; secx = new char[xlen + 1]; read_PDB(PDB_lines1[chain_i], xa, seqx, resi_vec1, read_resi); if (mirror_opt) for (r = 0; r < xlen; r++) xa[r][2] = -xa[r][2]; (mol_vec1[chain_i] > 0) ? make_sec(seqx, xa, xlen, secx, atom_opt) : make_sec(xa, xlen, secx); for (j = (dir_opt.size() > 0) * (i + 1); j < chain2_list.size(); j++) { if (dirpair_opt.size() && i != j) continue; if (PDB_lines2.size() == 0) { yname = chain2_list[j]; ychainnum = get_PDB_lines(yname, PDB_lines2, chainID_list2, mol_vec2, ter_opt, infmt2_opt, atom_opt, autojustify, split_opt, het_opt, chain2parse2, model2parse2); if (!ychainnum) continue; } for (chain_j = 0; chain_j < ychainnum; chain_j++) { ylen = PDB_lines2[chain_j].size(); if (mol_opt == "RNA") mol_vec2[chain_j] = 1; else if (mol_opt == "protein") mol_vec2[chain_j] = -1; if (ylen < 3) continue; NewArray(&ya, ylen, 3); seqy = new char[ylen + 1]; secy = new char[ylen + 1]; read_PDB(PDB_lines2[chain_j], ya, seqy, resi_vec2, read_resi); (mol_vec2[chain_j] > 0) ? make_sec(seqy, ya, ylen, secy, atom_opt) : make_sec(ya, ylen, secy); if (byresi_opt) extract_aln_from_resi(sequence, seqx, seqy, resi_vec1, resi_vec2, byresi_opt); // --- CORE DISPATCH LOGIC START --- if (mode == FLEX_USBCAT) { FlexAlignResult usbcat_res; bool force_fast_opt = (getmin(xlen, ylen) > 1500) ? true : fast_opt; usbcat_res.hingeNum = flexalign_usbcat_main( xa, ya, seqx, seqy, secx, secy, usbcat_res.t0, usbcat_res.u0, usbcat_res.tu_vec, usbcat_res.TM1, usbcat_res.TM2, usbcat_res.TM3, usbcat_res.TM4, usbcat_res.TM5, usbcat_res.d0_0, usbcat_res.TM_0, usbcat_res.d0A, usbcat_res.d0B, usbcat_res.d0u, usbcat_res.d0a, usbcat_res.d0_out, usbcat_res.seqM, usbcat_res.seqxA, usbcat_res.seqyA, usbcat_res.do_vec, usbcat_res.rmsd0, usbcat_res.L_ali, usbcat_res.Liden, usbcat_res.TM_ali, usbcat_res.rmsd_ali, usbcat_res.n_ali, usbcat_res.n_ali8, xlen, ylen, sequence, Lnorm_ass, d0_scale, i_opt, a_opt, u_opt, d_opt, force_fast_opt, mol_vec1[chain_i] + mol_vec2[chain_j], hinge_opt, ss_opt, 0, hinge_set, TMpass); if (outfmt_opt == 0) print_version(); output_flexalign_results( xname.substr(dir1_opt.size() + dir_opt.size() + dirpair_opt.size()), yname.substr(dir2_opt.size() + dir_opt.size() + dirpair_opt.size()), chainID_list1[chain_i], chainID_list2[chain_j], xlen, ylen, usbcat_res.t0, usbcat_res.u0, usbcat_res.tu_vec, usbcat_res.TM1, usbcat_res.TM2, usbcat_res.TM3, usbcat_res.TM4, usbcat_res.TM5, usbcat_res.rmsd0, usbcat_res.d0_out, usbcat_res.seqM.c_str(), usbcat_res.seqxA.c_str(), usbcat_res.seqyA.c_str(), usbcat_res.Liden, usbcat_res.n_ali8, usbcat_res.L_ali, usbcat_res.TM_ali, usbcat_res.rmsd_ali, usbcat_res.TM_0, usbcat_res.d0_0, usbcat_res.d0A, usbcat_res.d0B, Lnorm_ass, d0_scale, usbcat_res.d0a, usbcat_res.d0u, (m_opt ? fname_matrix : "").c_str(), outfmt_opt, ter_opt, false, split_opt, o_opt, fname_super, i_opt, a_opt, u_opt, d_opt, mirror_opt, resi_vec1, resi_vec2); } else { // === Standard & Best specific logic === FlexAlignResult best_res; double global_max_TM = -1.0; int start_ss = (mode == FLEX_BEST) ? 0 : ss_opt; int end_ss = (mode == FLEX_BEST) ? 1 : ss_opt; bool force_fast_opt = (getmin(xlen, ylen) > 1500) ? true : fast_opt; for (int cur_ss_opt = start_ss; cur_ss_opt <= end_ss; cur_ss_opt++) { FlexAlignResult cur_res; execute_flexalign_with_fallback( xa, ya, seqx, seqy, secx, secy, xlen, ylen, sequence, Lnorm_ass, d0_scale, i_opt, a_opt, u_opt, d_opt, force_fast_opt, mol_vec1[chain_i] + mol_vec2[chain_j], hinge_opt, cur_ss_opt, cur_res); double cur_max_TM = (cur_res.TM1 > cur_res.TM2) ? cur_res.TM1 : cur_res.TM2; if (cur_max_TM > global_max_TM) { global_max_TM = cur_max_TM; best_res = cur_res; } } if (outfmt_opt == 0) print_version(); output_flexalign_results( xname.substr(dir1_opt.size() + dir_opt.size() + dirpair_opt.size()), yname.substr(dir2_opt.size() + dir_opt.size() + dirpair_opt.size()), chainID_list1[chain_i], chainID_list2[chain_j], xlen, ylen, best_res.t0, best_res.u0, best_res.tu_vec, best_res.TM1, best_res.TM2, best_res.TM3, best_res.TM4, best_res.TM5, best_res.rmsd0, best_res.d0_out, best_res.seqM.c_str(), best_res.seqxA.c_str(), best_res.seqyA.c_str(), best_res.Liden, best_res.n_ali8, best_res.L_ali, best_res.TM_ali, best_res.rmsd_ali, best_res.TM_0, best_res.d0_0, best_res.d0A, best_res.d0B, Lnorm_ass, d0_scale, best_res.d0a, best_res.d0u, (m_opt ? fname_matrix : "").c_str(), outfmt_opt, ter_opt, false, split_opt, o_opt, fname_super, i_opt, a_opt, u_opt, d_opt, mirror_opt, resi_vec1, resi_vec2); } // --- CORE DISPATCH LOGIC END --- // Cleanup memory DeleteArray(&ya, ylen); delete[] seqy; delete[] secy; resi_vec2.clear(); } if (chain2_list.size() > 1) { yname.clear(); for (chain_j = 0; chain_j < ychainnum; chain_j++) PDB_lines2[chain_j].clear(); PDB_lines2.clear(); chainID_list2.clear(); mol_vec2.clear(); } } PDB_lines1[chain_i].clear(); DeleteArray(&xa, xlen); delete[] seqx; delete[] secx; resi_vec1.clear(); } xname.clear(); PDB_lines1.clear(); chainID_list1.clear(); mol_vec1.clear(); } if (chain2_list.size() == 1) { yname.clear(); for (chain_j = 0; chain_j < ychainnum; chain_j++) PDB_lines2[chain_j].clear(); PDB_lines2.clear(); resi_vec2.clear(); chainID_list2.clear(); mol_vec2.clear(); } return 0; } // ======================================================================= // Direct Drop-in Wrappers (No changes needed in main() bindings) // ======================================================================= // int flexalign(string &xname, string &yname, const string &fname_super, const string &fname_lign, const string &fname_matrix, vector &sequence, const double Lnorm_ass, const double d0_scale, const bool m_opt, const int i_opt, const int o_opt, const int a_opt, const bool u_opt, const bool d_opt, const double TMcut, const int infmt1_opt, const int infmt2_opt, const int ter_opt, const int split_opt, const int outfmt_opt, const bool fast_opt, const int mirror_opt, const int het_opt, const string &atom_opt, const bool autojustify, const string &mol_opt, const string &dir_opt, const string &dirpair_opt, const string &dir1_opt, const string &dir2_opt, const vector &chain2parse1, const vector &chain2parse2, const vector &model2parse1, const vector &model2parse2, const int byresi_opt, const vector &chain1_list, const vector &chain2_list, const int hinge_opt, const int ss_opt) // { // return flexalign_unified(xname, yname, fname_super, fname_lign, fname_matrix, sequence, Lnorm_ass, d0_scale, m_opt, i_opt, o_opt, a_opt, u_opt, d_opt, TMcut, infmt1_opt, infmt2_opt, ter_opt, split_opt, outfmt_opt, fast_opt, mirror_opt, het_opt, atom_opt, autojustify, mol_opt, dir_opt, dirpair_opt, dir1_opt, dir2_opt, chain2parse1, chain2parse2, model2parse1, model2parse2, byresi_opt, chain1_list, chain2_list, hinge_opt, ss_opt, FLEX_STANDARD); // } int flexalign_greedy(string &xname, string &yname, const string &fname_super, const string &fname_lign, const string &fname_matrix, vector &sequence, const double Lnorm_ass, const double d0_scale, const bool m_opt, const int i_opt, const int o_opt, const int a_opt, const bool u_opt, const bool d_opt, const double TMcut, const int infmt1_opt, const int infmt2_opt, const int ter_opt, const int split_opt, const int outfmt_opt, const bool fast_opt, const int mirror_opt, const int het_opt, const string &atom_opt, const bool autojustify, const string &mol_opt, const string &dir_opt, const string &dirpair_opt, const string &dir1_opt, const string &dir2_opt, const vector &chain2parse1, const vector &chain2parse2, const vector &model2parse1, const vector &model2parse2, const int byresi_opt, const vector &chain1_list, const vector &chain2_list, const int hinge_opt) { return flexalign_unified(xname, yname, fname_super, fname_lign, fname_matrix, sequence, Lnorm_ass, d0_scale, m_opt, i_opt, o_opt, a_opt, u_opt, d_opt, TMcut, infmt1_opt, infmt2_opt, ter_opt, split_opt, outfmt_opt, fast_opt, mirror_opt, het_opt, atom_opt, autojustify, mol_opt, dir_opt, dirpair_opt, dir1_opt, dir2_opt, chain2parse1, chain2parse2, model2parse1, model2parse2, byresi_opt, chain1_list, chain2_list, hinge_opt, 0 /* ss_opt is ignored in BEST mode */, FLEX_BEST); } int flexalign_usbcat(string &xname, string &yname, const string &fname_super, const string &fname_lign, const string &fname_matrix, vector &sequence, const double Lnorm_ass, const double d0_scale, const bool m_opt, const int i_opt, const int o_opt, const int a_opt, const bool u_opt, const bool d_opt, const double TMcut, const int infmt1_opt, const int infmt2_opt, const int ter_opt, const int split_opt, const int outfmt_opt, const bool fast_opt, const int mirror_opt, const int het_opt, const string &atom_opt, const bool autojustify, const string &mol_opt, const string &dir_opt, const string &dirpair_opt, const string &dir1_opt, const string &dir2_opt, const vector &chain2parse1, const vector &chain2parse2, const vector &model2parse1, const vector &model2parse2, const int byresi_opt, const vector &chain1_list, const vector &chain2_list, const int hinge_opt, bool hinge_set = false, double TMpass = 0.85) { return flexalign_unified(xname, yname, fname_super, fname_lign, fname_matrix, sequence, Lnorm_ass, d0_scale, m_opt, i_opt, o_opt, a_opt, u_opt, d_opt, TMcut, infmt1_opt, infmt2_opt, ter_opt, split_opt, outfmt_opt, fast_opt, mirror_opt, het_opt, atom_opt, autojustify, mol_opt, dir_opt, dirpair_opt, dir1_opt, dir2_opt, chain2parse1, chain2parse2, model2parse1, model2parse2, byresi_opt, chain1_list, chain2_list, hinge_opt, 0 /* ss_opt ignore */, FLEX_USBCAT, hinge_set, TMpass); } int main(int argc, char *argv[]) { if (argc < 2) print_help(); clock_t t1, t2; t1 = clock(); /**********************/ /* get argument */ /**********************/ string xname = ""; string yname = ""; string fname_super = ""; // file name for superposed structure string fname_lign = ""; // file name for user alignment string fname_matrix = ""; // file name for output matrix vector sequence; // get value from alignment file double Lnorm_ass, d0_scale; bool h_opt = false; // print full help message bool v_opt = false; // print version bool m_opt = false; // flag for -m, output rotation matrix int i_opt = 0; // 1 for -i, 3 for -I int o_opt = 0; // 1 for -o, 2 for -rasmol, 3 for -chimerax int a_opt = 0; // flag for -a, do not normalized by average length bool u_opt = false; // flag for -u, normalized by user specified length bool d_opt = false; // flag for -d, user specified d0 bool do_opt = false; // flag for -do, output distance of i-th aligned pair bool full_opt = false; // do not show chain level alignment double TMcut = -1; bool se_opt = false; int infmt1_opt = -1; // PDB or PDBx/mmCIF format for chain_1 int infmt2_opt = -1; // PDB or PDBx/mmCIF format for chain_2 int ter_opt = -1; // default change to 2 (END, or different chainID) int split_opt = -1; // default change to 2 (split each chains) int outfmt_opt = 0; // set -outfmt to full output bool fast_opt = false; // flags for -fast, fTM-align algorithm int cp_opt = 0; // do not check circular permutation int closeK_opt = -1; // number of atoms for SOI initial alignment. // 5 and 0 for -mm 5 and 6 int hinge_opt = 9; // maximum number of hinge allowed for flexible bool hinge_set = false; double TMpass_opt = 0.85; int mirror_opt = 0; // do not align mirror int het_opt = 0; // do not read HETATM residues int mm_opt = 0; // do not perform MM-align bool usbcat_opt = false; // flag for -afp, only valid with -mm 7 string atom_opt = "auto"; // use C alpha atom for protein and C3' for RNA string mol_opt = "auto"; // auto-detect the molecule type as protein/RNA string suffix_opt = ""; // set -suffix to empty string dir_opt = ""; // set -dir to empty string dirpair_opt = ""; // set -dirpair to empty string dir1_opt = ""; // set -dir1 to empty string dir2_opt = ""; // set -dir2 to empty string chainmapfile = ""; // chain mapping between two complexes int byresi_opt = 0; // set -byresi to 0 vector chain1_list; // only when -dir1 is set vector chain2_list; // only when -dir2 is set vector chain2parse1; vector chain2parse2; vector model2parse1; vector model2parse2; vector > chain_pair_list; // only when -dirpair is set for (int i = 1; i < argc; i++) { if (!strcmp(argv[i], "-o")) { if (i >= (argc - 1)) PrintErrorAndQuit("ERROR! Missing value for -o"); if (o_opt == 2) cerr << "Warning! -rasmol is already set. Ignore -o" << endl; else if (o_opt == 3) cerr << "Warning! -chimerax is already set. Ignore -o" << endl; else { fname_super = argv[i + 1]; o_opt = 1; } i++; } else if (!strcmp(argv[i], "-rasmol")) { if (i >= (argc - 1)) PrintErrorAndQuit("ERROR! Missing value for -rasmol"); if (o_opt == 1) cerr << "Warning! -o is already set. Ignore -rasmol" << endl; else if (o_opt == 3) cerr << "Warning! -chimerax is already set. Ignore -rasmol" << endl; else { fname_super = argv[i + 1]; o_opt = 2; } i++; } else if (!strcmp(argv[i], "-chimerax")) { if (i >= (argc - 1)) PrintErrorAndQuit("ERROR! Missing value for -chimerax"); if (o_opt == 1) cerr << "Warning! -o is already set. Ignore -chimerax" << endl; else if (o_opt == 2) cerr << "Warning! -rasmol is already set. Ignore -chimerax" << endl; else { fname_super = argv[i + 1]; o_opt = 3; } i++; } else if (!strcmp(argv[i], "-u") || !strcmp(argv[i], "-L")) { if (i >= (argc - 1)) PrintErrorAndQuit("ERROR! Missing value for -u or -L"); Lnorm_ass = atof(argv[i + 1]); u_opt = true; i++; if (Lnorm_ass <= 0) PrintErrorAndQuit( "ERROR! The value for -u or -L should be >0"); } else if (!strcmp(argv[i], "-a")) { if (i >= (argc - 1)) PrintErrorAndQuit("ERROR! Missing value for -a"); if (!strcmp(argv[i + 1], "T")) a_opt = true; else if (!strcmp(argv[i + 1], "F")) a_opt = false; else { a_opt = atoi(argv[i + 1]); if (a_opt != -2 && a_opt != -1 && a_opt != 1) PrintErrorAndQuit("-a must be -2, -1, 1, T or F"); } i++; } else if (!strcmp(argv[i], "-full")) { if (i >= (argc - 1)) PrintErrorAndQuit("ERROR! Missing value for -full"); if (!strcmp(argv[i + 1], "T")) full_opt = true; else if (!strcmp(argv[i + 1], "F")) full_opt = false; else PrintErrorAndQuit("-full must be T or F"); i++; } else if (!strcmp(argv[i], "-d")) { if (i >= (argc - 1)) PrintErrorAndQuit("ERROR! Missing value for -d"); d0_scale = atof(argv[i + 1]); d_opt = true; i++; } else if (!strcmp(argv[i], "-closeK")) { if (i >= (argc - 1)) PrintErrorAndQuit("ERROR! Missing value for -closeK"); closeK_opt = atoi(argv[i + 1]); i++; } else if (!strcmp(argv[i], "-hinge")) { if (i >= (argc - 1)) PrintErrorAndQuit("ERROR! Missing value for -hinge"); hinge_set = true; hinge_opt = atoi(argv[i + 1]); i++; } else if (!strcmp(argv[i], "-v")) { v_opt = true; } else if (!strcmp(argv[i], "-do")) { do_opt = true; } else if (!strcmp(argv[i], "-h")) { h_opt = true; } else if (!strcmp(argv[i], "-i")) { if (i >= (argc - 1)) PrintErrorAndQuit("ERROR! Missing value for -i"); if (i_opt == 3) PrintErrorAndQuit("ERROR! -i and -I cannot be used together"); fname_lign = argv[i + 1]; i_opt = 1; i++; } else if (!strcmp(argv[i], "-I")) { if (i >= (argc - 1)) PrintErrorAndQuit("ERROR! Missing value for -I"); if (i_opt == 1) PrintErrorAndQuit("ERROR! -I and -i cannot be used together"); fname_lign = argv[i + 1]; i_opt = 3; i++; } else if (!strcmp(argv[i], "-chainmap")) { if (i >= (argc - 1)) PrintErrorAndQuit("ERROR! Missing value for -chainmap"); chainmapfile = argv[i + 1]; i++; } else if (!strcmp(argv[i], "-chain1")) { if (i >= (argc - 1)) PrintErrorAndQuit("ERROR! Missing value for -chain1"); split(argv[i + 1], chain2parse1, ','); i++; } else if (!strcmp(argv[i], "-chain2")) { if (i >= (argc - 1)) PrintErrorAndQuit("ERROR! Missing value for -chain2"); split(argv[i + 1], chain2parse2, ','); i++; } else if (!strcmp(argv[i], "-model1")) { if (i >= (argc - 1)) PrintErrorAndQuit("ERROR! Missing value for -model1"); split(argv[i + 1], model2parse1, ','); i++; } else if (!strcmp(argv[i], "-model2")) { if (i >= (argc - 1)) PrintErrorAndQuit("ERROR! Missing value for -model2"); split(argv[i + 1], model2parse2, ','); i++; } else if (!strcmp(argv[i], "-m")) { if (i >= (argc - 1)) PrintErrorAndQuit("ERROR! Missing value for -m"); fname_matrix = argv[i + 1]; m_opt = true; i++; } // get filename for rotation matrix else if (!strcmp(argv[i], "-fast")) { fast_opt = true; } else if (!strcmp(argv[i], "-se")) { se_opt = true; } else if (!strcmp(argv[i], "-infmt1")) { if (i >= (argc - 1)) PrintErrorAndQuit("ERROR! Missing value for -infmt1"); infmt1_opt = atoi(argv[i + 1]); i++; if (infmt1_opt < -1 || infmt1_opt > 3) PrintErrorAndQuit("ERROR! -infmt1 can only be -1, 0, 1, 2, or 3"); } else if (!strcmp(argv[i], "-infmt2")) { if (i >= (argc - 1)) PrintErrorAndQuit("ERROR! Missing value for -infmt2"); infmt2_opt = atoi(argv[i + 1]); i++; if (infmt2_opt < -1 || infmt2_opt > 3) PrintErrorAndQuit("ERROR! -infmt2 can only be -1, 0, 1, 2, or 3"); } else if (!strcmp(argv[i], "-ter")) { if (i >= (argc - 1)) PrintErrorAndQuit("ERROR! Missing value for -ter"); ter_opt = atoi(argv[i + 1]); i++; } else if (!strcmp(argv[i], "-split")) { if (i >= (argc - 1)) PrintErrorAndQuit("ERROR! Missing value for -split"); split_opt = atoi(argv[i + 1]); i++; } else if (!strcmp(argv[i], "-atom")) { if (i >= (argc - 1)) PrintErrorAndQuit("ERROR! Missing value for -atom"); atom_opt = argv[i + 1]; i++; } else if (!strcmp(argv[i], "-mol")) { if (i >= (argc - 1)) PrintErrorAndQuit("ERROR! Missing value for -mol"); mol_opt = argv[i + 1]; i++; if (mol_opt == "prot") mol_opt = "protein"; else if (mol_opt == "DNA") mol_opt = "RNA"; if (mol_opt != "auto" && mol_opt != "protein" && mol_opt != "RNA") PrintErrorAndQuit("ERROR! Molecule type must be one of the " "following:\nauto, prot (the same as 'protein'), and " "RNA (the same as 'DNA')."); } else if (!strcmp(argv[i], "-dir")) { if (i >= (argc - 1)) PrintErrorAndQuit("ERROR! Missing value for -dir"); dir_opt = argv[i + 1]; i++; } else if (!strcmp(argv[i], "-dirpair")) { if (i >= (argc - 1)) PrintErrorAndQuit("ERROR! Missing value for -dirpair"); dirpair_opt = argv[i + 1]; i++; } else if (!strcmp(argv[i], "-dir1")) { if (i >= (argc - 1)) PrintErrorAndQuit("ERROR! Missing value for -dir1"); dir1_opt = argv[i + 1]; i++; } else if (!strcmp(argv[i], "-dir2")) { if (i >= (argc - 1)) PrintErrorAndQuit("ERROR! Missing value for -dir2"); dir2_opt = argv[i + 1]; i++; } else if (!strcmp(argv[i], "-suffix")) { if (i >= (argc - 1)) PrintErrorAndQuit("ERROR! Missing value for -suffix"); suffix_opt = argv[i + 1]; i++; } else if (!strcmp(argv[i], "-outfmt")) { if (i >= (argc - 1)) PrintErrorAndQuit("ERROR! Missing value for -outfmt"); outfmt_opt = atoi(argv[i + 1]); i++; } else if (!strcmp(argv[i], "-TMcut")) { if (i >= (argc - 1)) PrintErrorAndQuit("ERROR! Missing value for -TMcut"); TMcut = atof(argv[i + 1]); i++; } else if (!strcmp(argv[i], "-byresi") || !strcmp(argv[i], "-tmscore") || !strcmp(argv[i], "-TMscore")) { if (i >= (argc - 1)) PrintErrorAndQuit("ERROR! Missing value for -byresi"); byresi_opt = atoi(argv[i + 1]); i++; } else if (!strcmp(argv[i], "-seq")) { byresi_opt = 5; } else if (!strcmp(argv[i], "-cp")) { mm_opt = 3; } else if (!strcmp(argv[i], "-mirror")) { if (i >= (argc - 1)) PrintErrorAndQuit("ERROR! Missing value for -mirror"); mirror_opt = atoi(argv[i + 1]); i++; } else if (!strcmp(argv[i], "-het")) { if (i >= (argc - 1)) PrintErrorAndQuit("ERROR! Missing value for -het"); het_opt = atoi(argv[i + 1]); i++; if (het_opt != 0 && het_opt != 1 && het_opt != 2) PrintErrorAndQuit("-het must be 0, 1, or 2"); } else if (!strcmp(argv[i], "-mm")) { if (i >= (argc - 1)) PrintErrorAndQuit("ERROR! Missing value for -mm"); mm_opt = atoi(argv[i + 1]); i++; } else if (!strcmp(argv[i], "-afp")) { usbcat_opt = true; } else if (!strcmp(argv[i], "-TMpass")) // Parse the -TMpass argument { if (i >= (argc - 1)) PrintErrorAndQuit("ERROR! Missing value for -TMpass"); TMpass_opt = atof(argv[i + 1]); i++; } else if (xname.size() == 0) xname = argv[i]; else if (yname.size() == 0) yname = argv[i]; else PrintErrorAndQuit(string("ERROR! Undefined option ") + argv[i]); } if (xname.size() == 0 || (yname.size() && dir_opt.size()) || (yname.size() && dirpair_opt.size()) || (yname.size() == 0 && dir_opt.size() == 0 && dirpair_opt.size() == 0)) { if (h_opt) print_help(h_opt); if (v_opt) { print_version(); exit(EXIT_FAILURE); } if (xname.size() == 0) PrintErrorAndQuit("Please provide input structures"); else if (yname.size() == 0 && dir_opt.size() == 0 && dirpair_opt.size() == 0 && mm_opt != 4) PrintErrorAndQuit("Please provide structure B"); else if (yname.size() && dir_opt.size() + dirpair_opt.size()) PrintErrorAndQuit("Please provide only one file name if -dir is set"); } if (suffix_opt.size() && dir_opt.size() + dirpair_opt.size() + dir1_opt.size() + dir2_opt.size() == 0) PrintErrorAndQuit("-suffix is only valid if -dir, -dir1 or -dir2 is set"); if ((dir_opt.size() || dirpair_opt.size() || dir1_opt.size() || dir2_opt.size())) { if (mm_opt != 2 && mm_opt != 4) { if (o_opt) PrintErrorAndQuit("-o cannot be set with -dir, -dir1 or -dir2"); if (m_opt && fname_matrix != "-") PrintErrorAndQuit("-m can only be - or unset when using -dir, -dir1 or -dir2"); } else if ((dir_opt.size() || dirpair_opt.size()) && (dir1_opt.size() || dir2_opt.size())) PrintErrorAndQuit("-dir cannot be set with -dir1 or -dir2"); else if (dir_opt.size() && dirpair_opt.size()) PrintErrorAndQuit("-dir cannot be set with -dirpair"); } if (o_opt && (infmt1_opt != -1 && infmt1_opt != 0 && infmt1_opt != 3)) PrintErrorAndQuit("-o can only be used with -infmt1 -1, 0 or 3"); bool autojustify = (atom_opt == "auto" || atom_opt == "PC4'"); // auto re-pad atom name if (mol_opt == "protein" && atom_opt == "auto") atom_opt = " CA "; else if (mol_opt == "RNA" && atom_opt == "auto") atom_opt = " C3'"; if (atom_opt.size() != 4) { cerr << "ERROR! Atom name must have 4 characters, including space.\n" "For example, C alpha, C3' and P atoms should be specified by\n" "-atom \" CA \", -atom \" P \" and -atom \" C3'\", respectively." << endl; if (atom_opt.size() >= 5 || atom_opt.size() == 0) return 1; else if (atom_opt.size() == 1) atom_opt = " " + atom_opt + " "; else if (atom_opt.size() == 2) atom_opt = " " + atom_opt + " "; else if (atom_opt.size() == 3) atom_opt = " " + atom_opt; cerr << "Change -atom to \"" << atom_opt << "\"" << endl; } if (d_opt && d0_scale <= 0) PrintErrorAndQuit("Wrong value for option -d! It should be >0"); if (outfmt_opt >= 2 && (a_opt || u_opt || d_opt)) PrintErrorAndQuit("-outfmt 2 cannot be used with -a, -u, -L, -d"); if (byresi_opt != 0) { if (i_opt) PrintErrorAndQuit("-TMscore >=1 cannot be used with -i or -I"); if (byresi_opt < 0 || byresi_opt > 7) PrintErrorAndQuit("-TMscore can only be 0 to 7"); if ((byresi_opt == 2 || byresi_opt == 3 || byresi_opt == 6) && ter_opt >= 2) PrintErrorAndQuit("-TMscore 2 and 6 must be used with -ter <=1"); } // if (split_opt==1 && ter_opt!=0) // PrintErrorAndQuit("-split 1 should be used with -ter 0"); // else if (split_opt==2 && ter_opt!=0 && ter_opt!=1) // PrintErrorAndQuit("-split 2 should be used with -ter 0 or 1"); if (split_opt < 0) if (byresi_opt == 2 || byresi_opt == 3) split_opt = 0; else split_opt = 2; else if (split_opt > 2) PrintErrorAndQuit("-split can only be 0, 1 or 2"); if (mm_opt == 3) { cp_opt = true; mm_opt = 0; } if (cp_opt && i_opt) PrintErrorAndQuit("-mm 3 cannot be used with -i or -I"); if (mirror_opt && het_opt != 1) cerr << "WARNING! -mirror was not used with -het 1. " << "D amino acids may not be correctly aligned." << endl; if (ter_opt < 0) { if (mm_opt == 1 || mm_opt == 2 || byresi_opt == 2 || byresi_opt == 3 || byresi_opt == 6 || byresi_opt == 7) ter_opt = 1; else ter_opt = 2; } if (mm_opt) { if (i_opt) PrintErrorAndQuit("-mm cannot be used with -i or -I"); if (u_opt) PrintErrorAndQuit("-mm cannot be used with -u or -L"); // if (cp_opt) PrintErrorAndQuit("-mm cannot be used with -cp"); if (dir_opt.size() && mm_opt == 2) PrintErrorAndQuit("-mm 2 cannot be used with -dir"); if (byresi_opt) PrintErrorAndQuit("-mm cannot be used with -byresi"); if (ter_opt >= 2 && (mm_opt == 1 || mm_opt == 2)) PrintErrorAndQuit("-mm 1 or 2 must be used with -ter 0 or -ter 1"); if (mm_opt == 4 && (yname.size() || dir2_opt.size())) cerr << "WARNING! structure_2 is ignored for -mm 4" << endl; if (dirpair_opt.size() && (mm_opt == 2 || mm_opt == 4)) PrintErrorAndQuit("-mm 2 or 4 cannot be used with -dirpair"); } else if (full_opt) PrintErrorAndQuit("-full can only be used with -mm"); if (o_opt && ter_opt <= 1 && split_opt == 2) { if (mm_opt && o_opt == 2) cerr << "WARNING! -mm may generate incorrect" << " RasMol output due to limitations in PDB file format. " << "When -mm is used, -o is recommended over -rasmol" << endl; else if (mm_opt == 0) cerr << "WARNING! Only the superposition of the" << " last aligned structure pair will be generated" << endl; } if (closeK_opt < 0) { if (mm_opt == 5) closeK_opt = 5; else closeK_opt = 0; } if (mm_opt >= 7 && hinge_opt >= 10) PrintErrorAndQuit("ERROR! -hinge must be <10"); if (usbcat_opt && mm_opt != 7) PrintErrorAndQuit("ERROR! -afp parameter can only be used when -mm 7 is set"); if (chainmapfile.size() && mm_opt != 1) PrintErrorAndQuit("ERROR! -chainmap must be used with -mm 1"); /* read initial alignment file from 'align.txt' */ if (i_opt) read_user_alignment(sequence, fname_lign, i_opt); if (byresi_opt == 6 || byresi_opt == 7) mm_opt = 1; else if (byresi_opt) i_opt = 3; if (m_opt && fname_matrix == "") // Output rotation matrix: matrix.txt PrintErrorAndQuit("ERROR! Please provide a file name for option -m!"); /* parse file list */ int i; if (dirpair_opt.size()) file2chainpairlist(chain1_list, chain2_list, xname, dirpair_opt, suffix_opt); else { if (dir1_opt.size() + dir_opt.size() == 0) chain1_list.push_back(xname); else file2chainlist(chain1_list, xname, dir_opt + dir1_opt, suffix_opt); if (dir_opt.size()) for (i = 0; i < chain1_list.size(); i++) chain2_list.push_back(chain1_list[i]); else if (dir2_opt.size() == 0) chain2_list.push_back(yname); else file2chainlist(chain2_list, yname, dir2_opt, suffix_opt); } if (outfmt_opt == 2) { if (mm_opt == 2) cout << "#Query\tTemplate\tTM" << endl; else if (mm_opt >= 7) cout << "#PDBchain1\tPDBchain2\tTM1\tTM2\t" << "RMSD\tID1\tID2\tIDali\tL1\tL2\tLali\tNblk" << endl; else cout << "#PDBchain1\tPDBchain2\tTM1\tTM2\t" << "RMSD\tID1\tID2\tIDali\tL1\tL2\tLali" << endl; } /* real alignment. entry functions are MMalign_main and * TMalign_main */ if (mm_opt == 0) TMalign(xname, yname, fname_super, fname_lign, fname_matrix, sequence, Lnorm_ass, d0_scale, m_opt, i_opt, o_opt, a_opt, u_opt, d_opt, TMcut, infmt1_opt, infmt2_opt, ter_opt, split_opt, outfmt_opt, fast_opt, cp_opt, mirror_opt, het_opt, atom_opt, autojustify, mol_opt, dir_opt, dirpair_opt, dir1_opt, dir2_opt, chain2parse1, chain2parse2, model2parse1, model2parse2, byresi_opt, chain1_list, chain2_list, se_opt, do_opt); else if (mm_opt == 1) { if (dir_opt.size() > 0 || dir1_opt.size() > 0 || dir2_opt.size() > 0) { for (int ii = 0; ii < chain1_list.size(); ii++) { xname = chain1_list[ii]; vector tmp_vec1(1, xname); for (int jj = 0; jj < chain2_list.size(); jj++) { if (dir_opt.size() > 0 && jj <= ii) continue; yname = chain2_list[jj]; vector tmp_vec2(1, yname); MMalign(xname, yname, fname_super, fname_lign, fname_matrix, sequence, d0_scale, m_opt, o_opt, a_opt, d_opt, full_opt, TMcut, infmt1_opt, infmt2_opt, ter_opt, split_opt, outfmt_opt, fast_opt, mirror_opt, het_opt, atom_opt, autojustify, mol_opt, dir_opt + dir1_opt, dir_opt + dir2_opt, chain2parse1, chain2parse2, model2parse1, model2parse2, tmp_vec1, tmp_vec2, byresi_opt, chainmapfile, se_opt); vector().swap(tmp_vec2); } vector().swap(tmp_vec1); } } else if (dirpair_opt.size() == 0) MMalign(xname, yname, fname_super, fname_lign, fname_matrix, sequence, d0_scale, m_opt, o_opt, a_opt, d_opt, full_opt, TMcut, infmt1_opt, infmt2_opt, ter_opt, split_opt, outfmt_opt, fast_opt, mirror_opt, het_opt, atom_opt, autojustify, mol_opt, dir1_opt, dir2_opt, chain2parse1, chain2parse2, model2parse1, model2parse2, chain1_list, chain2_list, byresi_opt, chainmapfile, se_opt); else { vector tmp_vec1; vector tmp_vec2; for (i = 0; i < chain1_list.size(); i++) { xname = chain1_list[i]; yname = chain2_list[i]; tmp_vec1.push_back(xname); tmp_vec2.push_back(yname); MMalign(xname, yname, fname_super, fname_lign, fname_matrix, sequence, d0_scale, m_opt, o_opt, a_opt, d_opt, full_opt, TMcut, infmt1_opt, infmt2_opt, ter_opt, split_opt, outfmt_opt, fast_opt, mirror_opt, het_opt, atom_opt, autojustify, mol_opt, dirpair_opt, dirpair_opt, chain2parse1, chain2parse2, model2parse1, model2parse2, tmp_vec1, tmp_vec2, byresi_opt, chainmapfile, se_opt); tmp_vec1[0].clear(); tmp_vec1.clear(); tmp_vec2[0].clear(); tmp_vec2.clear(); } } chainmapfile.clear(); } else if (mm_opt == 2) MMdock(xname, yname, fname_super, fname_matrix, sequence, Lnorm_ass, d0_scale, m_opt, o_opt, a_opt, u_opt, d_opt, TMcut, infmt1_opt, infmt2_opt, ter_opt, split_opt, outfmt_opt, fast_opt, mirror_opt, het_opt, atom_opt, autojustify, mol_opt, dir1_opt, dir2_opt, chain2parse1, chain2parse2, model2parse1, model2parse2, chain1_list, chain2_list, do_opt); else if (mm_opt == 3) ; // should be changed to mm_opt=0, cp_opt=true else if (mm_opt == 4) mTMalign(xname, yname, fname_super, fname_matrix, sequence, Lnorm_ass, d0_scale, m_opt, i_opt, o_opt, a_opt, u_opt, d_opt, full_opt, TMcut, infmt1_opt, ter_opt, split_opt, outfmt_opt, fast_opt, het_opt, atom_opt, autojustify, mol_opt, dir_opt, byresi_opt, chain1_list, chain2parse1, model2parse1, se_opt); else if (mm_opt == 5 || mm_opt == 6) SOIalign(xname, yname, fname_super, fname_lign, fname_matrix, sequence, Lnorm_ass, d0_scale, m_opt, i_opt, o_opt, a_opt, u_opt, d_opt, TMcut, infmt1_opt, infmt2_opt, ter_opt, split_opt, outfmt_opt, fast_opt, cp_opt, mirror_opt, het_opt, atom_opt, autojustify, mol_opt, dir_opt, dirpair_opt, dir1_opt, dir2_opt, chain2parse1, chain2parse2, model2parse1, model2parse2, chain1_list, chain2_list, se_opt, closeK_opt, mm_opt); else if (mm_opt == 7) { if (usbcat_opt) flexalign_usbcat(xname, yname, fname_super, fname_lign, fname_matrix, sequence, Lnorm_ass, d0_scale, m_opt, i_opt, o_opt, a_opt, u_opt, d_opt, TMcut, infmt1_opt, infmt2_opt, ter_opt, split_opt, outfmt_opt, fast_opt, mirror_opt, het_opt, atom_opt, autojustify, mol_opt, dir_opt, dirpair_opt, dir1_opt, dir2_opt, chain2parse1, chain2parse2, model2parse1, model2parse2, byresi_opt, chain1_list, chain2_list, hinge_opt, hinge_set, TMpass_opt); else flexalign_greedy(xname, yname, fname_super, fname_lign, fname_matrix, sequence, Lnorm_ass, d0_scale, m_opt, i_opt, o_opt, a_opt, u_opt, d_opt, TMcut, infmt1_opt, infmt2_opt, ter_opt, split_opt, outfmt_opt, fast_opt, mirror_opt, het_opt, atom_opt, autojustify, mol_opt, dir_opt, dirpair_opt, dir1_opt, dir2_opt, chain2parse1, chain2parse2, model2parse1, model2parse2, byresi_opt, chain1_list, chain2_list, hinge_opt); } else cerr << "WARNING! -mm " << mm_opt << " not implemented" << endl; /* clean up */ vector().swap(chain1_list); vector().swap(chain2_list); vector().swap(chain2parse1); vector().swap(chain2parse2); vector().swap(model2parse1); vector().swap(model2parse2); vector().swap(sequence); vector >().swap(chain_pair_list); t2 = clock(); float diff = ((float)t2 - (float)t1) / CLOCKS_PER_SEC; if (outfmt_opt < 2) printf("#Total CPU time is %5.2f seconds\n", diff); return 0; }