~ubuntu-branches/ubuntu/wily/openms/wily

defaults_.setValue("charge_remote_threshold", 0.3, "If the probability for the proton at the N-terminus is lower than this value, enable charge-remote pathways");

defaults_.setValue("charge_directed_threshold", 0.3, "Limit the proton availability at the N-terminus to at least this value for charge-directed pathways");

defaults_.setValue("model_depth", 10, "The number of explicitly modeled backbone cleavages from N-terminus and C-terminus, would be 9 for the default value", StringList::create("advanced"));

defaults_.setValue("visible_model_depth", 50, "The maximal possible size of a peptide to be modeled", StringList::create("advanced"));

// tolerances

defaults_.setValue("precursor_mass_tolerance", 3.0, "Mass tolerance of the precursor peak, used to identify the precursor peak and its loss peaks for training");

defaults_.setValue("fragment_mass_tolerance", 0.3, "Peak mass tolerance of the product ions, used to identify the ions for training");

// modification parameters

vector<String> all_mods;

ModificationsDB::getInstance()->getAllSearchModifications(all_mods);

defaults_.setValue("variable_modifications", StringList::create("Oxidation (M),Carbamyl (K)"), "Modifications which should be included in the model, represented by PSI-MOD accessions.");

defaults_.setValidStrings("variable_modifications", all_mods);

defaults_.setValue("fixed_modifications", StringList::create(""), "Modifications which should replace the unmodified amino acid, represented by PSI-MOD accessions.");

defaults_.setValidStrings("fixed_modifications", all_mods);

defaults_.setValue("min_enhancement_factor", 2.0, "Minimal factor for bxyz backbone cleavages.", StringList::create("advanced"));

defaults_.setValue("min_y_ion_intensity", 0.0, "Minimal intensity for y ions.", StringList::create("advanced"));

defaults_.setValue("min_b_ion_intensity", 0.0, "Minimal intensity for b ions.", StringList::create("advanced"));

100

defaults_.setValue("min_a_ion_intensity", 0.0, "Minimal intensity for a ions.", StringList::create("advanced"));

101

defaults_.setValue("min_y_loss_intensity", 0.0, "Minimal intensity for y ions with neutral loss.", StringList::create("advanced"));

102

defaults_.setValue("min_b_loss_intensity", 0.0, "Minimal intensity for b ions with neutral loss.", StringList::create("advanced"));

103

104

defaults_.setValue("side_chain_activation", 0.2, "Additional activation of proton sitting at side chain, especially important at lysine and histidine residues");

105

defaults_.setValue("pseudo_counts", 1e-15, "Value which is added for every transition trained of the underlying hidden Markov model", StringList::create("advanced"));

106

defaults_.setValue("max_isotope", 2, "Maximal isotope peak which is reported by the model, 0 means all isotope peaks are reported");

107

108

defaults_.setValue("max_fragment_charge_training", 2, "Maximal allowed charge states for ions to be considered for training");

109

defaults_.setValue("max_fragment_charge", 4, "Maximal charge state allowed for fragment ions");

110

111

defaultsToParam_();

112

}

113

114

PILISModel::~PILISModel()

115

{

116

}

117

118

PILISModel::PILISModel(const PILISModel& model)

119

: DefaultParamHandler(model),

120

hmm_(model.hmm_),

121

prot_dist_(model.prot_dist_),

122

tsg_(model.tsg_),

123

valid_(model.valid_),

124

peaks_(model.peaks_),

125

spectra_aligner_(model.spectra_aligner_),

126

precursor_model_cr_(model.precursor_model_cr_),

127

precursor_model_cd_(model.precursor_model_cd_),

128

129

a_ion_losses_cr_(model.a_ion_losses_cr_),

130

a_ion_losses_cd_(model.a_ion_losses_cd_),

131

132

b_ion_losses_cr_(model.b_ion_losses_cr_),

133

b_ion_losses_cd_(model.b_ion_losses_cd_),

134

135

b2_ion_losses_cr_(model.b2_ion_losses_cr_),

136

b2_ion_losses_cd_(model.b2_ion_losses_cd_),

137

138

y_ion_losses_cr_(model.y_ion_losses_cr_),

139

y_ion_losses_cd_(model.y_ion_losses_cd_)

140

{

141

}

142

143

PILISModel& PILISModel::operator = (const PILISModel& model)

144

{

145

if (this != &model)

146

{

147

DefaultParamHandler::operator=(model);

148

hmm_ = model.hmm_;

149

prot_dist_ = model.prot_dist_;

150

tsg_ = model.tsg_;

151

valid_ = model.valid_;

152

peaks_ = model.peaks_;

153

spectra_aligner_ = model.spectra_aligner_;

154

precursor_model_cr_ = model.precursor_model_cr_;

155

precursor_model_cd_ = model.precursor_model_cd_;

156

157

a_ion_losses_cr_ = model.a_ion_losses_cr_;

158

a_ion_losses_cd_ = model.a_ion_losses_cd_;

159

160

b_ion_losses_cr_ = model.b_ion_losses_cr_;

161

b_ion_losses_cd_ = model.b_ion_losses_cd_;

162

163

b2_ion_losses_cr_ = model.b2_ion_losses_cr_;

164

b2_ion_losses_cd_ = model.b2_ion_losses_cd_;

165

166

y_ion_losses_cr_ = model.y_ion_losses_cr_;

167

y_ion_losses_cd_ = model.y_ion_losses_cd_;

168

}

169

return *this;

170

}

171

172

void PILISModel::init(bool generate_models)

173

{

174

if (generate_models)

175

{

176

PILISModelGenerator gen;

177

Param gen_param(gen.getParameters());

178

gen_param.setValue("variable_modifications", (StringList)param_.getValue("variable_modifications"));

179

gen_param.setValue("fixed_modifications", (StringList)param_.getValue("fixed_modifications"));

180

gen_param.setValue("model_depth", param_.getValue("model_depth"));

181

gen_param.setValue("visible_model_depth", param_.getValue("visible_model_depth"));

182

gen.setParameters(gen_param);

183

gen.getModel(hmm_);

184

}

185

186

187

Param pre_param(precursor_model_cr_.getParameters());

188

pre_param.setValue("fragment_mass_tolerance", (DoubleReal)param_.getValue("fragment_mass_tolerance"));

189

pre_param.setValue("variable_modifications", (StringList)param_.getValue("variable_modifications"));

190

pre_param.setValue("fixed_modifications", (StringList)param_.getValue("fixed_modifications"));

191

pre_param.setValue("pseudo_counts", (DoubleReal)param_.getValue("pseudo_counts"));

192

pre_param.setValue("C_term_H2O_loss", "true");

193

pre_param.setValue("ion_name", "p");

194

precursor_model_cr_.setParameters(pre_param);

195

precursor_model_cd_.setParameters(pre_param);

196

197

if (generate_models)

198

{

199

precursor_model_cr_.generateModel();

200

precursor_model_cd_.generateModel();

201

}

202

203

pre_param.setValue("C_term_H2O_loss", "false");

204

pre_param.setValue("enable_double_losses", "false");

205

pre_param.setValue("ion_name", "b");

206

b_ion_losses_cr_.setParameters(pre_param);

207

b_ion_losses_cd_.setParameters(pre_param);

208

if (generate_models)

209

{

210

b_ion_losses_cr_.generateModel();

211

b_ion_losses_cd_.generateModel();

212

}

213

214

pre_param.setValue("C_term_H2O_loss", "false");

215

pre_param.setValue("enable_double_losses", "false");

216

pre_param.setValue("ion_name", "b2");

217

b2_ion_losses_cr_.setParameters(pre_param);

218

b2_ion_losses_cd_.setParameters(pre_param);

219

if (generate_models)

220

{

221

b2_ion_losses_cr_.generateModel();

222

b2_ion_losses_cd_.generateModel();

223

}

224

225

pre_param.setValue("C_term_H2O_loss", "false");

226

pre_param.setValue("enable_double_losses", "false");

227

pre_param.setValue("ion_name", "a");

228

a_ion_losses_cr_.setParameters(pre_param);

229

a_ion_losses_cd_.setParameters(pre_param);

230

if (generate_models)

231

{

232

a_ion_losses_cr_.generateModel();

233

a_ion_losses_cd_.generateModel();

234

}

235

236

pre_param.setValue("C_term_H2O_loss", "true");

237

pre_param.setValue("ion_name", "y");

238

y_ion_losses_cr_.setParameters(pre_param);

239

y_ion_losses_cd_.setParameters(pre_param);

240

if (generate_models)

241

{

242

y_ion_losses_cr_.generateModel();

243

y_ion_losses_cd_.generateModel();

244

}

245

246

if (generate_models)

247

{

248

valid_ = true;

249

}

250

}

251

252

253

void PILISModel::readFromFile(const String& filename)

254

