color_dec.cpp

   1 /*
   2  * color_dec.cpp
   3  *
   4  *  Created on: October 15, 2009
   5  *      Author: Ben Langmead
   6  */
   7
   8 #include <iostream>
   9 #include <string>
  10 #include <stdlib.h>
  11 #include "alphabet.h"
  12 #include "color_dec.h"
  13 #include "color.h"
  14 #include "qual.h"
  15
  16 using namespace std;
  17
  18 // 4-bit pop count
  19 static int alts[] = {
  20         -1, 1, 1, 2, 1, 2, 2, 3,
  21          1, 2, 2, 3, 2, 3, 3, 4
  22 };
  23 static int firsts[] = {
  24         -1, 0, 1, 0, 2, 0, 1, 0,
  25          3, 0, 1, 0, 2, 0, 1, 0
  26 };
  27
  28 /**
  29  * Given a nucleotide mask, pick one matching nucleotide at random.
  30  */
  31 static int randFromMask(int mask) {
  32         assert_gt(mask, 0);
  33         if(alts[mask] == 1) return firsts[mask];
  34         assert_gt(mask, 0);
  35         assert_lt(mask, 16);
  36         int r = rand() % alts[mask];
  37         assert_geq(r, 0);
  38         assert_lt(r, alts[mask]);
  39         for(int i = 0; i < 4; i++) {
  40                 if((mask & (1 << i)) != 0) {
  41                         if(r == 0) return i;
  42                         r--;
  43                 }
  44         }
  45         cerr << "Shouldn't get here" << endl;
  46         throw 1;
  47         return -1;
  48 }
  49
  50 /**
  51  * Does a 2-bit-encoded base match any bit in a mask?
  52  */
  53 static inline bool matches(int i, int j) {
  54         return ((1 << i) & j) != 0;
  55 }
  56
  57 /**
  58  * Given the dynamic programming table, trace backwards from the last
  59  * column and populate the 's' and 'cmm' strings accordingly.  Whenever
  60  * there are multiple equally good ways of backtracking, choose one at
  61  * random.
  62  */
  63 static void backtrack(int table[4][6][1025], // filled-in table
  64                       const char *read, size_t readi, size_t readf,
  65                       const char *ref, size_t refi, size_t reff,
  66                       char *s,    // final nucleotide string
  67                       char *cmm,  // color mismatches
  68                       char *nmm,  // nucleotide mismatches
  69                       int& cmms,  // # color mismatches
  70                       int& nmms)  // # nucleotide mismatches
  71 {
  72         const size_t len = reff-refi;
  73         cmms = nmms = 0;
  74         int min = INT_MAX;
  75         int bests = 0;
  76         // Determine best base in final column of table
  77         for(int i = 0; i < 4; i++) {
  78                 // Install minimum and backtrack info
  79                 int m = table[i][4][len-1];
  80                 if(m < min) {
  81                         min = m;
  82                         bests = (1 << i);
  83                 } else if(m == min) {
  84                         bests |= (1 << i);
  85                 }
  86         }
  87         // i <- position of rightmost nucleotide
  88         int i = (int)len-1;
  89         // to <- rightmost nucleotide
  90         int to = randFromMask(bests);
  91         while(true) {
  92                 bests = table[to][5][i]; // get next best mask
  93                 s[i--] = to; // install best nucleotide
  94                 if(i < 0) break; // done
  95                 assert_gt(bests, 0);
  96                 assert_lt(bests, 16);
  97                 to = randFromMask(bests); // select
  98         }
  99         // Determine what reference nucleotides were matched against
 100         for(size_t i = 0; i < len; i++) {
 101                 if(matches(s[i], ref[refi+i])) {
 102                         assert_eq(1, alts[(int)ref[refi+i]]);
 103                         // Just plain matched
 104                         nmm[i] = 'M';
 105                 } else {
 106                         // If ref is ambiguous here, does it matter which one we
 107                         // choose?  I don't think so.
 108                         assert_eq(1, alts[(int)ref[refi+i]]);
 109                         // SNP here
 110                         nmm[i] = 'S';
 111                         nmms++;
 112                 }
 113         }
 114         for(size_t i = 0; i < len-1; i++) {
 115                 int c1 = (int)read[readi+i]; // actual
 116                 int c2 = dinuc2color[(int)s[i]][(int)s[i+1]]; // decoded
 117                 assert_leq(c1, 4); assert_geq(c1, 0);
 118                 if(c1 != c2 || c1 == 4) {
 119                         // Actual != decoded
 120                         assert_lt(c2, 4); assert_geq(c2, 0);
 121                         cmm[i] = "0123."[c2];
 122                         cmms++;
 123                 } else {
 124                         cmm[i] = 'M';
 125                 }
 126         }
 127         // done
 128 }
 129
 130 /**
 131  * Decode the colorspace read 'read' as aligned against the reference
 132  * string 'ref', assuming that it's a hit.
 133  */
 134 void decodeHit(
 135                 const char *read, // ASCII colors, '0', '1', '2', '3', '.'
 136                 const char *qual, // ASCII quals, Phred+33 encoded
 137                 size_t readi, // offset of first character within 'read' to consider
 138                 size_t readf, // offset of last char (exclusive) in 'read' to consider
 139                 const char *ref, // reference sequence, as masks
 140                 size_t refi, // offset of first character within 'ref' to consider
 141                 size_t reff, // offset of last char (exclusive) in 'ref' to consider
 142                 int snpPhred, // penalty incurred by a SNP
 143                 char *ns,  // decoded nucleotides are appended here
 144                 char *cmm, // where the color mismatches are in the string
 145                 char *nmm, // where nucleotide mismatches are in the string
 146                 int& cmms, // number of color mismatches
 147                 int& nmms) // number of nucleotide mismatches
 148 {
 149         assert_lt(refi, reff);
 150         assert_lt(readi, readf);
 151         assert_eq(reff-refi-1, readf-readi);
 152
 153         //
 154         // Dynamic programming table; good for colorspace reads up to 1024
 155         // colors in length.
 156         //
 157         int table[4][6][1025];
 158         // 0 -> A, 1 -> C, 2 -> G, 3 -> T, 4 -> min(A, C, G, T),
 159         // 5 -> backtrack mask, 6 -> min mismatches
 160
 161         // The first column of the table just considers the first
 162         // nucleotide and whether it matches the ref nucleotide.
 163         for(int to = 0; to < 4; to++) {
 164                 if(matches(to, ref[refi])) {
 165                         // The assigned subject nucleotide matches the reference;
 166                         // no penalty
 167                         table[to][0][0] = 0;
 168                         table[to][1][0] = 0;
 169                         table[to][2][0] = 0;
 170                         table[to][3][0] = 0;
 171                         table[to][4][0] = 0;
 172                         table[to][5][0] = 15;
 173                 } else {
 174                         // The assigned subject nucleotide does not match the
 175                         // reference nucleotide, so we add a SNP penalty
 176                         table[to][0][0] = snpPhred;
 177                         table[to][1][0] = snpPhred;
 178                         table[to][2][0] = snpPhred;
 179                         table[to][3][0] = snpPhred;
 180                         table[to][4][0] = snpPhred;
 181                         table[to][5][0] = 15;
 182                 }
 183         }
 184
 185         // Successive columns examine successive alignment positions
 186         int omin = INT_MAX, t = 0;
 187         int lastOmin = INT_MAX;
 188         for(size_t c = readi; c < readf; c++) {
 189                 const int readc = (int)read[c];
 190                 assert_leq(readc, 4);
 191                 assert_geq(readc, 0);
 192                 lastOmin = omin;
 193                 omin = INT_MAX;
 194                 // t <- index of column in dynamic programming table
 195                 t = c - readi + 1;
 196                 const int refc = ref[refi + t];
 197                 int from[] = { table[0][4][t-1], table[1][4][t-1],
 198                                table[2][4][t-1], table[3][4][t-1] };
 199                 // For each downstream nucleotide
 200                 for(int to = 0; to < 4; to++) {
 201                         // For each upstream nucleotide
 202                         int min = INT_MAX;
 203                         const int goodfrom = nuccol2nuc[to][readc];
 204                         int q = qual[c];
 205                         // Reward the preceding position
 206                         if(goodfrom < 4) from[goodfrom] -= q;
 207                         min = from[0];
 208                         table[to][5][t] = 1;
 209                         if(from[1] < min) {
 210                                 min = from[1];
 211                                 table[to][5][t] = 2;
 212                         } else if(from[1] == min) {
 213                                 table[to][5][t] |= 2;
 214                         }
 215                         if(from[2] < min) {
 216                                 min = from[2];
 217                                 table[to][5][t] = 4;
 218                         } else if(from[2] == min) {
 219                                 table[to][5][t] |= 4;
 220                         }
 221                         if(from[3] < min) {
 222                                 min = from[3];
 223                                 table[to][5][t] = 8;
 224                         } else if(from[3] == min) {
 225                                 table[to][5][t] |= 8;
 226                         }
 227                         min += q;
 228                         if(!matches(to, refc)) {
 229                                 min += snpPhred;
 230                         }
 231                         table[to][4][t] = min;
 232                         if(min < omin) omin = min;
 233                         if(goodfrom < 4) from[goodfrom] += q;
 234                 }
 235         }
 236
 237         t++;
 238         assert_eq(t, (int)(reff - refi));
 239         // Install the best backward path into ns, cmm, nmm
 240         backtrack(table,
 241                   read, readi, readi + t - 1,
 242                   ref, refi, refi + t,
 243                   ns, cmm, nmm, cmms, nmms);
 244 }
 245
 246 #ifdef MAIN_COLOR_DEC
 247
 248 #include <sstream>
 249 #include <getopt.h>
 250
 251 static const char *short_opts = "s:m:r:e:";
 252 static struct option long_opts[] = {
 253         {(char*)"snppen",  required_argument, 0, 's'},
 254         {(char*)"misspen", required_argument, 0, 'm'},
 255         {(char*)"seed",    required_argument, 0, 'r'},
 256         {(char*)"maxpen",  required_argument, 0, 'e'}
 257 };
 258
 259 template<typename T>
 260 T parse(const char *s) {
 261         T tmp;
 262         stringstream ss(s);
 263         ss >> tmp;
 264         return tmp;
 265 }
 266
 267 int main(int argc, char **argv) {
 268         int option_index = 0;
 269         int next_option;
 270         int snppen = 30;
 271         int misspen = 20;
 272         int maxPenalty = 70;
 273         unsigned seed = 0;
 274         do {
 275                 next_option = getopt_long(argc, argv, short_opts, long_opts, &option_index);
 276                 switch (next_option) {
 277                         case 's': snppen = parse<int>(optarg); break;
 278                         case 'm': misspen = parse<int>(optarg); break;
 279                         case 'r': seed = parse<unsigned>(optarg); break;
 280                         case 'e': maxPenalty = parse<int>(optarg); break;
 281                         case -1: break;
 282                         default: {
 283                                 cerr << "Unknown option: " << (char)next_option << endl;
 284                                 exit(1);
 285                         }
 286                 }
 287         } while(next_option != -1);
 288         srand(seed);
 289         if(argc - optind < 2) {
 290                 cerr << "Not enough options" << endl;
 291                 exit(1);
 292         }
 293         string read, ref;
 294         read = argv[optind];
 295         for(size_t i = 0; i < read.length(); i++) {
 296                 read[i] = asc2col[(int)read[i]];
 297                 assert_leq(read[i], 4);
 298                 assert_geq(read[i], 0);
 299         }
 300         ref = argv[optind+1];
 301         for(size_t i = 0; i < ref.length(); i++) {
 302                 int num = 0;
 303                 int alts[] = {4, 4, 4, 4};
 304                 decodeNuc(toupper(ref[i]), num, alts);
 305                 assert_leq(num, 4);
 306                 assert_gt(num, 0);
 307                 ref[i] = 0;
 308                 for(int j = 0; j < num; j++) {
 309                         ref[i] |= (1 << alts[j]);
 310                 }
 311         }
 312         string ns;
 313         string quals;
 314         quals.resize(read.length(), misspen);
 315         string cmm, nmm;
 316         int score = decode(read, quals, 0, read.length(),
 317                            ref, 0, ref.length(), maxPenalty,
 318                            snppen, ns, cmm, nmm);
 319         cout << " Score: " << score << " (max: " << maxPenalty << ")" << endl;
 320         cout << "   MMs:  ";
 321         for(size_t i = 0; i < cmm.length(); i++) {
 322                 cout << cmm[i] << " ";
 323         }
 324         cout << endl;
 325         cout << "Colors:  ";
 326         for(size_t i = 0; i < read.length(); i++) {
 327                 printColor((int)read[i]);
 328                 cout << " ";
 329         }
 330         cout << endl;
 331         cout << " Bases: ";
 332         for(size_t i = 0; i < ns.length(); i++) {
 333                 cout << "ACGTN"[(int)ns[i]] << " ";
 334         }
 335         cout << endl;
 336         cout << "   Ref: ";
 337         for(size_t i = 0; i < ref.length(); i++) {
 338                 cout << mask2iupac[(int)ref[i]] << " ";
 339         }
 340         cout << endl;
 341         cout << "   MMs: ";
 342         for(size_t i = 0; i < ref.length(); i++) {
 343                 cout << nmm[i] << " ";
 344         }
 345         cout << endl;
 346 }
 347 #endif