QUASAR-Scoring-Matrizen in C++ lesen

English Deutsch

Problem:

Du möchtest Alignments-Matrizen wie BLOSUM62 im QUASAR-Format lesen. Die Lösung muss leicht in C++-Code integrierbar sein.

Lösung

Leider unterstützt das offizielle Toolkit nur Java. Um QUASAR-Matrizen effizient zu lesen, habe ich ein kleines Programm geschrieben, das eine abstrakte AlignmentMatrix-Klasse verwendet, die die Matrix in einem einzelnen Array speichert. Aktuell werden nur quadratische Matrizen mit denselben Zeilen-/Spaltenbeschriftungen unterstützt (das ist die einzige Option, die in realen Bioinformatik-Anwendungen normalerweise Sinn macht).

Die Hauptfunktion testet derzeit nur, ob eine Matrixdatei im QUASAR-Format fehlerfrei gelesen werden kann.

quasar_reader.cpp
/**
 * QUASARReader.cpp
 * Ein effizienter Leser für das QUASAR-Substitutionsmatrix-Format.
 * Siehe http://www.bio.ifi.lmu.de/QUASAR/
 *
 * Benötigt Boost String Algorithms (nur Header)
 *
 * Kompilieren wie folgt: clang++ -o quasar quasar.cpp
 *
 * Version 1.0.3
 *
 * Copyright (c) 2013 Uli Koehler
 * Diese Datei ist als Public Domain veröffentlicht.
 * Bitte erwähne den ursprünglichen Autor, wenn in deiner Software verwendet.
 */
#include <fstream>
#include <iostream>
#include <cassert>
#include <string>
#include <cstring>
#include <map>
#include <sstream>
#include <algorithm>
#include <vector>
#include <boost/algorithm/string.hpp>

/**
 * Eine quadratische Substitutionsmatrix.
 * Diese Klasse unterstützt nur symmetrische Substitutionsmatrizen,
 * die gleiche Zeichensätze mappen.
 */
class SubstitutionMatrix {
public:
    /**
     * Konstruiere eine neue Substitutionsmatrix mit Zeichenindex und Standard-Score.
     */
    SubstitutionMatrix(const std::string& index, double defaultScore = 0) : index(index) {
        size = index.size();
        matrix = new double[size * size];
        for(size_t i = 0; i < size * size; i++) {
            matrix[i] = defaultScore;
        }
    }
    ~SubstitutionMatrix() {
        if(matrix != nullptr) {
            delete[] matrix;
        }
    }
    /**
     * Für ein gegebenes Zeichen, bestimme den Index
     * in der Matrix.
     * @return Den Index, oder string::npos falls nicht gefunden
     */
    int getIndexForChar(char c) {
        return index.find(c);
    }
    /**
     * Setze den Score für das Ersetzen von Zeichen a durch Zeichen b.
     */
    void setScore(char a, char b, double score) {
        setScore(getIndexForChar(a), getIndexForChar(b), score);
    }
    /**
     * Setze einen Score anhand des Index.
     * @param a Der Index von Zeichen a
     * @param b Der Index von Zeichen b
     * @param score Der Score für das Ersetzen von a durch b (oder b durch a)
     */
    void setScore(int a, int b, double score) {
        assert(a >= 0);
        assert(b >= 0);
        assert(a < size);
        assert(b < size);
        //Symmetrische Matrix - a->b == b->a
        matrix[a + size*b] = score;
        matrix[b + size*a] = score;
    }
    /**
     * Hole den Score für das Ersetzen
     * von Zeichen a durch Zeichen b.
     * @return Den Score
     */
    double getScore(char a, char b) {
        //Die Reihenfolge von a und b spielt keine Rolle (symmetrische Matrix)
        return getScore(getIndexForChar(a), getIndexForChar(b));
    }
    /**
     * Hole den Score für das Ersetzen
     * des Zeichens an Position a durch
     * das Zeichen an Position b
     * @return Den Score
     */
    double getScore(int a, int b) {
        assert(a >= 0);
        assert(b >= 0);
        assert(a < size);
        assert(b < size);
        //Die Reihenfolge von a und b spielt keine Rolle (symmetrische Matrix)
        return matrix[b + size * a];
    }
    /**
     * Hole die Matrixgröße.
     * Die Matrix hat getMatrixSize()**2 Einträge.
     */
    size_t getMatrixSize() {
        return size * size;
    }
    /**
     * Hole einen String aller Zeichen im Matrix-Alphabet
     */
    std::string getAlphabet() {
        return index;
    }
    /**
     * Hole den rohen Matrix-Puffer.
     * Die gültige Länge ist sizeof(double) * size * size
     */
    double* getMatrix() {
        return matrix;
    }
private:
    /**
     * Die zweidimensionale Matrix, die die Scores enthält.
     * Zeilen- oder Spalten-Erster-Index spielt keine Rolle,
     * da die Matrix symmetrisch ist.
     */
    double* matrix;
    size_t size;
    /**
     * Die Liste der Alphabet-Zeichen,
     * die den Array-Indizes der Matrix zugeordnet sind.
     * Das erste Zeichen ist Index 0 zugeordnet
     */
    std::string index;
};