{

255

// read the model

256

if (!File::readable(filename))

257

{

258

throw Exception::FileNotReadable(__FILE__, __LINE__, __PRETTY_FUNCTION__, filename);

259

}

260

if (File::empty(filename))

261

{

262

throw Exception::FileEmpty(__FILE__, __LINE__, __PRETTY_FUNCTION__, filename);

263

}

264

265

init(false);

266

267

TextFile file;

268

file.load(filename, true);

269

270

TextFile::Iterator it_begin(file.begin()), it_end(file.begin());

271

it_begin = file.search(it_begin, "BASE_MODEL_BEGIN");

272

it_end = file.search(it_begin, "BASE_MODEL_END");

273

parseHMMModel_(++it_begin, it_end, hmm_, param_);

274

275

// seek to next interval

276

it_begin = file.search(it_end, "PRECURSOR_MODEL_CR_BEGIN");

277

it_end = file.search(it_begin, "PRECURSOR_MODEL_CR_END");

278

Param p;

279

HiddenMarkovModel precursor_model_cr_hmm;

280

parseHMMModel_(++it_begin, it_end, precursor_model_cr_hmm, p);

281

precursor_model_cr_.setHMM(precursor_model_cr_hmm);

282

precursor_model_cr_.setParameters(p);

283

284

it_begin = file.search(it_end, "PRECURSOR_MODEL_CD_BEGIN");

285

it_end = file.search(it_begin, "PRECURSOR_MODEL_CD_END");

286

p.clear();

287

HiddenMarkovModel precursor_model_cd_hmm;

288

parseHMMModel_(++it_begin, it_end, precursor_model_cd_hmm, p);

289

precursor_model_cd_.setHMM(precursor_model_cd_hmm);

290

precursor_model_cd_.setParameters(p);

291

292

293

// b loss models

294

it_begin = file.search(it_end, "BION_LOSS_MODEL_CR_BEGIN");

295

it_end = file.search(it_begin, "BION_LOSS_MODEL_CR_END");

296

p.clear();

297

HiddenMarkovModel b_ion_loss_hmm_cr;

298

parseHMMModel_(++it_begin, it_end, b_ion_loss_hmm_cr, p);

299

b_ion_losses_cr_.setHMM(b_ion_loss_hmm_cr);

300

b_ion_losses_cr_.setParameters(p);

301

302

it_begin = file.search(it_end, "BION_LOSS_MODEL_CD_BEGIN");

303

it_end = file.search(it_begin, "BION_LOSS_MODEL_CD_END");

304

p.clear();

305

HiddenMarkovModel b_ion_loss_hmm_cd;

306

parseHMMModel_(++it_begin, it_end, b_ion_loss_hmm_cd, p);

307

b_ion_losses_cd_.setHMM(b_ion_loss_hmm_cd);

308

b_ion_losses_cd_.setParameters(p);

309

310

it_begin = file.search(it_end, "B2ION_LOSS_MODEL_CR_BEGIN");

311

it_end = file.search(it_begin, "B2ION_LOSS_MODEL_CR_END");

312

p.clear();

313

HiddenMarkovModel b2_ion_loss_hmm_cr;

314

parseHMMModel_(++it_begin, it_end, b2_ion_loss_hmm_cr, p);

315

b2_ion_losses_cr_.setHMM(b2_ion_loss_hmm_cr);

316

b2_ion_losses_cr_.setParameters(p);

317

318

it_begin = file.search(it_end, "B2ION_LOSS_MODEL_CD_BEGIN");

319

it_end = file.search(it_begin, "B2ION_LOSS_MODEL_CD_END");

320

p.clear();

321

HiddenMarkovModel b2_ion_loss_hmm_cd;

322

parseHMMModel_(++it_begin, it_end, b2_ion_loss_hmm_cd, p);

323

b2_ion_losses_cd_.setHMM(b2_ion_loss_hmm_cd);

324

b2_ion_losses_cd_.setParameters(p);

325

326

// a-ion loss model

327

it_begin = file.search(it_end, "AION_LOSS_MODEL_CR_BEGIN");

328

it_end = file.search(it_begin, "AION_LOSS_MODEL_CR_END");

329

p.clear();

330

HiddenMarkovModel a_ion_loss_hmm_cr;

331

parseHMMModel_(++it_begin, it_end, a_ion_loss_hmm_cr, p);

332

a_ion_losses_cr_.setHMM(a_ion_loss_hmm_cr);

333

a_ion_losses_cr_.setParameters(p);

334

335

it_begin = file.search(it_end, "AION_LOSS_MODEL_CD_BEGIN");

336

it_end = file.search(it_begin, "AION_LOSS_MODEL_CD_END");

337

p.clear();

338

HiddenMarkovModel a_ion_loss_hmm_cd;

339

parseHMMModel_(++it_begin, it_end, a_ion_loss_hmm_cd, p);

340

a_ion_losses_cd_.setHMM(a_ion_loss_hmm_cd);

341

a_ion_losses_cd_.setParameters(p);

342

343

// y-ion loss model

344

it_begin = file.search(it_end, "YION_LOSS_MODEL_CR_BEGIN");

345

it_end = file.search(it_begin, "YION_LOSS_MODEL_CR_END");

346

p.clear();

347

HiddenMarkovModel y_ion_loss_hmm_cr;

348

parseHMMModel_(++it_begin, it_end, y_ion_loss_hmm_cr, p);

349

y_ion_losses_cr_.setHMM(y_ion_loss_hmm_cr);

350

y_ion_losses_cr_.setParameters(p);

351

352

it_begin = file.search(it_end, "YION_LOSS_MODEL_CD_BEGIN");

353

it_end = file.search(it_begin, "YION_LOSS_MODEL_CD_END");

354

p.clear();

355

HiddenMarkovModel y_ion_loss_hmm_cd;

356

parseHMMModel_(++it_begin, it_end, y_ion_loss_hmm_cd, p);

357

y_ion_losses_cd_.setHMM(y_ion_loss_hmm_cd);

358

y_ion_losses_cd_.setParameters(p);

359

360

valid_ = true;

361

return;

362

}

363

364

void PILISModel::writeGraphMLFile(const String& filename)

365

{

366

hmm_.writeGraphMLFile(filename);

367

return;

368

}

369

370

void PILISModel::writeToFile(const String& filename)

371

{

372

#ifdef SIM_DEBUG

373

cerr << "writing to file '" << filename << "'" << endl;

374

#endif

375

376

ofstream out(filename.c_str());

377

out << "BASE_MODEL_BEGIN" << endl;

378

writeParameters_(out, param_);

379

hmm_.write(out);

380

out << "BASE_MODEL_END" << endl;

381

382

out << "PRECURSOR_MODEL_CR_BEGIN" << endl;

383

writeParameters_(out, precursor_model_cr_.getParameters());

384

precursor_model_cr_.getHMM().write(out);

385

out << "PRECURSOR_MODEL_CR_END" << endl;

386

387

out << "PRECURSOR_MODEL_CD_BEGIN" << endl;

388

writeParameters_(out, precursor_model_cd_.getParameters());

389

precursor_model_cd_.getHMM().write(out);

390

out << "PRECURSOR_MODEL_CD_END" << endl;

391

392

out << "BION_LOSS_MODEL_CR_BEGIN" << endl;

393

writeParameters_(out, b_ion_losses_cr_.getParameters());

394

b_ion_losses_cr_.getHMM().write(out);

395

out << "BION_LOSS_MODEL_CR_END" << endl;

396

397

out << "BION_LOSS_MODEL_CD_BEGIN" << endl;

398

writeParameters_(out, b_ion_losses_cd_.getParameters());

399

b_ion_losses_cd_.getHMM().write(out);

400

out << "BION_LOSS_MODEL_CD_END" << endl;

401

402

out << "B2ION_LOSS_MODEL_CR_BEGIN" << endl;

403

writeParameters_(out, b2_ion_losses_cr_.getParameters());

404

b2_ion_losses_cr_.getHMM().write(out);

405

out << "B2ION_LOSS_MODEL_CR_END" << endl;

406

407

out << "B2ION_LOSS_MODEL_CD_BEGIN" << endl;

408

writeParameters_(out, b2_ion_losses_cd_.getParameters());

409

b2_ion_losses_cd_.getHMM().write(out);

410

out << "B2ION_LOSS_MODEL_CD_END" << endl;

411

412

out << "AION_LOSS_MODEL_CR_BEGIN" << endl;

413

writeParameters_(out, a_ion_losses_cr_.getParameters());

414

a_ion_losses_cr_.getHMM().write(out);

415

out << "AION_LOSS_MODEL_CR_END" << endl;

416

417

out << "AION_LOSS_MODEL_CD_BEGIN" << endl;

418

writeParameters_(out, a_ion_losses_cd_.getParameters());

419

a_ion_losses_cd_.getHMM().write(out);

420

out << "AION_LOSS_MODEL_CD_END" << endl;

421

422

out << "YION_LOSS_MODEL_CR_BEGIN" << endl;

423

writeParameters_(out, y_ion_losses_cr_.getParameters());

424

y_ion_losses_cr_.getHMM().write(out);

425

out << "YION_LOSS_MODEL_CR_END" << endl;

426

427

out << "YION_LOSS_MODEL_CD_BEGIN" << endl;

428

writeParameters_(out, y_ion_losses_cd_.getParameters());

429

y_ion_losses_cd_.getHMM().write(out);

430

out << "YION_LOSS_MODEL_CD_END" << endl;

431

432

return;

433

}

434

435

void PILISModel::train(const RichPeakSpectrum& in_spec, const AASequence& peptide, UInt charge)

436

