~ubuntu-branches/debian/sid/libvcflib/sid

« back to all changes in this revision

Viewing changes to src/bFst.cpp

Committer: Package Import Robot
Author(s): Andreas Tille
Date: 2016-09-16 15:52:29 UTC
Revision ID: package-import@ubuntu.com-20160916155229-24mxrntfylvsshsg

Tags: upstream-1.0.0~rc1+dfsg

Import upstream version 1.0.0~rc1+dfsg

files added:

.gitignore

.gitmodules

.travis.yml

1KG_Index

1KG_Index/1KG_IDs.txt

1KG_Index/1kG_ethnic_index.py

1KG_Index/README

LICENSE

Makefile

Makefile.old

README.md

bin/plotBfst.R

bin/plotHapLrt.R

bin/plotHaplotypes.R

bin/plotPfst.R

bin/plotSmoothed.R

bin/plotWCfst.R

bin/plotXPEHH.R

examples

examples/612.bFst_SatApr12_12_29_17.png

examples/612.counts_FriApr11_13_16_42.png

examples/612.nocounts.smoothed_TueApr15_13_25_53.png

examples/612.nocounts_FriApr11_13_16_38.png

examples/612.wcfst.txt_FriApr11_13_17_02.png

examples/example-ehh.png

examples/headCrest.haps.txt_ThuMay29_15_08_53.pdf

examples/headCrest.haps.txt_ThuMay29_15_08_53.png

examples/phasing-diff-res.png

examples/phasing-results.png

examples/scaffold612.d-stat.10kb.txt_abba-baba_WedOct29_10_17_10.png

examples/t.hapass.txt_TueJun24_11_52_08.png

examples/xp-phased.txt

examples/xp-phased.txt_TueApr29_13_20_47.png

examples/xp-unphased.txt

examples/xp-unphased.txt_TueApr29_13_20_58.png

examples/xpehh_WedApr23_11_30_10.png

fastahack

filevercmp

fsom

intervaltree

logos

logos/websiteLogo.pdf

logos/websiteLogo.png

multichoose

samples

samples/sample.vcf

samples/scaffold612.phased.vcf

samples/scaffold612.vcf

smithwaterman

src/BedReader.h

src/LD.cpp

src/Variant.cpp

src/Variant.h

src/abba-baba.cpp

src/bFst.cpp

src/cdflib.cpp

src/cdflib.hpp

src/convert.h

src/dumpContigsFromHeader.cpp

src/genotypeSummary.cpp

src/gl-XPEHH.cpp

src/gpatInfo.hpp

src/hapLrt.cpp

src/iHS.cpp

src/join.h

src/meltEHH.cpp

src/mt19937ar.h

src/normalize-iHS.cpp

src/pFst.cpp

src/pdflib.cpp

src/pdflib.hpp

src/permuteGPAT++.cpp

src/permuteGPATsmoother.cpp

src/permuteRegions.cpp

src/permuteSmooth.cpp

src/plotHaps.cpp

src/popStats.cpp

src/rnglib.cpp

src/rnglib.hpp

src/segmentFst.cpp

src/segmentIhs.cpp

src/sequenceDiversity.cpp

src/smoother.cpp

src/split.cpp

src/split.h

src/splitUniqStarts.cpp

src/ssw.c

src/ssw.h

src/ssw_cpp.cpp

src/ssw_cpp.h

src/var.cpp

src/var.hpp

src/vcf2dag.cpp

src/vcf2fasta.cpp

src/vcf2tsv.cpp

src/vcfToHap.cpp

src/vcfaddinfo.cpp

src/vcfafpath.cpp

src/vcfallelicprimitives.cpp

src/vcfaltcount.cpp

src/vcfannotate.cpp

src/vcfannotategenotypes.cpp

src/vcfbreakmulti.cpp

src/vcfcat.cpp

src/vcfcheck.cpp

src/vcfclassify.cpp

src/vcfcleancomplex.cpp

src/vcfcombine.cpp

src/vcfcommonsamples.cpp

src/vcfcountalleles.cpp

src/vcfcreatemulti.cpp

src/vcfdistance.cpp

src/vcfecho.cpp

src/vcfentropy.cpp