SubstitutionMatrix* readSubstitutionMatrix(const char* filename) {
    static const char rowIndexPrefix[] = "ROWINDEX";
    static const char colIndexPrefix[] = "COLINDEX";
    static const char matrixPrefix[] = "MATRIX";
    std::ifstream infile(filename);
    std::string line = "", rowIndex = "", colIndex = "";
    int matrixRowIndex = 0; //Die Anzahl der bereits verarbeiteten MATRIX-Zeilen
    SubstitutionMatrix* matrix = nullptr;
    while(std::getline(infile, line)) {
        if(line.compare(0, sizeof(rowIndexPrefix) - 1, rowIndexPrefix) == 0) {
            rowIndex = line.substr(sizeof(rowIndexPrefix));
            boost::algorithm::trim(rowIndex);
        } else if(line.compare(0, sizeof(colIndexPrefix) - 1, colIndexPrefix) == 0) {
            colIndex = line.substr(sizeof(colIndexPrefix));
            boost::algorithm::trim(colIndex);
        } else if(line.compare(0, sizeof(matrixPrefix) - 1, matrixPrefix) == 0) {
            //Matrix initialisieren, falls noch nicht geschehen
            if(matrix == nullptr) {
                //Matrizen mit unterschiedlichen Zeilen- und Spaltenindizes werden aktuell nicht unterstützt
                assert(!rowIndex.empty());
                assert(!colIndex.empty());
                assert(rowIndex == colIndex);
                matrix = new SubstitutionMatrix(rowIndex);
            }
            //Elemente aufteilen
            int matrixColIndex = 0; //Die Anzahl der bereits verarbeiteten Spalten in der aktuellen Zeile
            std::vector<std::string> elements;
            std::string lineContent = boost::trim_copy(line.substr(sizeof(matrixPrefix)));
            boost::split(elements, lineContent, boost::is_any_of("\t "), boost::token_compress_on);
            for(std::string elem : elements) {
                boost::trim(elem);
                double score = atof(elem.c_str());
                //Zeilen- und Spaltenindex sind austauschbar
                matrix->setScore(matrixRowIndex, matrixColIndex, score);
                matrixColIndex++;
            }
            matrixRowIndex++;
        }
    }
    infile.close();
    return matrix;
}

using std::cout;
using std::cerr;
using std::endl;

/**
 * Die Hauptfunktion testet nur, ob der Matrix-Leser für eine gegebene Matrix nicht fehlschlägt.
 */
int main(int argc, char** argv) {
    if(argc < 2) {
        cerr << "Verwendung: " << argv[0] << " <QUASAR-Eingabedatei>" << endl;
        return 1;
    }
    SubstitutionMatrix* matrix = readSubstitutionMatrix(argv[1]);
    delete matrix;
    return 0;
}

Check out similar posts by category: Bioinformatics, C/C++