{

437

if (!valid_)

438

{

439

cerr << "PILISModel: cannot train, initialize model from file first, e.g. data/PILIS/PILIS_model_default.dat" << endl;

440

return;

441

}

442

443

if (peptide.size() >= (Size)param_.getValue("visible_model_depth"))

444

{

445

cerr << "PILISModel: cannot train peptide " << peptide << " of length " << peptide.size() << " (max for this model is " << param_.getValue("visible_model_depth") << ", as defined by parameter \"visible_model_depth\")" << endl;

446

return;

447

}

448

449

RichPeakSpectrum train_spec = in_spec;

450

train_spec.sortByPosition();

451

452

#ifdef TRAINING_DEBUG

453

cout << "peptide: " << peptide << "(z=" << charge << ")" << endl;

454

#endif

455

456

// get proton distribution

457

vector<DoubleReal> bb_charge_full, sc_charge_full;

458

prot_dist_.getProtonDistribution(bb_charge_full, sc_charge_full, peptide, charge, Residue::YIon);

459

prot_dist_.setPeptideProtonDistribution(bb_charge_full, sc_charge_full);

460

461

// get start probabilities

462

vector<DoubleReal> bb_init, sc_init, cr_init;

463

DoubleReal precursor_init(0);

464

bool is_charge_remote = getInitialTransitionProbabilities_(bb_init, cr_init, sc_init, precursor_init, bb_charge_full, sc_charge_full, peptide);

465

466

// clear the main Hidden Markov Model

467

hmm_.clearInitialTransitionProbabilities();

468

hmm_.clearTrainingEmissionProbabilities();

469

470

//DoubleReal harge_sum(0);

471

vector<AASequence> prefixes, suffixes;

472

473

DoubleReal pep_weight(peptide.getMonoWeight());

474

DoubleReal peptide_mz((pep_weight + charge)/(DoubleReal)charge);

475

476

// for each site:

477

// 1. set proton distribution,

478

// 2. initial training intensities,

479

// 3. train the model

480

for (Size i = 0; i != peptide.size() - 1; ++i)

481

{

482

String pos_name, prefix_size(i + 1), suffix_size(peptide.size() - 1 - i);

483

484

if (i < floor((peptide.size() - 1.0)/2.0))

485

{

486

pos_name = prefix_size;

487

}

488

else

489

{

490

pos_name = "k-"+suffix_size;

491

}

492

493

AASequence prefix(peptide.getPrefix(i + 1)), suffix(peptide.getSuffix(peptide.size() - i - 1));

494

AASequence aa1_seq, aa2_seq;

495

aa1_seq += &peptide[i];

496

aa2_seq += &peptide[i + 1];

497

String aa1(aa1_seq.toString()), aa2(aa2_seq.toString());

498

499

// calc PAs and get b/y ratios for bxyz pathway

500

vector<DoubleReal> b_cr_ints(charge, 0), y_cr_ints(charge, 0), b_sc_ints(charge, 0), y_sc_ints(charge, 0), b_ints(charge, 0), y_ints(charge, 0), a_ints(charge, 0), ay_ints(charge, 0);

501

prot_dist_.getChargeStateIntensities(peptide, prefix, suffix, charge, Residue::BIon, b_ints, y_ints, ProtonDistributionModel::ChargeDirected);

502

503

prot_dist_.getChargeStateIntensities(peptide, prefix, suffix, charge, Residue::AIon, a_ints, ay_ints, ProtonDistributionModel::ChargeDirected);

504

505

if (!is_charge_remote)

506

{

507

hmm_.setInitialTransitionProbability("BB"+pos_name, bb_init[i]);

508

hmm_.setInitialTransitionProbability(aa1+aa2+"_bxyz"+pos_name, bb_init[i]);

509

hmm_.setInitialTransitionProbability(aa1+aa2+"_axyz"+pos_name, bb_init[i]);

510

}

511

512

prefixes.push_back(prefix);

513

suffixes.push_back(suffix);

514

515

if ((aa1 == "D" || aa1 == "E" || pos_name == "k-1" || pos_name == "k-2") && is_charge_remote)

516

{

517

hmm_.setInitialTransitionProbability("CR"+pos_name, cr_init[i]);

518

prot_dist_.getChargeStateIntensities(peptide, prefix, suffix, charge, Residue::BIon, b_cr_ints, y_cr_ints, ProtonDistributionModel::ChargeRemote);

519

//cerr << "ChargeStats: CR=" << CR_charges[i] << ", " << peptide.getPrefix(i+1) << "-" << peptide.getSuffix(peptide.size() - 1 - i) << ", " << charge << endl;

520

}

521

522

if ((aa1 == "K" || aa1 == "H" || aa1 == "R") && is_charge_remote)

523

{

524

prot_dist_.getChargeStateIntensities(peptide, prefix, suffix, charge, Residue::BIon, b_sc_ints, y_sc_ints, ProtonDistributionModel::SideChain);

525

hmm_.setInitialTransitionProbability("SC"+pos_name, /*sc_charge_full[i]*/ sc_init[i]);

526

//cerr << "ChargeStats: SC=" << SC_charges[i] << ", " << peptide.getPrefix(i+1) << "-" << peptide.getSuffix(peptide.size() - 1 - i) << ", " << charge << endl;

527

}

528

529

530

UInt max_fragment_charge_training = param_.getValue("max_fragment_charge_training");

531

Map<int, DoubleReal> sum_a, sum_b, sum_y;

532

DoubleReal sum_a_ints(0.0), sum_b_ints(0.0), sum_y_ints(0.0);

533

for (UInt z = 1; z <= max_fragment_charge_training && z <= charge; ++z)

534

{

535

sum_a[z] = 0;

536

sum_b[z] = 0;

537

sum_y[z] = 0;

538

DoubleReal charge_remote_threshold((DoubleReal)param_.getValue("charge_remote_threshold"));

539

DoubleReal avail_bb_sum_prefix = getAvailableBackboneCharge_(prefix, Residue::BIon, z);

540

DoubleReal avail_bb_sum_suffix = getAvailableBackboneCharge_(suffix, Residue::YIon, z);

541

542

if (prefix.size() != 2)

543

{

544

DoubleReal pre_weight = prefix.getMonoWeight(Residue::BIon);

545

546

if (avail_bb_sum_prefix <= charge_remote_threshold)

547

{

548

//cerr << "Train b-ion losses, CR (avail=" << avail_bb_sum_prefix << ", z=" << z << ", prefix=" << prefix << endl;

549

sum_b[z] += b_ion_losses_cr_.train(in_spec, prefix, pre_weight, z, pep_weight);

550

}

551

else

552

{

553

//cerr << "Train b-ion losses, CD (avail=" << avail_bb_sum_prefix << ", z=" << z << ", prefix=" << prefix << endl;

554

sum_b[z] += b_ion_losses_cd_.train(in_spec, prefix, pre_weight, z, pep_weight);

555

}

556

}

557

else

558

{

559

DoubleReal pre_weight = prefix.getMonoWeight(Residue::BIon);

560

561

if (avail_bb_sum_prefix <= charge_remote_threshold)

562

{

563

sum_b[z] += b2_ion_losses_cr_.train(in_spec, prefix, pre_weight, z, pep_weight);

564

}

565

else

566

{

567

sum_b[z] += b2_ion_losses_cd_.train(in_spec, prefix, pre_weight, z, pep_weight);

568

}

569

}

570

571

DoubleReal a_pre_weight = prefix.getMonoWeight(Residue::AIon);

572

if (avail_bb_sum_prefix <= charge_remote_threshold)

573

{

574

sum_a[z] += a_ion_losses_cr_.train(in_spec, prefix, a_pre_weight, z, pep_weight);

575

}

576

else

577

{

578

sum_a[z] += a_ion_losses_cd_.train(in_spec, prefix, a_pre_weight, z, pep_weight);

579

}

580

581

DoubleReal suf_weight = suffix.getMonoWeight(Residue::YIon);

582

if (avail_bb_sum_suffix <= charge_remote_threshold)

583

{

584

//cerr << "Train y-ion losses, CR (avail=" << avail_bb_sum_suffix << ", z=" << z << ", suffix=" << suffix << endl;

585

sum_y[z] += y_ion_losses_cr_.train(in_spec, suffix, suf_weight, z, pep_weight);

586

}

587

else

588

{

589

//cerr << "Train y-ion losses, CD (avail=" << avail_bb_sum_suffix << ", z=" << z << ", suffix=" << suffix << endl;

590

sum_y[z] += y_ion_losses_cd_.train(in_spec, suffix, suf_weight, z, pep_weight);

591

}

592

593

sum_a_ints += sum_a[z];

594

sum_b_ints += sum_b[z];

595

sum_y_ints += sum_y[z];

596

}

597

598

// end state is always needed

599

hmm_.setTrainingEmissionProbability("end"+pos_name, 0.5/(DoubleReal(peptide.size() - 1)));

600

if (!is_charge_remote)

601

{

602

//hmm_.setTrainingEmissionProbability("AA"+pos_name, sum_a_ints + sum_b_ints + sum_y_ints);

603

hmm_.enableTransition("AA" + pos_name, aa1 + aa2 + "_bxyz" + pos_name);

604

hmm_.enableTransition("AA" + pos_name, aa1 + aa2 + "_axyz" + pos_name);

605

606

// now enable the states

607

hmm_.enableTransition("BB"+pos_name, "AA"+pos_name);

608

hmm_.enableTransition("BB"+pos_name, "end"+pos_name);

609

610

// bxyz path

611

hmm_.enableTransition(aa1+aa2+"_bxyz" + pos_name, "bxyz"+pos_name);

612

hmm_.enableTransition(aa1+aa2+"_bxyz" + pos_name, "end"+pos_name);

613

614

DoubleReal pre_emission_prob(0), suf_emission_prob(0);

615

for (UInt z = 1; z <= max_fragment_charge_training && z <= charge; ++z)

616

{

617

pre_emission_prob += b_ints[z - 1] * sum_b[z];

618

suf_emission_prob += y_ints[z - 1] * sum_y[z];

619

}

620

hmm_.setTransitionProbability("bxyz" + pos_name, "bxyz_" + pos_name + "-prefix", 1.0);

621

hmm_.setTransitionProbability("bxyz" + pos_name, "bxyz_" + pos_name + "-suffix", 1.0);

622

hmm_.enableTransition("bxyz_" + pos_name + "-prefix", "bxyz_" + pos_name + "-prefix-ions");

623

hmm_.enableTransition("bxyz_" + pos_name + "-prefix", "end" + pos_name);

624

hmm_.enableTransition("bxyz_" + pos_name + "-suffix", "bxyz_" + pos_name + "-suffix-ions");

625

hmm_.enableTransition("bxyz_" + pos_name + "-suffix", "end" + pos_name);

626

627

hmm_.setTrainingEmissionProbability("bxyz_" + pos_name + "-prefix-ions", pre_emission_prob);

628

hmm_.setTrainingEmissionProbability("bxyz_" + pos_name + "-suffix-ions", suf_emission_prob);

629

630

// axyz path

631

hmm_.enableTransition(aa1+aa2+"_axyz"+pos_name, "axyz"+pos_name);

632

hmm_.enableTransition(aa1+aa2+"_axyz"+pos_name, "end"+pos_name);

633

634

pre_emission_prob = 0;

635

suf_emission_prob = 0;

636

for (UInt z = 1; z <= max_fragment_charge_training && z <= charge; ++z)

637

{

638

pre_emission_prob += a_ints[z - 1] * sum_a[z];

639

//suf_emission_prob += ay_ints[z] * sum_y[z]; // TODO test if this can be done, or ruins the a-ion intensities

640

}

641

hmm_.setTransitionProbability("axyz" + pos_name, "axyz_" + pos_name + "-prefix", 1.0);

642

hmm_.enableTransition("axyz_" + pos_name + "-prefix", "axyz_" + pos_name + "-prefix-ions");

643

hmm_.enableTransition("axyz_" + pos_name + "-prefix", "end" + pos_name);

644

hmm_.setTrainingEmissionProbability("axyz_" + pos_name + "-prefix-ions", pre_emission_prob);

645

646

}

647

648

const String cr_aa("DE");

649

for (String::ConstIterator it = cr_aa.begin(); it != cr_aa.end(); ++it)

650

{

651

if (is_charge_remote && aa1 == String(*it))

652

{

653

hmm_.enableTransition("CR"+pos_name, "A"+pos_name);

654

hmm_.enableTransition("CR"+pos_name, "end"+pos_name);

655

656

DoubleReal pre_emission_prob(0), suf_emission_prob(0);

657

for (UInt z = 1; z <= max_fragment_charge_training && z <= charge; ++z)

658

{

659

pre_emission_prob += b_cr_ints[z - 1] * sum_b[z];

660

suf_emission_prob += y_cr_ints[z - 1] * sum_y[z];

661

}

662

663

//hmm_.setTrainingEmissionProbability("A"+pos_name, pre_emission_prob + suf_emission_prob);

664

hmm_.enableTransition("A" + pos_name, aa2 + "_" + String(*it) + pos_name);

665

666

hmm_.setTransitionProbability(String(*it) + pos_name, String(*it) + "_" + pos_name + "-prefix", 1.0);

667

hmm_.enableTransition(String(*it) + "_" + pos_name + "-prefix", String(*it) + "_" + pos_name + "-prefix-ions");

668

hmm_.enableTransition(String(*it) + "_" + pos_name + "-prefix", "end" + pos_name);

669

hmm_.setTrainingEmissionProbability(String(*it) + "_" + pos_name + "-prefix-ions", pre_emission_prob);

670

671

hmm_.setTransitionProbability(String(*it) + pos_name, String(*it) + "_" + pos_name + "-suffix", 1.0);

672

hmm_.enableTransition(String(*it) + "_" + pos_name + "-suffix", String(*it) + "_" + pos_name + "-suffix-ions");

673

hmm_.enableTransition(String(*it) + "_" + pos_name + "-suffix", "end" + pos_name);

674

hmm_.setTrainingEmissionProbability(String(*it) + "_" + pos_name + "-suffix-ions", suf_emission_prob);

675

676

hmm_.setInitialTransitionProbability(aa2+"_" + String(*it) + pos_name, cr_init[i]);

677

hmm_.enableTransition(aa2+"_" + String(*it) + pos_name, String(*it) + pos_name);

678

hmm_.enableTransition(aa2+"_" + String(*it) + pos_name, "end" + pos_name);

679

}

680

}

681

682

if (is_charge_remote && !peptide.has("D") && !peptide.has("E"))

683

{

684

685

if (pos_name == "k-1")

686

{

687

hmm_.enableTransition("CRk-1", "Ak-1");

688

hmm_.enableTransition("CRk-1", "endk-1");

689

690

DoubleReal pre_emission_prob(0), suf_emission_prob(0);

691

for (UInt z = 1; z <= max_fragment_charge_training && z <= charge; ++z)

692

{

693

pre_emission_prob += b_cr_ints[z - 1] * sum_b[z];

694

suf_emission_prob += y_cr_ints[z - 1] * sum_y[z];

695

}

696

697

//hmm_.setTrainingEmissionProbability("Ak-1", pre_emission_prob + suf_emission_prob);

698

hmm_.enableTransition("Ak-1", aa1 + "_bk-1");

699

700

701

hmm_.setTransitionProbability("bk-1", "bk-1_-prefix", 1.0);

702

hmm_.enableTransition("bk-1_-prefix", "bk-1_-prefix-ions");

703

hmm_.enableTransition("bk-1_-prefix", "endk-1");

704

hmm_.setTrainingEmissionProbability("bk-1_-prefix-ions", pre_emission_prob);

705

706

hmm_.setTransitionProbability("bk-1", "bk-1_-suffix", 1.0);

707

hmm_.enableTransition("bk-1_-suffix", "bk-1_-suffix-ions");

708

hmm_.enableTransition("bk-1_-suffix", "endk-1");

709

hmm_.setTrainingEmissionProbability("bk-1_-suffix-ions", suf_emission_prob);

710

711

hmm_.setTransitionProbability("bk-1", "bk-1_-ions", 1.0);

712

713

hmm_.setInitialTransitionProbability(aa1+"_bk-1", cr_init[i]);

714

hmm_.enableTransition(aa1+"_bk-1", "bk-1");

715

hmm_.enableTransition(aa1+"_bk-1", "endk-1");

716

}

717

718

if (pos_name == "k-2")

719

{

720

hmm_.enableTransition("CRk-2", "Ak-2");

721

hmm_.enableTransition("CRk-2", "endk-2");

722

723

DoubleReal pre_emission_prob(0), suf_emission_prob(0);

724

for (UInt z = 1; z <= max_fragment_charge_training && z <= charge; ++z)

725

{

726

pre_emission_prob += b_cr_ints[z - 1] * sum_b[z];

727

suf_emission_prob += y_cr_ints[z - 1] * sum_y[z];

728

}

729

730

//hmm_.setTrainingEmissionProbability("Ak-2", pre_emission_prob + suf_emission_prob);

731

hmm_.enableTransition("Ak-2", aa1 + "_bk-2");

732

733

734

hmm_.setTransitionProbability("bk-2", "bk-2_-prefix", 1.0);

735

hmm_.enableTransition("bk-2_-prefix", "bk-2_-prefix-ions");

736

hmm_.enableTransition("bk-2_-prefix", "endk-2");

737

hmm_.setTrainingEmissionProbability("bk-2_-prefix-ions", pre_emission_prob);

738

739

hmm_.setTransitionProbability("bk-2", "bk-2_-suffix", 1.0);

740

hmm_.enableTransition("bk-2_-suffix", "bk-2_-suffix-ions");

741

hmm_.enableTransition("bk-2_-suffix", "endk-2");

742

hmm_.setTrainingEmissionProbability("bk-2_-suffix-ions", suf_emission_prob);

743

744

hmm_.setTransitionProbability("bk-2", "bk-2_-ions", 1.0);

745

746

hmm_.setInitialTransitionProbability(aa1+"_bk-2", cr_init[i]);

747

hmm_.enableTransition(aa1+"_bk-2", "bk-2");

748

hmm_.enableTransition(aa1+"_bk-2", "endk-2");

749

}

750

}

751

752

const String sc_aa("KRH");

753

for (String::ConstIterator it = sc_aa.begin(); it != sc_aa.end(); ++it)

754

{

755

if (is_charge_remote && aa1 == String(*it))

756

{

757

758

hmm_.enableTransition("SC"+pos_name, "ASC"+pos_name);

759

hmm_.enableTransition("SC"+pos_name, "end"+pos_name);

760

761

DoubleReal pre_emission_prob(0), suf_emission_prob(0);

762

for (UInt z = 1; z <= max_fragment_charge_training && z <= charge; ++z)

763

{

764

pre_emission_prob += b_sc_ints[z - 1] * sum_b[z];

765

suf_emission_prob += y_sc_ints[z - 1] * sum_y[z];

766

}

767

768

//hmm_.setTrainingEmissionProbability("ASC"+pos_name, pre_emission_prob + suf_emission_prob);

769

hmm_.enableTransition("ASC" + pos_name, aa2 + "_" + String(*it) + pos_name);

770

771

hmm_.setTransitionProbability(String(*it) + pos_name, String(*it) + "_" + pos_name + "-prefix", 1.0);

772

hmm_.setTransitionProbability(String(*it) + pos_name, String(*it) + "_" + pos_name + "-suffix", 1.0);

773

774

hmm_.enableTransition(String(*it) + "_" + pos_name + "-prefix", String(*it) + "_" + pos_name + "-prefix-ions");

775

hmm_.enableTransition(String(*it) + "_" + pos_name + "-prefix", "end" + pos_name);

776

hmm_.enableTransition(String(*it) + "_" + pos_name + "-suffix", String(*it) + "_" + pos_name + "-suffix-ions");

777

hmm_.enableTransition(String(*it) + "_" + pos_name + "-suffix", "end" + pos_name);

778

779

hmm_.setTrainingEmissionProbability(String(*it) + "_" + pos_name + "-prefix-ions", pre_emission_prob);

780

hmm_.setTrainingEmissionProbability(String(*it) + "_" + pos_name + "-suffix-ions", suf_emission_prob);

781

782

783

hmm_.setInitialTransitionProbability(aa2 + "_" + String(*it) + pos_name, sc_init[i]);

784

hmm_.enableTransition(aa2+"_" + String(*it) + pos_name, String(*it) + pos_name);

785

hmm_.enableTransition(aa2+"_" + String(*it) + pos_name, "end"+pos_name);

786

787

}

788

}

789

}

790

791

if (is_charge_remote)

792

{

793

precursor_model_cr_.train(in_spec, peptide, peptide_mz, charge, pep_weight);

794

}

795

else

796

{

797

precursor_model_cd_.train(in_spec, peptide, peptide_mz, charge, pep_weight);

798

}

799

800

// now train the model with the data set

801

hmm_.train();

802

hmm_.disableTransitions();

803

804

return;

805

}

