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++
If this post helped you, please consider buying me a coffee or donating via PayPal to support research & publishing of new posts on TechOverflow