[ccp4bb] B-factor statistics around an ion - Summary

[Tim Gruene]

Dear all,

I received a few replies to my request, all but Pavel's privately.
Those included suggestions of using moleman2, pymol, baverage (from ccp4).

Several required manual intervention like visual selection of the
responding atoms, which made them less appealing to me.

I used Robbie Joosten's suggestion who pointed out that the PDB updated
all PDB entries to contain link records and that these link records
match the atom description on the ATOM/HETATM cards.
As my perl script had a bug I could not resolve within an acceptable
period of time I wrote a c++ program which makes use of the STL map to
create per-Atom statistics of the atoms within the link record.

I attach the program in case anyone is interested. It is not debugged
and utterly slow (21s for 50PDB files), but it served my purpose.

Thanks to everyone who responded.

Best,
Tim

On 01/23/2014 01:36 PM, Tim Gruene wrote:
> Dear all,
> 
> could anyone suggest a way to get the average B-factor from a PDB-file
> of those atoms a specific ion binds to (e.g. as judged by header LINK
> records or a distance interval)?
> 
> I would like to get this number for all K-ions from a set of
> PDB-files, and I hope there is a quicker way than using coot to click
> on the surrounding atoms and transfer the numbers into a chart.
> 
> Best,
> Tim
> 
> 

-- 
Dr Tim Gruene
Institut fuer anorganische Chemie
Tammannstr. 4
D-37077 Goettingen

GPG Key ID = A46BEE1A

#include <iostream>
#include <fstream>
#include <vector>
#include <map>
#include <string>

class KData {
public:
    double BK_;
    double enviB_;
    int numenvi_;

    KData() : BK_(0.0),
    enviB_(0.0), numenvi_(0) {
    }
};


void printKs(const std::map<std::string, std::vector<std::string> > & Ks) {
    std::map<std::string, std::vector<std::string> >::const_iterator it;
    for (it = Ks.begin(); it != Ks.end(); ++it) {
        std::cout << "K: " << it->first << ":\n";
        std::vector<std::string>::const_iterator it2;
        for (it2 = it->second.begin(); it2 != it->second.end(); ++it2) {
            std::cout << " ---> " << *it2 << "\n";
        }
    }
}

void printData (const std::map<std::string, KData>& data) {
    std::map<std::string, KData>::const_iterator it;
    for (it = data.begin(); it != data.end(); ++it) {
        std::cout << "# ---> \"" << it->first << "\": "
                << "     number of LINKed atoms: " 
                << it->second.numenvi_ << ";"
                << "  B/   average B:\n" 
                << it->second.BK_ << " "
                << it->second.enviB_ /  it->second.numenvi_<< "\n";
    }
}

void checkline(const std::string& line,
        const std::map<std::string, std::vector<std::string> > & Ks,
        std::map<std::string, KData>& data) {
    std::map<std::string, std::vector<std::string> >::const_iterator it;

    for (it = Ks.begin(); it != Ks.end(); ++it) {
        if (line.find(it->first) != std::string::npos) {
            std::string B(line.substr(60, 6));
            double b(std::stod(B));
            data[it->first].BK_ = b;
        }
        std::vector<std::string>::const_iterator it2;
        for (it2 = it->second.begin(); it2 != it->second.end(); ++it2) {
            if (line.find(*it2) != std::string::npos) {
                std::string B(line.substr(60, 6));
                double b(std::stod(B));
                data[it->first].enviB_ += b;
                ++(data[it->first].numenvi_);
            }
        }
    }
}


int main(int argc, char* argv[]) {
    std::map<std::string, std::vector<std::string> > Kenvi;
    std::map<std::string, KData> data;
    
    if (argc < 2) {
        std::cout << "*** Error: please provide one of more PDB filenames.\n";
    }

    char* line;
    int linewidth(100);
    line = new char[linewidth + 1];

    for (int i = 1; i < argc; ++i) {
        std::ifstream inp(argv[i]);
        if (!inp.is_open()) {
            std::cout << "---> Cannot open file " << argv[i] << "\n";
        }

        // extract link records
        while (true) {
            inp.getline(line, linewidth);
            if (!inp.good()) break;
            std::string myline(line);
            
            // break on CRYST1
            if (myline.substr(0, 6) == "CRYST1") {
                break;
            }
            if (myline.substr(0, 4) == "LINK") {
                std::string atom1 = myline.substr(12, 14);
                std::string atom2 = myline.substr(42, 14);
                if (atom1.find("   K") != std::string::npos) {
                    Kenvi[atom1].push_back(atom2);
                } else if (atom2.find("   K") != std::string::npos) {
                    Kenvi[atom2].push_back(atom1);
                }
            }

        }
        
        // map is set, find atoms
        while (true) {
            inp.getline(line, linewidth);
            if (!inp.good()) break;
            std::string myline(line);
            
            if (myline.substr(0, 6) == "ATOM  " || 
                    myline.substr (0, 6) == "HETATM") {
                // check each line against map content to extract K-values
                checkline (myline, Kenvi, data);
            }
        }
    }
    
    delete [] line;
#ifdef DEBUG
    printKs(Kenvi);
#endif
    printData(data);


    return 0;
}

signature.asc
Description: OpenPGP digital signature