src/vcfevenregions.cpp

src/vcffilter.cpp

src/vcffixup.cpp

src/vcfflatten.cpp

src/vcfgeno2alleles.cpp

src/vcfgeno2haplo.cpp

src/vcfgenosamplenames.cpp

src/vcfgenosummarize.cpp

src/vcfgenotypecompare.cpp

src/vcfgenotypes.cpp

src/vcfglbound.cpp

src/vcfglxgt.cpp

src/vcfhetcount.cpp

src/vcfhethomratio.cpp

src/vcfindex.cpp

src/vcfinfo2qual.cpp

src/vcfinfosummarize.cpp

src/vcfintersect.cpp

src/vcfkeepgeno.cpp

src/vcfkeepinfo.cpp

src/vcfkeepsamples.cpp

src/vcfleftalign.cpp

src/vcflength.cpp

src/vcfnumalt.cpp

src/vcfoverlay.cpp

src/vcfparsealts.cpp

src/vcfprimers.cpp

src/vcfqual2info.cpp

src/vcfrandom.cpp

src/vcfrandomsample.cpp

src/vcfremap.cpp

src/vcfremoveaberrantgenotypes.cpp

src/vcfremovesamples.cpp

src/vcfroc.cpp

src/vcfsample2info.cpp

src/vcfsamplediff.cpp

src/vcfsamplenames.cpp

src/vcfsamplestats.cpp

src/vcfsitesummarize.cpp

src/vcfsom.cpp

src/vcfstats.cpp

src/vcfstreamsort.cpp

src/vcfuniq.cpp

src/vcfuniqalleles.cpp

src/wcFst.cpp

src/xpEHH.cpp

tabixpp

tests

tests/lib

tests/lib/Local

tests/lib/Local/vcflib

tests/lib/Local/vcflib/Test.pm

tests/vcfdistance.t

Show diffs side-by-side

added added

removed removed

src/bFst.cpp

#include "Variant.h"

#include "split.h"

#include "pdflib.hpp"

#include <string>

#include <iostream>

#include <math.h>

#include <cmath>

#include <stdlib.h>

#include <time.h>

#include <stdio.h>

#include <getopt.h>

#include "gpatInfo.hpp"

using namespace std;

using namespace vcflib;

struct pop{

double nalt ;

double nref ;

double af ;

double nhomr;

double nhoma;

double nhet ;

double ngeno;

double fis ;

vector<int> questionable;

vector<int> geno_index ;

vector< vector< double > > unphred_p;

};

double unphred(string phred){

double unphred = atof(phred.c_str());

unphred = unphred / -10;

return unphred;

}

void initPop(pop & population){

population.nalt = 0;

population.nref = 0;

population.af = 0;

population.nhomr = 0;

population.nhoma = 0;

population.nhet = 0;

population.ngeno = 0;

population.fis = 0;

}

