test_sidechain_reconstructor.cc

// Copyright (c) 2013-2020, SIB - Swiss Institute of Bioinformatics and
//                          Biozentrum - University of Basel
// 
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
// 
//   http://www.apache.org/licenses/LICENSE-2.0
// 
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.



#include <promod3/modelling/sidechain_reconstructor.hh>
#include <promod3/core/message.hh>
#define BOOST_TEST_DYN_LINK
#include <boost/test/unit_test.hpp>

#include <ost/io/mol/pdb_reader.hh>
#include <ost/conop/heuristic.hh>
#include <ost/seq/sequence_op.hh>

BOOST_AUTO_TEST_SUITE( sidechain );

using namespace promod3;
using namespace promod3::modelling;

namespace {

ost::mol::EntityHandle LoadTestStructure(const String& pdb_name) {
  ost::io::PDBReader reader(pdb_name, ost::io::IOProfile());
  ost::mol::EntityHandle test_ent = ost::mol::CreateEntity();
  reader.Import(test_ent);
  ost::conop::HeuristicProcessor heu_proc;
  heu_proc.Process(test_ent);

  return test_ent;
}

ost::seq::SequenceList GetSeqList(ost::mol::ChainHandle chain) {
  BOOST_CHECK(chain.IsValid());
  ost::seq::SequenceList seqres_list = ost::seq::CreateSequenceList();
  seqres_list.AddSequence(ost::seq::SequenceFromChain(chain.GetName(), chain));
  return seqres_list;
}

SidechainReconstructorPtr GetScRec(const ost::mol::EntityHandle& test_ent) {
  ost::seq::SequenceList seqres_list = GetSeqList(test_ent.FindChain("A"));
  loop::AllAtomEnv env(seqres_list);
  env.SetInitialEnvironment(test_ent);
  SidechainReconstructorPtr sc_rec(new SidechainReconstructor);
  sc_rec->AttachEnvironment(env);
  return sc_rec;
}

// assume: single chain and res_idx = resnum-1
void CheckRecData(SidechainReconstructionDataPtr sc_data,
                  const std::vector<uint>& start_resnums,
                  const std::vector<uint>& loop_lengths,
                  uint num_residues) {
  // check all_pos size
  BOOST_CHECK_EQUAL(sc_data->all_pos->GetNumResidues(), num_residues);

  // check env_pos consistency
  const uint num_rec_residues = sc_data->env_pos->res_indices.size();
  BOOST_CHECK_EQUAL(sc_data->env_pos->all_pos->GetNumResidues(),
                    num_rec_residues);
  for (uint i = 0; i < num_rec_residues; ++i) {
    const uint res_idx = sc_data->env_pos->res_indices[i];
    BOOST_CHECK(sc_data->env_pos->all_pos->IsAllSet(i));
    BOOST_CHECK(res_idx < num_residues);
    BOOST_CHECK_EQUAL(sc_data->env_pos->all_pos->GetAA(i),
                      sc_data->all_pos->GetAA(res_idx));
  }

  // generate expected data for loops -> set to be ok with flips
  const uint num_loops = start_resnums.size();
  std::set<uint> exp_loop_lengths(loop_lengths.begin(), loop_lengths.end());
  std::set<uint> exp_res_indices;
  for (uint i_loop = 0; i_loop < num_loops; ++i_loop) {
    for (uint idx = 0; idx < loop_lengths[i_loop]; ++idx) {
      exp_res_indices.insert(start_resnums[i_loop] + idx - 1);
    }
  }
  // check loops
  BOOST_CHECK_EQUAL(sc_data->loop_start_indices.size(), num_loops);
  BOOST_CHECK_EQUAL(sc_data->loop_lengths.size(), num_loops);
  std::set<uint> sc_loop_lengths(sc_data->loop_lengths.begin(),
                                 sc_data->loop_lengths.end());
  BOOST_CHECK(sc_loop_lengths == exp_loop_lengths);
  std::set<uint> sc_res_indices;
  for (uint i_loop = 0; i_loop < num_loops; ++i_loop) {
    // all loop res. assumed to be contiguous
    const uint start_idx = sc_data->loop_start_indices[i_loop];
    for (uint idx = 0; idx < sc_data->loop_lengths[i_loop]; ++idx) {
      sc_res_indices.insert(sc_data->env_pos->res_indices[start_idx + idx]);
    }
  }
  BOOST_CHECK(sc_res_indices == exp_res_indices);

  // check N-ter / C-ter info
  BOOST_CHECK_EQUAL(sc_data->env_pos->res_indices.size(),
                    sc_data->is_n_ter.size());
  BOOST_CHECK_EQUAL(sc_data->env_pos->res_indices.size(),
                    sc_data->is_c_ter.size());
  for (uint i = 0; i < sc_data->env_pos->res_indices.size(); ++i) {
    const uint res_idx = sc_data->env_pos->res_indices[i];
    const bool is_nter = (res_idx == 0);
    const bool is_cter = (res_idx == num_residues-1);
    BOOST_CHECK(is_nter == sc_data->is_n_ter[i]);
    BOOST_CHECK(is_cter == sc_data->is_c_ter[i]);
  }

  // check rotamer_res_indices
  for (uint i = 0; i < sc_data->rotamer_res_indices.size(); ++i) {
    const uint i_res = sc_data->rotamer_res_indices[i];
    BOOST_CHECK(i_res < num_rec_residues);
    const uint res_idx = sc_data->env_pos->res_indices[i_res];
    BOOST_CHECK(!(sc_data->all_pos->IsAllSet(res_idx)));
  }

  // check disulfid_bridges
  for (uint i = 0; i < sc_data->disulfid_bridges.size(); ++i) {
    const uint i_res1 = sc_data->disulfid_bridges[i].first;
    const uint i_res2 = sc_data->disulfid_bridges[i].second;
    BOOST_CHECK(i_res1 < num_rec_residues);
    BOOST_CHECK(i_res2 < num_rec_residues);
    const uint res_idx1 = sc_data->env_pos->res_indices[i_res1];
    const uint res_idx2 = sc_data->env_pos->res_indices[i_res2];
    BOOST_CHECK_EQUAL(sc_data->all_pos->GetAA(res_idx1), ost::conop::CYS);
    BOOST_CHECK_EQUAL(sc_data->all_pos->GetAA(res_idx2), ost::conop::CYS);
  }
}

void CheckLoopRec(const SidechainReconstructorPtr sc_rec,
                  uint start_resnum, uint loop_length, uint num_residues) {
  // reconstruct
  SidechainReconstructionDataPtr sc_data;
  sc_data = sc_rec->Reconstruct(start_resnum, loop_length);
  // DEBUG-OUT
  // std::cout << "SRN: " << start_resnum << ", LL: " << loop_length
  //           << ", NR: " << num_residues
  //           << ", RIS: " << sc_data->env_pos->res_indices.size()
  //           << ", RRIS: " << sc_data->rotamer_res_indices.size()
  //           << ", DBS: " << sc_data->disulfid_bridges.size() << "\n";
  // turn into vector for check and to call other one
  std::vector<uint> start_resnums(1, start_resnum);
  std::vector<uint> loop_lengths(1, loop_length);
  // check
  CheckRecData(sc_data, start_resnums, loop_lengths, num_residues);
  // call other one
  std::vector<uint> chain_indices(1, 0);
  sc_data = sc_rec->Reconstruct(start_resnums, loop_lengths, chain_indices);
  // check
  CheckRecData(sc_data, start_resnums, loop_lengths, num_residues);
}

void CheckEmptyLoop(const SidechainReconstructionDataPtr sc_data,
                    uint num_residues) {
  // must all be valid and empty
  BOOST_CHECK_EQUAL(sc_data->env_pos->res_indices.size(), 0u);
  BOOST_CHECK_EQUAL(sc_data->env_pos->all_pos->GetNumResidues(), 0u);
  BOOST_CHECK_EQUAL(sc_data->loop_start_indices.size(), 0u);
  BOOST_CHECK_EQUAL(sc_data->loop_lengths.size(), 0u);
  BOOST_CHECK_EQUAL(sc_data->rotamer_res_indices.size(), 0u);
  BOOST_CHECK_EQUAL(sc_data->disulfid_bridges.size(), 0u);
  BOOST_CHECK_EQUAL(sc_data->is_n_ter.size(), 0u);
  BOOST_CHECK_EQUAL(sc_data->is_c_ter.size(), 0u);
  BOOST_CHECK_EQUAL(sc_data->rigid_frame_indices.size(), 0u);
  BOOST_CHECK_EQUAL(sc_data->all_pos->GetNumResidues(), num_residues);
}

} // anon ns