806

807

void PILISModel::evaluate()

808

{

809

hmm_.evaluate();

810

precursor_model_cr_.evaluate();

811

precursor_model_cd_.evaluate();

812

a_ion_losses_cr_.evaluate();

813

a_ion_losses_cd_.evaluate();

814

b_ion_losses_cr_.evaluate();

815

b_ion_losses_cd_.evaluate();

816

b2_ion_losses_cr_.evaluate();

817

b2_ion_losses_cd_.evaluate();

818

y_ion_losses_cr_.evaluate();

819

y_ion_losses_cd_.evaluate();

820

821

hmm_.setVariableModifications((StringList)param_.getValue("variable_modifications"));

822

hmm_.estimateUntrainedTransitions();

823

}

824

825

void PILISModel::getSpectrum(RichPeakSpectrum& spec, const AASequence& peptide, UInt charge)

826

{

827

if (!valid_)

828

{

829

cerr << "PILISModel: cannot simulate, initialize model from file first, e.g. data/PILIS/PILIS_model_default.dat" << endl;

830

return;

831

}

832

833

if (peptide.size() > (Size)param_.getValue("visible_model_depth"))

834

{

835

cerr << "PILISModel: cannot generate spectra of peptide '" << peptide << "' of length " << peptide.size() << " (max of this model is " << param_.getValue("visible_model_depth") << ", as defined by \"visible_model_depth\")" << endl;

836

return;

837

}

838

839

// calc proton distribution

840

vector<DoubleReal> sc_charge_full, bb_charge_full;

841

//cerr << "getProtonDistribution: " << peptide << " " << charge << endl;

842

prot_dist_.getProtonDistribution(bb_charge_full, sc_charge_full, peptide, charge, Residue::YIon);

843

prot_dist_.setPeptideProtonDistribution(bb_charge_full, sc_charge_full);

844

845

#ifdef INIT_CHARGE_DEBUG

846

for (Size i = 0; i != bb_charge_full.size() - 1; ++i)

847

{

848

cerr << "i: bb=" << bb_charge_full[i] << ", sc=" << sc_charge_full[i] << endl;

849

}

850

#endif

851

852

hmm_.clearInitialTransitionProbabilities();

853

hmm_.clearTrainingEmissionProbabilities();

854

855

// set charges

856

vector<DoubleReal> bb_init, sc_init, cr_init;

857

DoubleReal precursor_init(0);

858

bool is_charge_remote = getInitialTransitionProbabilities_(bb_init, cr_init, sc_init, precursor_init, bb_charge_full, sc_charge_full, peptide);

859

#ifdef INIT_CHARGE_DEBUG

860

cerr << "is_charge_remote=" << is_charge_remote << endl;

861

for (Size i = 0; i != bb_init.size(); ++i)

862

{

863

cerr << "bb=" << bb_init[i] << ", cr=" << cr_init[i] << ", sc=" << sc_init[i] << endl;

864

}

865

cerr << "precursor_init=" << precursor_init << endl;

866

#endif

867

868

vector<AASequence> suffixes, prefixes;

869

vector<String> pos_names;

870

871

vector<vector<DoubleReal> > all_b_cr_ints, all_y_cr_ints, all_b_sc_ints, all_y_sc_ints, all_b_ints, all_y_ints, all_a_ints, all_ay_ints;

872

// get the paths

873

for (Size i = 0; i != peptide.size() - 1; ++i)

874

{

875

String pos_name, y_name1, b_name1, a_name1, y_name, b_name;

876

String y_name2, b_name2, a_name2;

877

String i_plus1(i+1), pep_size_i(peptide.size() - 1 - i);

878

879

if (i < floor((peptide.size() - 1.0)/2.0))

880

{

881

pos_name = i_plus1;

882

}

883

else

884

{

885

pos_name = "k-"+pep_size_i;

886

}

887

888

AASequence prefix(peptide.getPrefix(i + 1)), suffix(peptide.getSuffix(peptide.size() - 1 - i));

889

suffixes.push_back(suffix);

890

prefixes.push_back(prefix);

891

pos_names.push_back(pos_name);

892

893

AASequence aa1_seq, aa2_seq;

894

aa1_seq += &peptide[i];

895

aa2_seq += &peptide[i + 1];

896

String aa1(aa1_seq.toString()), aa2(aa2_seq.toString());

897

898

hmm_.setInitialTransitionProbability("BB"+pos_name, bb_init[i]);

899

900

vector<DoubleReal> b_cr_ints(charge, 0), y_cr_ints(charge, 0), b_sc_ints(charge, 0), y_sc_ints(charge, 0), b_ints(charge, 0), y_ints(charge, 0), a_ints(charge, 0), ay_ints(charge, 0);

901

if ((aa1 == "D" || aa1 == "E" || pos_name == "k-1" || pos_name == "k-2") && is_charge_remote)

902

{

903

hmm_.setInitialTransitionProbability("CR"+pos_name, cr_init[i]);

904

prot_dist_.getChargeStateIntensities(peptide, prefix, suffix, charge, Residue::BIon, b_cr_ints, y_cr_ints, ProtonDistributionModel::ChargeRemote);

905

}

906

907

if ((aa1 == "K" || aa1 == "H" || aa1 == "R") && is_charge_remote)

908

{

909

hmm_.setInitialTransitionProbability("SC"+pos_name, sc_init[i]);

910

prot_dist_.getChargeStateIntensities(peptide, prefix, suffix, charge, Residue::BIon, b_sc_ints, y_sc_ints, ProtonDistributionModel::SideChain);

911

}

912

913

prot_dist_.getChargeStateIntensities(peptide, prefix, suffix, charge, Residue::BIon, b_ints, y_ints, ProtonDistributionModel::ChargeDirected);

914

prot_dist_.getChargeStateIntensities(peptide, prefix, suffix, charge, Residue::AIon, a_ints, ay_ints, ProtonDistributionModel::ChargeDirected);

915

916

//if (charge == 3)

917

//{

918

#ifdef INIT_CHARGE_DEBUG

919

cerr << prefix << " - " << suffix << " ";

920

for (Size ions_charge = 0; ions_charge != b_ints.size(); ++ions_charge)

921

{

922

cerr << b_ints[ions_charge] << " - " << y_ints[ions_charge] << " | ";

923

}

924

cerr << endl;

925

#endif

926

//}

927

928

all_a_ints.push_back(a_ints);

929

all_ay_ints.push_back(ay_ints);

930

all_b_ints.push_back(b_ints);

931

all_b_cr_ints.push_back(b_cr_ints);

932

all_y_cr_ints.push_back(y_cr_ints);

933

all_b_sc_ints.push_back(b_sc_ints);

934

all_y_sc_ints.push_back(y_sc_ints);

935

all_y_ints.push_back(y_ints);

936

937

// now enable the states

938

hmm_.enableTransition("BB"+pos_name, "AA"+pos_name);

939

hmm_.enableTransition("BB"+pos_name, "end"+pos_name);

940

941

hmm_.enableTransition("AA"+pos_name, aa1+aa2+"_bxyz"+pos_name);

942

hmm_.enableTransition(aa1+aa2+"_bxyz"+pos_name, "bxyz"+pos_name);

943

hmm_.enableTransition(aa1+aa2+"_bxyz"+pos_name, "end"+pos_name);

944

945

hmm_.setTransitionProbability("bxyz" + pos_name, "bxyz_" + pos_name + "-prefix", 1.0);

946

hmm_.setTransitionProbability("bxyz" + pos_name, "bxyz_" + pos_name + "-suffix", 1.0);

947

948

hmm_.enableTransition("bxyz_" + pos_name + "-prefix", "bxyz_" + pos_name + "-prefix-ions");

949

hmm_.enableTransition("bxyz_" + pos_name + "-prefix", "end" + pos_name);

950

hmm_.enableTransition("bxyz_" + pos_name + "-suffix", "bxyz_" + pos_name + "-suffix-ions");

951

hmm_.enableTransition("bxyz_" + pos_name + "-suffix", "end" + pos_name);

952

953

954

// axyz

955

hmm_.enableTransition("AA"+pos_name, aa1+aa2+"_axyz"+pos_name);

956

hmm_.enableTransition(aa1+aa2+"_axyz"+pos_name, "axyz"+pos_name);

957

hmm_.enableTransition(aa1+aa2+"_axyz"+pos_name, "end"+pos_name);

958

hmm_.setTransitionProbability("axyz" + pos_name, "axyz_" + pos_name + "-prefix", 1.0);

959

hmm_.enableTransition("axyz_" + pos_name + "-prefix", "axyz_" + pos_name + "-prefix-ions");

960

hmm_.enableTransition("axyz_" + pos_name + "-prefix", "end" + pos_name);

961

962

const String cr_aa("DE");

963

for (String::ConstIterator it = cr_aa.begin(); it != cr_aa.end(); ++it)

964

{

965

if (is_charge_remote && aa1 == String(*it))

966

{

967

hmm_.setTransitionProbability(String(*it) + pos_name, String(*it) + "_" + pos_name + "-prefix", 1.0);

968

hmm_.setTransitionProbability(String(*it) + pos_name, String(*it) + "_" + pos_name + "-suffix", 1.0);

969

970

hmm_.enableTransition(String(*it) + "_" + pos_name + "-prefix", String(*it) + "_" + pos_name + "-prefix-ions");

971

hmm_.enableTransition(String(*it) + "_" + pos_name + "-prefix", "end" + pos_name);

972

hmm_.enableTransition(String(*it) + "_" + pos_name + "-suffix", String(*it) + "_" + pos_name + "-suffix-ions");

973

hmm_.enableTransition(String(*it) + "_" + pos_name + "-suffix", "end" + pos_name);

974

975

hmm_.enableTransition("CR"+pos_name, "A"+pos_name);

976

hmm_.enableTransition("CR"+pos_name, "end"+pos_name);

977

hmm_.enableTransition("A"+pos_name, aa2+"_" + String(*it) + pos_name);

978

hmm_.enableTransition(aa2+"_" + String(*it) + pos_name, String(*it) + pos_name);

979

hmm_.enableTransition(aa2+"_" + String(*it) + pos_name, "end"+pos_name);

980

}

981

}

982

983

if (!peptide.has("D") && !peptide.has("E") && is_charge_remote)

984

{

985

if (pos_name == "k-1")

986

{

987

hmm_.setTransitionProbability("bk-1", "bk-1_-prefix", 1.0);

988

hmm_.setTransitionProbability("bk-1", "bk-1_-suffix", 1.0);

989

990

hmm_.enableTransition("bk-1_-prefix", "bk-1_-prefix-ions");

991

hmm_.enableTransition("bk-1_-prefix", "endk-1" );

992

hmm_.enableTransition("bk-1_-suffix", "bk-1_-suffix-ions");

993

hmm_.enableTransition("bk-1_-suffix", "endk-1");

994

995

hmm_.enableTransition("CRk-1", "Ak-1");

996

hmm_.enableTransition("CRk-1", "endk-1");

997

hmm_.enableTransition("Ak-1", aa1+"_bk-1");

998

999

hmm_.enableTransition(aa1 + "_bk-1", "bk-1");

1000

hmm_.enableTransition(aa1 + "_bk-1", "endk-1");

1001

}

1002

1003

if (pos_name == "k-2")

1004

{

1005

hmm_.setTransitionProbability("bk-2", "bk-2_-prefix", 1.0);

1006

hmm_.setTransitionProbability("bk-2", "bk-2_-suffix", 1.0);

1007

1008

hmm_.enableTransition("bk-2_-prefix", "bk-2_-prefix-ions");

1009

hmm_.enableTransition("bk-2_-prefix", "endk-2" );

1010

hmm_.enableTransition("bk-2_-suffix", "bk-2_-suffix-ions");

1011

hmm_.enableTransition("bk-2_-suffix", "endk-2");

1012

1013

hmm_.enableTransition("CRk-2", "Ak-2");

1014

hmm_.enableTransition("CRk-2", "endk-2");

1015

hmm_.enableTransition("Ak-2", aa1+"_bk-2");

1016

1017

hmm_.enableTransition(aa1 + "_bk-2", "bk-2");

1018

hmm_.enableTransition(aa1 + "_bk-2", "endk-2");

1019

}

1020

}

1021

1022

const String sc_aa("HKR");

1023

for (String::ConstIterator it = sc_aa.begin(); it != sc_aa.end(); ++it)

1024

{

1025

if (is_charge_remote && aa1 == String(*it))

1026

{

1027

hmm_.setTransitionProbability(String(*it) + pos_name, String(*it) + "_" + pos_name + "-prefix", 1.0);

1028

hmm_.setTransitionProbability(String(*it) + pos_name, String(*it) + "_" + pos_name + "-suffix", 1.0);

1029

1030

hmm_.enableTransition(String(*it) + "_" + pos_name + "-prefix", String(*it) + "_" + pos_name + "-prefix-ions");

1031

hmm_.enableTransition(String(*it) + "_" + pos_name + "-prefix", "end" + pos_name);

1032

hmm_.enableTransition(String(*it) + "_" + pos_name + "-suffix", String(*it) + "_" + pos_name + "-suffix-ions");

1033

hmm_.enableTransition(String(*it) + "_" + pos_name + "-suffix", "end" + pos_name);

1034

1035

1036

1037

hmm_.enableTransition("SC"+pos_name, "ASC"+pos_name);

1038

hmm_.enableTransition("SC"+pos_name, "end"+pos_name);

1039

hmm_.enableTransition("ASC"+pos_name, aa2+"_" + String(*it) + pos_name);

1040

hmm_.enableTransition(aa2+"_" + String(*it) + pos_name, String(*it) + pos_name);

1041

hmm_.enableTransition(aa2+"_" + String(*it) + pos_name, "end"+pos_name);

1042

}

1043

}

1044

}

1045

1046

//cerr << "Collecting peaks..." << endl;

1047

1048

Map<HMMState*, DoubleReal> tmp;

1049

hmm_.calculateEmissionProbabilities(tmp);

1050

1051

// clear peaks from last spectrum

1052

peaks_.clear();

1053

1054

//stringstream peptide_ss;

1055

//peptide_ss << peptide;

1056

//hmm_.writeGraphMLFile(String("model_graph_train_"+peptide_ss.str()+"_"+String(charge)+".graphml").c_str());

1057

1058

UInt max_isotope = (UInt)param_.getValue("max_isotope");

1059

UInt max_fragment_charge = (UInt)param_.getValue("max_fragment_charge");

1060

DoubleReal charge_remote_threshold = (DoubleReal)param_.getValue("charge_remote_threshold");

1061

for (Size i = 0; i != prefixes.size(); ++i)

1062

{

1063

String aa;

1064

AASequence aa_seq;

1065

aa_seq += peptide[i].getOneLetterCode();

1066

aa = aa_seq.toString();

1067

1068

vector<DoubleReal> pre_path_intensities, suf_path_intensities;

1069

pre_path_intensities.push_back(tmp[hmm_.getState("bxyz_" + pos_names[i] + "-prefix-ions")]);

1070

suf_path_intensities.push_back(tmp[hmm_.getState("bxyz_" + pos_names[i] + "-suffix-ions")]);

1071

#ifdef INIT_CHARGE_DEBUG

1072

cerr << prefixes[i] << " - " << suffixes[i] << ": bxyz=" << pre_path_intensities.back() << "|" << suf_path_intensities.back() << " ";

1073

#endif

1074

1075

const String cr_and_sc("DEHKR");

1076

for (String::ConstIterator it = cr_and_sc.begin(); it != cr_and_sc.end(); ++it)

1077

{

1078

if (aa == String(*it))

1079

{

1080

pre_path_intensities.push_back(tmp[hmm_.getState(String(*it) + "_" + pos_names[i] + "-prefix-ions")]);

1081

suf_path_intensities.push_back(tmp[hmm_.getState(String(*it) + "_" + pos_names[i] + "-suffix-ions")]);

1082

#ifdef INIT_CHARGE_DEBUG

1083

cerr << ": " << aa << "=" << pre_path_intensities.back() << "|" << suf_path_intensities.back() << " ";

1084

#endif

1085

}

1086

}

1087

1088

if (i == peptide.size() - 2)

1089

{

1090

pre_path_intensities.push_back(tmp[hmm_.getState("bk-1_-prefix-ions")]);

1091

suf_path_intensities.push_back(tmp[hmm_.getState("bk-1_-suffix-ions")]);

1092

}

1093

if (i == peptide.size() - 3)

1094

{

1095

pre_path_intensities.push_back(tmp[hmm_.getState("bk-2_-prefix-ions")]);

1096

suf_path_intensities.push_back(tmp[hmm_.getState("bk-2_-suffix-ions")]);

1097

}

1098

1099

pre_path_intensities.push_back(tmp[hmm_.getState("axyz_" + pos_names[i] + "-prefix-ions")]);

1100

suf_path_intensities.push_back(tmp[hmm_.getState("axyz_" + pos_names[i] + "-suffix-ions")]);

1101

1102

vector<vector<vector<DoubleReal> > > prefix_intensities, suffix_intensities;

1103

prefix_intensities.push_back(all_b_ints);

1104

suffix_intensities.push_back(all_y_ints);

1105

const String sc_aa("HKR");

1106

if (sc_aa.hasSubstring(aa))

1107

{

1108

prefix_intensities.push_back(all_b_sc_ints);

1109

suffix_intensities.push_back(all_y_sc_ints);

1110

}

1111

const String cr_aa("DE");

1112

if (cr_aa.hasSubstring(aa))

1113

{

1114

prefix_intensities.push_back(all_b_cr_ints);

1115

suffix_intensities.push_back(all_y_cr_ints);

1116

}

1117

1118

if (i == peptide.size() - 2 || i == peptide.size() - 3)

1119

{

1120

prefix_intensities.push_back(all_b_cr_ints);

1121

suffix_intensities.push_back(all_y_cr_ints);

1122

}

1123

1124

prefix_intensities.push_back(all_a_ints);

1125

suffix_intensities.push_back(all_ay_ints);

1126

1127

for (Size int_it = 0; int_it != pre_path_intensities.size(); ++int_it)

1128

{

1129

if (pre_path_intensities[int_it] < MIN_DECIMAL_VALUE && suf_path_intensities[int_it] < MIN_DECIMAL_VALUE)

1130

{

1131

continue;

1132

}

1133

DoubleReal prefix_weight = prefixes[i].getMonoWeight(Residue::BIon);

1134

DoubleReal suffix_weight = suffixes[i].getMonoWeight(Residue::YIon);

1135

IsotopeDistribution prefix_id;

1136

1137

if (int_it != pre_path_intensities.size() - 1)

1138

{

1139

prefix_id = prefixes[i].getFormula(Residue::BIon).getIsotopeDistribution(max_isotope);

1140

}

1141

else

1142

{

1143

prefix_id = prefixes[i].getFormula(Residue::AIon).getIsotopeDistribution(max_isotope);

1144

}

1145

IsotopeDistribution suffix_id = suffixes[i].getFormula(Residue::YIon).getIsotopeDistribution(max_isotope);

1146

1147

for (UInt z = 1; z <= charge && z <= max_fragment_charge; ++z)

1148

{

1149

if (prefix_intensities[int_it][i][z - 1] > MIN_DECIMAL_VALUE)

1150

{

1151

vector<RichPeak1D> b_loss_peaks;

1152

DoubleReal avail_bb_sum_prefix(0);

1153

1154

if (int_it != pre_path_intensities.size() - 1)

1155

{

1156

avail_bb_sum_prefix = getAvailableBackboneCharge_(prefixes[i], Residue::BIon, z);

1157

if (avail_bb_sum_prefix <= charge_remote_threshold)

1158

{

1159

if (i == 1) // second BB position

1160

{

1161

b2_ion_losses_cr_.getIons(b_loss_peaks, prefixes[i], prefix_intensities[int_it][i][z - 1] * pre_path_intensities[int_it]);

1162

}

1163

else

1164

{

1165

b_ion_losses_cr_.getIons(b_loss_peaks, prefixes[i], prefix_intensities[int_it][i][z - 1] * pre_path_intensities[int_it]);

1166

}

1167

}

1168

else

1169

{

1170

if (i == 1)

1171

{

1172

b2_ion_losses_cd_.getIons(b_loss_peaks, prefixes[i], prefix_intensities[int_it][i][z - 1] * pre_path_intensities[int_it]);

1173

}

1174

else

1175

{

1176

b_ion_losses_cd_.getIons(b_loss_peaks, prefixes[i], prefix_intensities[int_it][i][z - 1] * pre_path_intensities[int_it]);

1177

}

1178

}

1179

}

1180

else

1181

{

1182

avail_bb_sum_prefix = getAvailableBackboneCharge_(prefixes[i], Residue::AIon, z);

1183

if (avail_bb_sum_prefix <= charge_remote_threshold)

1184

{

1185

a_ion_losses_cr_.getIons(b_loss_peaks, prefixes[i], prefix_intensities[int_it][i][z - 1] * pre_path_intensities[int_it]);

1186

}

1187

else

1188

{

1189

a_ion_losses_cd_.getIons(b_loss_peaks, prefixes[i], prefix_intensities[int_it][i][z - 1] * pre_path_intensities[int_it]);

1190

}

1191

}

1192

1193

for (vector<RichPeak1D>::const_iterator it = b_loss_peaks.begin(); it != b_loss_peaks.end(); ++it)

1194

{

1195

String b_ion_name = it->getMetaValue("IonName");

1196

#ifdef INIT_CHARGE_DEBUG

1197

cerr << b_ion_name << " " << it->getMZ() << " " << it->getIntensity() << endl;

1198

#endif

1199

vector<String> split;

1200

b_ion_name.split('-', split);

1201

if (split.empty())

1202

{

1203

b_ion_name += String(i + 1);

1204

}

1205

else

1206

{

1207

b_ion_name = split[0] + String(i + 1) + "-";

1208

for (Size j = 1; j != split.size(); ++j)

1209

{

1210

b_ion_name += split[j];

1211

}

1212

}

1213

1214

b_ion_name += String(z, '+');

1215

if (int_it != pre_path_intensities.size() - 1)

1216

{

1217

addPeaks_(prefix_weight, z, it->getMZ(), it->getIntensity(), spec, prefix_id, b_ion_name);

1218

}

1219

else

1220

{

1221

addPeaks_(prefix_weight - 28.0, z, it->getMZ(), it->getIntensity(), spec, prefix_id, b_ion_name);

1222

}

1223

}

1224

}

1225

1226

if (suffix_intensities[int_it][i][z - 1] > MIN_DECIMAL_VALUE && int_it != suf_path_intensities.size() - 1) // no ay

1227

{

1228

vector<RichPeak1D> y_loss_peaks;

1229

DoubleReal avail_bb_sum_suffix = getAvailableBackboneCharge_(suffixes[i], Residue::YIon, z);

1230

if (avail_bb_sum_suffix <= charge_remote_threshold)

1231

{

1232

y_ion_losses_cr_.getIons(y_loss_peaks, suffixes[i], suffix_intensities[int_it][i][z - 1] * suf_path_intensities[int_it]);

1233

}

1234

else

1235

{

1236

y_ion_losses_cd_.getIons(y_loss_peaks, suffixes[i], suffix_intensities[int_it][i][z - 1] * suf_path_intensities[int_it]);

1237

}

1238

1239

for (vector<RichPeak1D>::const_iterator it = y_loss_peaks.begin(); it != y_loss_peaks.end(); ++it)

1240

{

1241

String y_ion_name = it->getMetaValue("IonName");

1242

#ifdef INIT_CHARGE_DEBUG

1243

cerr << y_ion_name << " " << it->getMZ() << " " << it->getIntensity() << endl;

1244

#endif

1245

vector<String> split;

1246

y_ion_name.split('-', split);

1247

if (split.empty())

1248

{

1249

y_ion_name += String(peptide.size() - i - 1);

1250

}

1251

else

1252

{

1253

y_ion_name = split[0] + String(peptide.size() - i - 1) + "-";

1254

for (Size j = 1; j != split.size(); ++j)

1255

{

1256

y_ion_name += split[j];

1257

}

1258

}

1259

1260

y_ion_name += String(z, '+');

1261

addPeaks_(suffix_weight, z, it->getMZ(), it->getIntensity(), spec, suffix_id, y_ion_name);

1262

}

1263

}

1264

}

1265

}

1266

}

1267

1268

hmm_.disableTransitions();

1269

1270

// precursor intensities

1271

vector<RichPeak1D> pre_peaks;

1272

if (is_charge_remote)

1273

{

1274

precursor_model_cr_.getIons(pre_peaks, peptide, precursor_init);

1275

}

1276

else

1277

{

1278

precursor_model_cd_.getIons(pre_peaks, peptide, precursor_init);

1279

}

1280

1281

DoubleReal weight = peptide.getMonoWeight();

1282

IsotopeDistribution id(max_isotope);

1283

id.estimateFromPeptideWeight(weight);

1284

for (vector<RichPeak1D>::const_iterator it = pre_peaks.begin(); it != pre_peaks.end(); ++it)

1285

{

1286

addPeaks_(weight, charge, it->getMZ(), it->getIntensity(), spec, id, it->getMetaValue("IonName"));

1287

}

1288

1289

// now build the spectrum with the peaks

1290

DoubleReal intensity_max(0);

1291

for (Map<DoubleReal, vector<RichPeak1D> >::ConstIterator it = peaks_.begin(); it != peaks_.end(); ++it)

1292

{

1293

if (it->second.size() == 1 && it->second.begin()->getIntensity() != 0)

1294

{

1295

spec.push_back(*it->second.begin());

1296

if (intensity_max < spec.back().getIntensity())

1297

{

1298

intensity_max = spec.back().getIntensity();

1299

}

1300

}

1301

else

1302

{

1303

RichPeak1D p;

1304

DoubleReal int_sum(0);

1305

p = *it->second.begin();

1306

set<String> names;

1307

for (vector<RichPeak1D>::const_iterator pit = it->second.begin(); pit != it->second.end(); ++pit)

1308

{

1309

int_sum += pit->getIntensity();

1310

if (String(pit->getMetaValue("IonName")) != "")

1311

{

1312

names.insert(pit->getMetaValue("IonName"));

1313

}

1314

}

1315

1316

String name;

1317

for (set<String>::const_iterator nit = names.begin(); nit != names.end(); ++nit)

1318

{

1319

name += *nit + "/";

1320

}

1321

p.setMetaValue("IonName", name);

1322

p.setIntensity(int_sum);

1323

spec.push_back(p);

1324

if (intensity_max < int_sum)

1325

{

1326

intensity_max = int_sum;

1327

}

1328

}

1329

}

1330

1331

spec.sortByPosition();

1332

1333

RichPeakSpectrum new_spec;

1334

// add up peak intensities which are close together

1335

vector<RichPeak1D> close_peaks;

1336

for (RichPeakSpectrum::ConstIterator it = spec.begin(); it != spec.end(); ++it)

1337

{

1338

// empty

1339

if (close_peaks.empty())

1340

{

1341

close_peaks.push_back(*it);

1342

continue;

1343

}

1344

1345

// include peak in cluster

1346

if (it->getPosition() - close_peaks.begin()->getMZ() < 0.1)

1347

{

1348

close_peaks.push_back(*it);

1349

continue;

1350

}

1351

else

1352

{

1353

// move cluster to

1354

if (close_peaks.size() > 1)

1355

{

1356

DoubleReal mz(0), intensity(0);

1357

String name;

1358

for (vector<RichPeak1D>::const_iterator pit = close_peaks.begin(); pit != close_peaks.end(); ++pit)

1359

{

1360

mz += pit->getMZ();

1361

intensity += pit->getIntensity();

1362

String p_name = pit->getMetaValue("IonName");

1363

if (p_name != "")

1364

{

1365

name += p_name;

1366

if (pit != close_peaks.end() - 1)

1367

{

1368

name += "/";

1369

}

1370

}

1371

}

1372

RichPeak1D peak;

1373

peak.setPosition(mz / (DoubleReal)close_peaks.size());

1374

peak.setIntensity(intensity);

1375

peak.setMetaValue("IonName", name);

1376

new_spec.push_back(peak);

1377

1378

// clear the actual cluster and actual peak

1379

close_peaks.clear();

1380

close_peaks.push_back(*it);

1381

continue;

1382

}

1383

// cluster has only one peak, add it to the spec and start new cluster

1384

new_spec.push_back(*close_peaks.begin());

1385

close_peaks.clear();

1386

close_peaks.push_back(*it);

1387

}

1388

}

1389

spec = new_spec;

1390

1391

1392

DoubleReal min_y_int((DoubleReal)param_.getValue("min_y_ion_intensity"));

1393

DoubleReal min_b_int((DoubleReal)param_.getValue("min_b_ion_intensity"));

1394

DoubleReal min_a_int((DoubleReal)param_.getValue("min_a_ion_intensity"));

1395

DoubleReal min_y_loss_int((DoubleReal)param_.getValue("min_y_loss_intensity"));

1396

DoubleReal min_b_loss_int((DoubleReal)param_.getValue("min_b_loss_intensity"));

1397

1398

1399

for (RichPeakSpectrum::Iterator it = spec.begin(); it != spec.end(); ++it)

1400

{

1401

it->setIntensity(it->getIntensity() / intensity_max);

1402

1403

String ion_name(it->getMetaValue("IonName"));

1404

if (ion_name != "")

1405

{

1406

if (ion_name.hasSubstring("y") && (charge > 2 || !ion_name.hasSubstring("++")))

1407

{

1408

if (ion_name.hasSubstring("H2O") || ion_name.hasSubstring("NH3"))

1409

{

1410

if (it->getIntensity() < min_y_loss_int)

1411

{

1412

it->setIntensity(min_y_loss_int);

1413

}

1414

}

1415

else

1416

{

1417

if (it->getIntensity() < min_y_int)

1418

{

1419

it->setIntensity(min_y_int);

1420

}

1421

}

1422

}

1423

1424

if (ion_name.hasSubstring("b") && (charge > 2 || !ion_name.hasSubstring("++")))

1425

{

1426

1427

if (ion_name.hasSubstring("H2O") || ion_name.hasSubstring("NH3"))

1428

{

1429

if (it->getIntensity() < min_b_loss_int)

1430

{

1431

it->setIntensity(min_b_loss_int);

1432

}

1433

}

1434

else

1435

{

1436

if (it->getIntensity() < min_b_int)

1437

{

1438

it->setIntensity(min_b_int);

1439

}

1440

}

1441

}

1442

1443

if (ion_name.hasSubstring("a") && (charge > 2 || !ion_name.hasSubstring("++")))

1444

{

1445

if (it->getIntensity() < min_a_int)

1446

{

1447

it->setIntensity(min_a_int);

1448

}

1449

}

1450

}

1451

}

1452

1453

return;

1454

}