void loadPop( vector< map< string, vector<string> > >& group, pop & population){

vector< map< string, vector<string> > >::iterator targ_it = group.begin();

int index = 0;

for(; targ_it != group.end(); targ_it++){

string genotype = (*targ_it)["GT"].front();

vector<double> phreds;

phreds.push_back( unphred((*targ_it)["PL"][0]));

phreds.push_back( unphred((*targ_it)["PL"][1]));

phreds.push_back( unphred((*targ_it)["PL"][2]));

double scaled ;

double norm = log(exp(phreds[0]) + exp(phreds[1]) + exp(phreds[2]));

population.unphred_p.push_back(phreds);

while(1){

if(genotype == "0/0"){

population.ngeno += 1;

population.nhomr += 1;

population.nref += 2;

population.geno_index.push_back(0);

scaled = exp(phreds[0] - norm);

break;

}

if(genotype == "0/1"){

population.ngeno += 1;

population.nhet += 1;

population.nref += 1;

population.nalt += 1;

population.geno_index.push_back(1);

scaled = exp(phreds[1] - norm);

break;

}

if(genotype == "1/1"){

population.ngeno += 1;

population.nhoma += 1;

population.nalt += 2;

population.geno_index.push_back(2);

scaled = exp(phreds[2] - norm);

break;

}

if(genotype == "0|0"){

100

population.ngeno += 1;

101

population.nhomr += 1;

102

population.nref += 2;

103

population.geno_index.push_back(0);

104

scaled = exp(phreds[0] - norm);

105

break;

106

}

107

if(genotype == "0|1"){

108

population.ngeno += 1;

109

population.nhet += 1;

110

population.nref += 1;

111

population.nalt += 1;

112

population.geno_index.push_back(1);

113

scaled = exp(phreds[1] - norm);

114

break;

115

}

116

if(genotype == "1|1"){

117

population.ngeno += 1;

118

population.nhoma += 1;

119

population.nalt += 2;

120

population.geno_index.push_back(2);

121

scaled = exp(phreds[2] - norm);

122

break;

123

}

124

cerr << "FATAL: unknown genotype" << endl;

125

exit(1);

126

}

127

128

if(scaled < 0.75){

129

population.questionable.push_back(index);

130

}

131

index += 1;

132

}

133

134

if(population.nalt == 0 && population.nref == 0){

135

population.af = -1;

136

}

137

else{

138

139

population.af = (population.nalt / (population.nref + population.nalt));

140

141

if(population.nhet > 0){

142

population.fis = ( 1 - ((population.nhet/population.ngeno) / (2*population.af*(1 - population.af))));

143

}

144

else{

145

population.fis = 1;

146

}

147

if(population.fis < 0){

148

population.fis = 0.00001;

149

}

150

}

151

}

152

153

double bound(double v){

154

if(v <= 0.00001){

155

return 0.00001;

156

}

157

if(v >= 0.99999){

158

return 0.99999;

159

}

160

return v;

161

}

162

163

164

void phardy(vector<double>& results, double af, double fis){

165

166

double p0 = pow((1 - af),2) + ((1 - af)*af*fis);

167

double p1 = 2*(1 - af)*af*(1 - fis);

168

double p2 = pow(af,2) + ((1 - af)*af*fis);

169

170

results.push_back(p0);

171

results.push_back(p1);

172

results.push_back(p2);

173

174

}

175

176

double likelihood(pop & population, double af, double fis){

177

178

af = bound(af);

179

fis = bound(fis);

180

181

double het;

182

double all;

183

double ref;

184

double alt;

185

double loglikelihood = 0;

186

187

vector<double> genotypeProbs;

188

189

phardy(genotypeProbs, af, fis);

190

191

vector<int>::iterator it = population.geno_index.begin();

192

193

int geno_indx = 0;

194

195

for(; it != population.geno_index.end(); it++){

196

197

double aa = population.unphred_p[geno_indx][0] + log(genotypeProbs[0]);

198

double ab = population.unphred_p[geno_indx][1] + log(genotypeProbs[1]);

199

double bb = population.unphred_p[geno_indx][2] + log(genotypeProbs[2]);

200

double norm = exp(aa) + exp(ab) + exp(bb);

201

202

double prop = population.unphred_p[geno_indx][*it] + log(genotypeProbs[*it]);

203

204

loglikelihood += (prop - norm);

205

206

geno_indx++;

207

}

208

209

return loglikelihood;

210

211

}

212

213

double FullProb(pop & target, pop & back, vector<double>& p)

214

{

215

216

// parameters targetAf backgroundAf targetFis backgroundFis totalAf fst

217

218

double alpha = ( (1-p[5])/p[5] ) * p[4];

219

double beta = ( (1-p[5])/p[5] ) * (1 - p[4]);

220

221

222

double afprior = log( r8_normal_pdf (p[6], 0.1, p[4]));

223

double afpriorT = log( r8_normal_pdf (target.af, 0.05, p[0]));

224

double afpriorB = log( r8_normal_pdf (back.af, 0.05, p[1]));

225

226

if(std::isinf(afprior) || std::isnan(afprior)){

227

return -100000;

228

}

229

230

double ptaf = log( r8_beta_pdf(alpha, beta, p[0]) );

231

double pbaf = log( r8_beta_pdf(alpha, beta, p[1]) );

232

233

if( std::isinf(ptaf) || std::isnan(ptaf) || std::isinf(pbaf) || std::isnan(pbaf) ){

234

return -100000;

235

}

236

237

238

double llt = likelihood(target, p[0], p[2]);

239

double llb = likelihood(back, p[1], p[3]);

240

double full = llt + llb + ptaf + pbaf + afprior + afpriorT + afpriorB;

241

242

243

return full;

244

245

}

