From 9a6ea8e71c68b8b4f6d82b919527015845f874cc Mon Sep 17 00:00:00 2001 From: Marco Biasini <marco.biasini@unibas.ch> Date: Wed, 17 Nov 2010 09:38:09 +0100 Subject: [PATCH] added executable to calculate LDT score --- modules/mol/alg/src/CMakeLists.txt | 3 + modules/mol/alg/src/ldt.cc | 94 ++++++++++++++++++++++++++++++ 2 files changed, 97 insertions(+) create mode 100644 modules/mol/alg/src/ldt.cc diff --git a/modules/mol/alg/src/CMakeLists.txt b/modules/mol/alg/src/CMakeLists.txt index e613a9c73..fef02349e 100644 --- a/modules/mol/alg/src/CMakeLists.txt +++ b/modules/mol/alg/src/CMakeLists.txt @@ -29,6 +29,9 @@ if (ENABLE_IMG) set(MOL_ALG_DEPS ${MOL_ALG_DEPS} img img_alg) endif() +executable(NAME ldt SOURCES ldt.cc + DEPENDS_ON io mol_alg STATIC) + module(NAME mol_alg SOURCES ${OST_MOL_ALG_SOURCES} HEADERS ${OST_MOL_ALG_HEADERS} HEADER_OUTPUT_DIR ost/mol/alg diff --git a/modules/mol/alg/src/ldt.cc b/modules/mol/alg/src/ldt.cc new file mode 100644 index 000000000..589bd125f --- /dev/null +++ b/modules/mol/alg/src/ldt.cc @@ -0,0 +1,94 @@ +#include <ost/mol/alg/local_dist_test.hh> +#include <ost/io/mol/pdb_reader.hh> +#include <ost/io/io_exception.hh> +using namespace ost; +using namespace ost::io; +using namespace ost::mol; + +EntityHandle load(const String& file) +{ + try { + PDBReader reader(file); + if (reader.HasNext()) { + EntityHandle ent=CreateEntity(); + reader.Import(ent); + return ent; + } + std::cerr << "ERROR: '" << file << "' does not contain any ATOM records. " + << "Are you sure this is a PDB file?" << std::endl; + return EntityHandle(); + } catch (io::IOException& e) { + std::cerr << "ERROR: " << e.what() << std::endl; + return EntityHandle(); + } +} + +void usage() +{ + std::cerr << "usage: ldt [options] <mod1> [mod1 [mod2]] <ref>" << std::endl; + std::cerr << " -s selection performed on ref" << std::endl; + std::cerr << " -c use Calphas only" << std::endl; + std::cerr << " -f filter clashes (not implemented yet)" << std::endl; + std::cerr << " -t fault tolerant parsing" << std::endl; +} + +int main (int argc, char* const *argv) +{ + // parse options + String sel; + int flags=0; + bool filter_clashes=false; + char ch=0; + while((ch=getopt(argc, argv, "ftcs:"))!=-1) { + switch (ch) { + case 'c': + flags|=PDB::CALPHA_ONLY; + break; + case 't': + flags|=PDB::SKIP_FAULTY_RECORDS; + break; + case 'f': + filter_clashes=true; + break; + case 's': + sel=optarg; + break; + case '?': + default: + usage(); + return -1; + } + } + argc-=optind; + argv+=optind; + if (argc<2) { + usage(); + return -1; + } + PDB::PushFlags(flags); + + String ref_file=argv[argc-1]; + EntityHandle ref=load(ref_file); + if (!ref) { + return -1; + } + EntityView ref_view=ref.Select(sel); + for (int i=0; i<argc-1; ++i) { + EntityHandle model=load(argv[i]); + if (!model) { + if (!(flags&PDB::SKIP_FAULTY_RECORDS)) { + return -1; + } + continue; + } + EntityView v=model.CreateFullView(); + float cutoffs[]={0.5,1,2,4}; + float ldt=0.0; + for (int n=0; n<4; ++n) { + ldt+=alg::LocalDistTest(v, ref_view, cutoffs[n], 8.0); + } + ldt/=4.0; + std::cout << argv[i] << " " << ldt << std::endl; + } + return 0; +} \ No newline at end of file -- GitLab