1455

1456

DoubleReal PILISModel::getAvailableBackboneCharge_(const AASequence& ion, Residue::ResidueType res_type, int charge)

1457

{

1458

DoubleReal bb_sum(0);

1459

vector<DoubleReal> bb_charges, sc_charges;

1460

prot_dist_.getProtonDistribution(bb_charges, sc_charges, ion, charge, res_type);

1461

1462

for (vector<DoubleReal>::const_iterator it = bb_charges.begin(); it != bb_charges.end(); ++it)

1463

{

1464

bb_sum += *it;

1465

}

1466

1467

// activation of protons sitting at lysine and histidine side chains

1468

DoubleReal side_chain_activation(param_.getValue("side_chain_activation"));

1469

for (Size i = 0; i != ion.size(); ++i)

1470

{

1471

if (ion[i].getOneLetterCode() != "R")

1472

{

1473

bb_sum += side_chain_activation * sc_charges[i];

1474

}

1475

}

1476

1477

if (bb_sum > 1)

1478

{

1479

bb_sum = 1;

1480

}

1481

1482

if (bb_sum < (DoubleReal)param_.getValue("charge_directed_threshold") * charge)

1483

{

1484

bb_sum = (DoubleReal)param_.getValue("charge_directed_threshold") * charge;

1485

}

1486

return bb_sum;

1487

}