246

247

248

249

void updateParameters(pop & target, pop & background, vector<double>& parameters, int pindx){

250

251

// parameters targetAf backgroundAf targetFis backgroundFis totalAf fst

252

253

double origpar = parameters[pindx];

254

255

double accept = ((double)rand() / (double)(RAND_MAX));

256

double up = ((double)rand() / (double)(RAND_MAX))/10 - 0.05;

257

double updatep = parameters[pindx] + up;

258

259

// cerr << accept << "\t" << up << endl;

260

261

if(updatep >= 1 || updatep <= 0){

262

return;

263

}

264

265

double llB = FullProb(target, background, parameters);

266

parameters[pindx] = updatep;

267

double llT = FullProb(target, background, parameters);

268

269

if((llT - llB) > accept){

270

return;

271

}

272

else{

273

parameters[pindx] = origpar;

274

}

275

}

276

277

void updateGenotypes(pop & target, pop & background, vector<double>& parameters, int gindex, int tbindex){

278

279

// tbindex indicates if the subroutine will update the target or background genotype;

280

281

double accept = ((double)rand() / (double)(RAND_MAX));

282

int newGindex = rand() % 3;

283

284

//cerr << newGindex << endl;

285

//cerr << "gindex " << gindex << endl;

286

//cerr << "gsize t:" << target.geno_index.size() << endl;

287

//cerr << "gsize b:" << background.geno_index.size() << endl;

288

289

290

291

int oldtindex = target.geno_index[gindex] ;

292

int oldbindex = background.geno_index[gindex] ;

293

294

double llB = FullProb(target, background, parameters);

295

296

if(tbindex == 0){

297

//udate target

298

target.geno_index[gindex] = newGindex;

299

}

300

else{

301

// update background

302

background.geno_index[gindex] = newGindex;

303

}

304

double llT = FullProb(target, background, parameters);

305

306

if((llT - llB) > accept){

307

return;

308

}

309

else{

310

if(tbindex == 0){

311

target.geno_index[gindex] = oldtindex;

312

}

313

else{

314

target.geno_index[gindex] = oldbindex;

315

}

316

}

317

}

318

319

320

int cmp(const void *x, const void *y)

321

{

322

double xx = *(double*)x, yy = *(double*)y;

323

if (xx < yy) return -1;

324

if (xx > yy) return 1;

325

return 0;

326

}

327

328

329

void loadIndices(map<int, int> & index, string set){

330

331

vector<string> indviduals = split(set, ",");

332

333

vector<string>::iterator it = indviduals.begin();

334

335

for(; it != indviduals.end(); it++){

336

index[ atoi( (*it).c_str() ) ] = 1;

337

}

338

}

339

340

341