BOOST_AUTO_TEST_CASE(test_sidechain_reconstructor_loops) {
  // setup
  ost::mol::EntityHandle test_ent = LoadTestStructure("data/1eye_sc_test.pdb");
  // -> missing at least 1 sidechain for each AA (total: 270 res., 30 sc cut)
  // -> sc cut: ['PRO1', 'VAL2', 'GLN3', 'MET5', 'LEU8', 'ASN9', 'THR11',
  //             'ASP12', 'SER14', 'PHE15', 'CYS20', 'TYR21', 'LYS29', 'HIS30',
  //             'ILE41', 'GLU47', 'SER48', 'SER49', 'ARG50', 'PRO51', 'THR54',
  //             'ARG55', 'VAL56', 'ASP57', 'PRO58', 'VAL60', 'GLU61', 'THR62',
  //             'ARG64', 'TRP123']
  uint num_residues = test_ent.GetResidueCount();
  SidechainReconstructorPtr sc_rec = GetScRec(test_ent);
  
  // try to reconstruct some loops
  CheckLoopRec(sc_rec, 1, 4, num_residues);
  CheckLoopRec(sc_rec, 100, 3, num_residues);
  CheckLoopRec(sc_rec, num_residues-4, 5, num_residues);

  // try empty loop
  SidechainReconstructionDataPtr sc_data = sc_rec->Reconstruct(2, 0);
  CheckEmptyLoop(sc_data, num_residues);
  // same with "multi-loop" signature
  std::vector<uint> start_resnums(1, 2);
  std::vector<uint> loop_lengths(1, 0);
  std::vector<uint> chain_indices(1, 0);
  sc_data = sc_rec->Reconstruct(start_resnums, loop_lengths, chain_indices);
  CheckEmptyLoop(sc_data, num_residues);

  // try three loops (note: overlap resolution is tested with IdxHandler)
  loop_lengths[0] = 3;
  start_resnums.push_back(20);
  start_resnums.push_back(30);
  loop_lengths.push_back(2);
  loop_lengths.push_back(1);
  chain_indices.push_back(0);
  chain_indices.push_back(0);
  sc_data = sc_rec->Reconstruct(start_resnums, loop_lengths, chain_indices);
  CheckRecData(sc_data, start_resnums, loop_lengths, num_residues);
}

BOOST_AUTO_TEST_CASE(test_sidechain_reconstructor_loops_crn) {
  // setup crn to test disulfid bridges
  ost::mol::EntityHandle test_ent = LoadTestStructure("data/1crn_sc_test.pdb");
  // -> sc missing for THR2 & CYS3, 3 disulfid bridges expected...
  uint num_residues = test_ent.GetResidueCount();
  SidechainReconstructorPtr sc_rec = GetScRec(test_ent);
  
  // try to reconstruct some loops
  CheckLoopRec(sc_rec, 1, 4, num_residues);
}

BOOST_AUTO_TEST_SUITE_END();