1488

1489

1490

bool PILISModel::getInitialTransitionProbabilities_(std::vector<DoubleReal>& bb_init,

1491

std::vector<DoubleReal>& cr_init,

1492

std::vector<DoubleReal>& sc_init,

1493

DoubleReal& precursor_init,

1494

const vector<DoubleReal>& bb_charges,

1495

const vector<DoubleReal>& sc_charges,

1496

const AASequence& peptide)

1497

{

1498

bool is_charge_remote(false);

1499

1500

DoubleReal bb_sum(0), bb_sum_orig(0);

1501

for (vector<DoubleReal>::const_iterator it = bb_charges.begin(); it != bb_charges.end(); ++it)

1502

{

1503

bb_sum += *it;

1504

}

1505

1506

if (bb_sum > 1)

1507

{

1508

bb_sum = 1;

1509

}

1510

1511

#ifdef INIT_CHARGE_DEBUG

1512

cerr << "bb_sum=" << bb_sum << endl;

1513

#endif

1514

bb_sum_orig = bb_sum;

1515

1516

if (bb_sum < (DoubleReal)param_.getValue("charge_remote_threshold"))

1517

{

1518

is_charge_remote = true;

1519

}

1520

1521

if (bb_sum < (DoubleReal)param_.getValue("charge_directed_threshold"))

1522

{

1523

bb_sum = (DoubleReal)param_.getValue("charge_directed_threshold");

1524

}

1525

1526

// side-chain activiation

1527

DoubleReal side_chain_activation(param_.getValue("side_chain_activation"));

1528

for (Size i = 0; i != peptide.size(); ++i)

1529

{

1530

if (peptide[i].getOneLetterCode() != "R")

1531

{

1532

bb_sum += side_chain_activation * sc_charges[i];

1533

}

1534

}

1535

1536

if (bb_sum > 1)

1537

{

1538

bb_sum = 1;

1539

}

1540

1541

#ifdef INIT_CHARGE_DEBUG

1542

cerr << "bb_sum after side-chain-activiation=" << bb_sum << endl;

1543

#endif

1544

1545

vector<DoubleReal> bb_charges_all = bb_charges;

1546

sort(bb_charges_all.begin(), bb_charges_all.end());

1547

DoubleReal bb_charges_median = bb_charges_all[(Size)(bb_charges_all.size()/2.0)];

1548

#ifdef INIT_CHARGE_DEBUG

1549

cerr << "bb_charges_median=" << bb_charges_median << endl;

1550

#endif

1551

1552

DoubleReal min_enhancement_factor = param_.getValue("min_enhancement_factor");

1553

DoubleReal blocker_sum(1.0);

1554

for (Size i = 0; i != peptide.size() - 1; ++i)

1555

{

1556

DoubleReal bb_enhance_factor(max(min_enhancement_factor, sqrt(bb_charges[i+1] / bb_charges_median)));

1557

#ifdef INIT_CHARGE_DEBUG

1558

cerr << "bb_enhance_factor=" << bb_enhance_factor << endl;

1559

#endif

1560

1561

if (sc_charges[i] != 0)

1562

{

1563

blocker_sum += 10.0 * sc_charges[i];

1564

}

1565

1566

bb_init.push_back(/*bb_charges[i+1] **/ bb_sum * bb_enhance_factor / blocker_sum);

1567

String aa(peptide[i].getOneLetterCode());

1568

if ((aa == "K" || aa == "R" || aa == "H"))

1569

{

1570

sc_init.push_back(sc_charges[i] /* * bb_charges_median*/);

1571

}

1572

else

1573

{

1574

sc_init.push_back(0.0);

1575

}

1576

1577

if (is_charge_remote && (aa == "D" || aa == "E" || i == peptide.size() - 2 || i == peptide.size() - 3))

1578

{

1579

cr_init.push_back(((1 - bb_sum_orig) /** bb_charges_median*/));

1580

}

1581

else

1582

{

1583

cr_init.push_back(0.0);

1584

}

1585

}

1586

precursor_init = (1 - bb_sum) /** bb_charges_median*//*bb_avg*/ / 10.0;

1587

1588

// normalize the initial probability values

1589

DoubleReal init_prob_sum(0);

1590

for (Size i = 0; i != bb_init.size(); ++i)

1591

{

1592

init_prob_sum += bb_init[i] + sc_init[i] + cr_init[i];

1593

}

1594

init_prob_sum += precursor_init;

1595

1596

for (Size i = 0; i != bb_init.size(); ++i)

1597

{

1598

bb_init[i] /= init_prob_sum;

1599

sc_init[i] /= init_prob_sum;

1600

cr_init[i] /= init_prob_sum;

1601

}

1602

precursor_init /= init_prob_sum;

1603

1604

return is_charge_remote;

1605

}