int main(int argc, char** argv) {

342

343

// set the random seed for MCMC

344

345

srand((unsigned)time(NULL));

346

347

// the filename

348

349

string filename = "NA";

350

351

// using vcflib; thanks to Erik Garrison

352

353

VariantCallFile variantFile ;

354

355

// zero based index for the target and background indivudals

356

357

map<int, int> it, ib;

358

359

// deltaaf is the difference of allele frequency we bother to look at

360

361

string deltaaf ;

362

double daf = -1;

363

364

const struct option longopts[] =

365

{

366

{"version" , 0, 0, 'v'},

367

{"help" , 0, 0, 'h'},

368

{"file" , 1, 0, 'f'},

369

{"target" , 1, 0, 't'},

370

{"background", 1, 0, 'b'},

371

{"deltaaf" , 1, 0, 'd'},

372

{0,0,0,0}

373

};

374

375

int index;

376

int iarg = 0;

377

378

while(iarg != -1)

379

{

380

iarg = getopt_long(argc, argv, "d:t:b:f:hv", longopts, &index);

381

382

switch (iarg)

383

{

384

case 0:

385

break;

386

case 'h':

387

cerr << endl;

388

389

cerr << "INFO: help: " << endl << endl;

390

391

cerr << " bFst is a Bayesian approach to Fst. Importantly bFst account for genotype uncertainty in the model using genotype likelihoods." << endl;

392

cerr << " For a more detailed description see: Holsinger et al. Molecular Ecology Vol 11, issue 7 2002. The likelihood function has been " << endl;

393

cerr << " modified to use genotype likelihoods provided by variant callers. There are five free parameters estimated in the model: each " << endl;

394

cerr << " subpopulation's allele frequency and Fis (fixation index, within each subpopulation), a free parameter for the total population\'s " << endl;

395

cerr << " allele frequency, and Fst. " << endl << endl;

396

397

cerr << "Output : 11 columns : " << endl;

398

cerr << " 1. Seqid " << endl;

399

cerr << " 2. Position " << endl;

400

cerr << " 3. Observed allele frequency in target. " << endl;

401

cerr << " 4. Estimated allele frequency in target. " << endl;

402

cerr << " 5. Observed allele frequency in background. " << endl;

403

cerr << " 6. Estimated allele frequency in background. " << endl;

404

cerr << " 7. Observed allele frequency combined. " << endl;

405

cerr << " 8. Estimated allele frequency in combined. " << endl;

406

cerr << " 9. ML estimate of Fst (mean) " << endl;

407

cerr << " 10. Lower bound of the 95% credible interval " << endl;

408

cerr << " 11. Upper bound of the 95% credible interval " << endl << endl;

409

410

411

cerr << "INFO: usage: bFst --target 0,1,2,3,4,5,6,7 --background 11,12,13,16,17,19,22 --file my.vcf --deltaaf 0.1" << endl;

412

cerr << endl;

413

cerr << "INFO: required: t,target -- a zero bases comma separated list of target individuals corrisponding to VCF columns" << endl;

414

cerr << "INFO: required: b,background -- a zero bases comma separated list of background individuals corrisponding to VCF columns" << endl;

415

cerr << "INFO: required: f,file a -- a proper formatted VCF file. the FORMAT field MUST contain \"PL\"" << endl;

416

cerr << "INFO: required: d,deltaaf -- skip sites were the difference in allele frequency is less than deltaaf" << endl;

417

cerr << endl;

418

printVersion();

419

cerr << endl << endl;

420

return 0;

421

422

case 'v':

423

printVersion();

424

return 0;

425

426

case 't':

427

loadIndices(ib, optarg);

428

cerr << "INFO: There are " << ib.size() << " individuals in the target" << endl;

429

break;

430

431

case 'b':

432

loadIndices(it, optarg);

433

cerr << "INFO: There are " << it.size() << " individuals in the background" << endl;

434

break;

435

436

case 'f':

437

cerr << "INFO: File: " << optarg << endl;

438

filename = optarg;

439

break;

440

441

case 'd':

442

cerr << "INFO: difference in allele frequency : " << optarg << endl;

443

deltaaf = optarg;

444

daf = atof(deltaaf.c_str());

445

break;

446

default:

447

break;

448

cerr << endl;

449

cerr << "FATAL: unknown command line option " << optarg << endl << endl ;

450

cerr << "INFO: please use bFst --help " << endl;

451

cerr << endl;

452

return(1);

453

}

454

455

}

456

457

if(daf == -1){

458

cerr << endl;

459

cerr << "FATAL: did not specify deltaaf" << endl;

460

cerr << "INFO: please use bFst --help " << endl;

461

cerr << endl;

462

return(1);

463

}

464

465

if(filename == "NA"){

466

cerr << endl;

467

cerr << "FATAL: did not specify VCF file" << endl;

468

cerr << "INFO: please use bFst --help " << endl;

469

cerr << endl;

470

return(1);

471

}

472

473

variantFile.open(filename);

474

475

476

if (!variantFile.is_open()) {

477

cerr << endl;

478

cerr << "FATAL: could not open VCF file" << endl;

479

cerr << "INFO: please use bFst --help" << endl;

480

cerr << endl;

481

return(1);

482

}

483

if(it.size() < 2){

484

cerr << endl;

485

cerr << "FATAL: target not specified or less than two indviduals" << endl;

486

cerr << "INFO: please use bFst --help " << endl;

487

cerr << endl;

488

}

489

if(ib.size() < 2){

490

cerr << endl;

491

cerr << "FATAL: target not specified or less than two indviduals"<< endl;

492

cerr << "INFO: please use bFst --help " << endl;

493

cerr << endl;

494

}

495

496

Variant var(variantFile);

497

498

vector<string> samples = variantFile.sampleNames;

499

int nsamples = samples.size();

500

501

while (variantFile.getNextVariant(var)) {

502

503

// biallelic sites naturally

504

505

if(var.alt.size() > 1){

506

continue;

507

}

508

509

510

vector < map< string, vector<string> > > target, background, total;

511

512

int index = 0;

513

514

for(int nsamp = 0; nsamp < nsamples; nsamp++){

515

516

map<string, vector<string> > sample = var.samples[ samples[nsamp]];

517

518

if(sample["GT"].front() != "./."){

519

if(it.find(index) != it.end() ){

520

target.push_back(sample);

521

total.push_back(sample);

522

523

}

524

if(ib.find(index) != ib.end()){

525

background.push_back(sample);

526

total.push_back(sample);

527

}

528

}

529

530

index += 1;

531

}

532

533

if(target.size() < 2 || background.size() < 2 ){

534

continue;

535

}

536

537

pop popt, popb, popTotal;

538

539

initPop(popt);

540

initPop(popb);

541

initPop(popTotal);

542

543

loadPop(target, popt);

544

loadPop(background, popb);

545

loadPop(total, popTotal);

546

547

if(popt.af == -1 || popb.af == -1){

548

continue;

549

}

550

if(popt.af == 1 && popb.af == 1){

551

continue;

552

}

553

if(popt.af == 0 && popb.af == 0){

554

continue;

555

}

556

557

double afdiff = abs(popt.af - popb.af);

558

559

if(afdiff < daf){

560

continue;

561

}

562

563

564

cerr << "INFO: target has " << popt.questionable.size() << " questionable genotypes " << endl;

565

cerr << "INFO: background has " << popb.questionable.size() << " questionable genotypes " << endl;

566

567

// Parameters- targetAf backgroundAf targetFis backgroundFis totalAf fst

568

vector<double> parameters;

569

parameters.push_back(popt.af);

570

parameters.push_back(popb.af);

571

parameters.push_back(popt.fis);

572

parameters.push_back(popb.fis);

573

parameters.push_back(popTotal.af);

574

parameters.push_back(0.1);

575

parameters.push_back(popTotal.af);

576

577

double sums [6] = {0};

578

double fsts [10000] ;

579

580

for(int i = 0; i < 15000; i++){

581

582

// update each of j parameters

583

584

for(int j = 0; j < 6; j++ ){

585

586

updateParameters(popt, popb, parameters, j);

587

if(i > 4999){

588

sums[j] += parameters[j];

589

}

590

}

591

if(i > 4999){

592

fsts[i - 5000] = parameters[5];

593

}

594

for(vector<int>::iterator itt = popt.questionable.begin(); itt != popt.questionable.end(); itt++){

595

updateGenotypes(popt, popb, parameters, (*itt), 0);

596

597

}

598

for(vector<int>::iterator itb = popb.questionable.begin(); itb != popb.questionable.end(); itb++){

599

updateGenotypes(popt, popb, parameters, (*itb) , 1);

600

}

601

}

602

603

qsort (fsts, sizeof(fsts)/sizeof(fsts[0]), sizeof(fsts[0]), cmp );

604

605

double lcredint = fsts[500];

606

double hcredint = fsts[9500];

607

608

cout << var.sequenceName << "\t" << var.position

609

<< "\t" << popt.af

610

<< "\t" << sums[0]/10000

611

<< "\t" << popb.af

612

<< "\t" << sums[1]/10000

613

<< "\t" << popTotal.af

614

<< "\t" << sums[4]/10000

615

<< "\t" << sums[5]/10000

616

<< "\t" << lcredint

617

<< "\t" << hcredint

618

<< endl;

619

}

620

return 0;

621

}

Older »