1606

1607

void PILISModel::addPeaks_(DoubleReal mz, int charge, DoubleReal offset, DoubleReal intensity, RichPeakSpectrum& /*spectrum*/, const IsotopeDistribution& id, const String& name)

1608

{

1609

if (intensity < MIN_DECIMAL_VALUE)

1610

{

1611

return;

1612

}

1613

static RichPeak1D p;

1614

UInt i = 0;

1615

for (IsotopeDistribution::ConstIterator it = id.begin(); it != id.end(); ++it, ++i)

1616

{

1617

DoubleReal pos = (mz + i + charge + offset) / (DoubleReal)charge;

1618

p.setPosition(pos);

1619

if (it == id.begin())

1620

{

1621

p.setMetaValue("IonName", String(name.c_str()));

1622

}

1623

1624

if (pos >= (DoubleReal)param_.getValue("lower_mz") && pos <= (DoubleReal)param_.getValue("upper_mz"))

1625

{

1626

p.setIntensity(intensity * it->second);

1627

peaks_[p.getMZ()].push_back(p);

1628

}

1629

1630

if (it == id.begin())

1631

{

1632

p.setMetaValue("IonName", String(""));

1633

}

1634

}

1635

1636

return;

1637

}

1638

1639

void PILISModel::parseHMMModel_(const TextFile::ConstIterator& begin, const TextFile::ConstIterator& end, HiddenMarkovModel& hmm, Param& param)

1640

{

1641

if (begin == end)

1642

{

1643

return;

1644

}

1645

1646

//Size num_syn(0);

1647

for (TextFile::ConstIterator it = begin; it != end; ++it)

1648

{

1649

String line = *it;

1650

// comment?

1651

//cerr << line << endl;

1652

if (line[0] == '#')

1653

{

1654

continue;

1655

}

1656

1657

vector<String> split;

1658

line.split(' ', split, true);

1659

1660

if ( !split.empty() )

1661

{

1662

String id = split[0];

1663

1664

if (id == "State")

1665

{

1666

bool hidden(true);

1667

if (split.size() != 2 && split[2] == "false")

1668

{

1669

hidden = false;

1670

}

1671

hmm.addNewState(new HMMState(split[1], hidden));

1672

//cerr << "added new state: '" << split[1] << "', " << hidden << endl;

1673

continue;

1674

}

1675

1676

if (id == "Synonym")

1677

{

1678

//++num_syn;

1679

//hmm.addSynonymTransition(split[1], split[2], split[3], split[4]);

1680

hmm.addSynonymTransition(split[3], split[4], split[1], split[2]);

1681

continue;

1682

}

1683

1684

if (id == "Transition")

1685

{

1686

hmm.setTransitionProbability(split[1], split[2], split[3].toFloat());

1687

continue;

1688

}

1689

1690

if (id == "Parameter")

1691

{

1692

// Parameter charge_remote_threshold 0.3 float

1693

if (split[split.size() - 1] == "float")

1694

{

1695

param.setValue(split[1], split[2].toDouble());

1696

}

1697

else if (split[split.size() - 1] == "int")

1698

{

1699

param.setValue(split[1], split[2].toInt());

1700

}

1701

else if (split[split.size() - 1] == "string_list")

1702

{

1703

String tmp_list;

1704

for (Size i = 2; i < split.size() - 1; ++i)

1705

{

1706

tmp_list += split[i];

1707

}

1708

param.setValue(split[1], StringList::create(tmp_list));

1709

}

1710

else if (split[split.size() - 1] == "string")

1711

{

1712

param.setValue(split[1], split[2]);

1713

}

1714

else

1715

{

1716

throw Exception::InvalidParameter(__FILE__, __LINE__, __PRETTY_FUNCTION__, line);

1717

}

1718

}

1719

}

1720

}

1721

//hmm_.disableTransitions();

1722

//hmm.buildSynonyms();

1723

hmm.disableTransitions();

1724

1725

//cerr << hmm_.getNumberOfStates() << endl;

1726

1727

return;

1728

}

1729

1730

void PILISModel::writeParameters_(ostream& os, const Param& param)

1731

{

1732

for (Param::ParamIterator it = param.begin(); it != param.end(); ++it)

1733

{

1734

os << "Parameter ";

1735

if (it->value.valueType() == DataValue::DOUBLE_VALUE)

1736

{

1737

os << it->name << " \"" << it->value << "\" float" << endl;

1738

continue;

1739

}

1740

if (it->value.valueType() == DataValue::INT_VALUE)

1741

{

1742

os << it->name << " \"" << it->value << "\" int" << endl;

1743

continue;

1744

}

1745

if (it->value.valueType() == DataValue::STRING_LIST)

1746

{

1747

StringList value = (StringList)it->value;

1748

String tmp;

1749

tmp.concatenate(value.begin(), value.end(), ",");

1750

os << it->name << " \"" << tmp << "\" string_list" << endl;

1751

continue;

1752

}

1753

if (it->value.valueType() == DataValue::STRING_VALUE)

1754

{

1755

os << it->name << " \"" << it->value << "\" string" << endl;

1756

continue;

1757

}

1758

throw Exception::InvalidParameter(__FILE__, __LINE__, __PRETTY_FUNCTION__, "Value-type of parameter " + it->name + " not implemented!");

1759

}

1760

return;

1761

}

1762

1763

void PILISModel::updateMembers_()

1764

{

1765

DoubleReal pseudo_counts = (DoubleReal)param_.getValue("pseudo_counts");

1766

hmm_.setPseudoCounts(pseudo_counts);

1767

1768

// pass parameters to precursor model

1769

Param precursor_param_cr(precursor_model_cr_.getParameters());

1770

precursor_param_cr.setValue("pseudo_counts", pseudo_counts);

1771

precursor_model_cr_.setParameters(precursor_param_cr);

1772

1773

Param precursor_param_cd = precursor_model_cd_.getParameters();

1774

precursor_param_cd.setValue("pseudo_counts", pseudo_counts);

1775

precursor_model_cd_.setParameters(precursor_param_cd);

1776

1777

Param b_ion_losses_param_cr = b_ion_losses_cr_.getParameters();

1778

b_ion_losses_param_cr.setValue("pseudo_counts", pseudo_counts);

1779

b_ion_losses_cr_.setParameters(b_ion_losses_param_cr);

1780

1781

Param b_ion_losses_param_cd = b_ion_losses_cd_.getParameters();

1782

b_ion_losses_param_cd.setValue("pseudo_counts", pseudo_counts);

1783

b_ion_losses_cd_.setParameters(b_ion_losses_param_cd);

1784

1785

1786

Param b2_ion_losses_param_cr = b2_ion_losses_cr_.getParameters();

1787

b2_ion_losses_param_cr.setValue("pseudo_counts", pseudo_counts);

1788

b2_ion_losses_cr_.setParameters(b2_ion_losses_param_cr);

1789

1790

Param b2_ion_losses_param_cd = b2_ion_losses_cd_.getParameters();

1791

b2_ion_losses_param_cd.setValue("pseudo_counts", pseudo_counts);

1792

b2_ion_losses_cd_.setParameters(b2_ion_losses_param_cd);

1793

1794

1795

Param y_ion_losses_param_cr = y_ion_losses_cr_.getParameters();

1796

y_ion_losses_param_cr.setValue("pseudo_counts", pseudo_counts);

1797

y_ion_losses_cr_.setParameters(y_ion_losses_param_cr);

1798

1799

Param y_ion_losses_param_cd = y_ion_losses_cd_.getParameters();

1800

y_ion_losses_param_cd.setValue("pseudo_counts", pseudo_counts);

1801

y_ion_losses_cd_.setParameters(y_ion_losses_param_cd);

1802

1803

1804

Param a_ion_losses_param_cr = a_ion_losses_cr_.getParameters();

1805

a_ion_losses_param_cr.setValue("pseudo_counts", pseudo_counts);

1806

a_ion_losses_cr_.setParameters(a_ion_losses_param_cr);

1807

1808

Param a_ion_losses_param_cd = a_ion_losses_cd_.getParameters();

1809

a_ion_losses_param_cd.setValue("pseudo_counts", pseudo_counts);

1810

a_ion_losses_cd_.setParameters(a_ion_losses_param_cd);

1811

}

1812

} // namespace OpenMS

1813

1814

Older »