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

« back to all changes in this revision

Viewing changes to src/meltEHH.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/meltEHH.cpp

#include "Variant.h"

#include "split.h"

#include "cdflib.hpp"

#include "pdflib.hpp"

#include "var.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"

// maaas speed

#include <omp.h>

// print lock

omp_lock_t lock;

struct opts{

int threads ;

int pos ;

std::string filename ;

std::string mapFile ;

std::string seqid ;

std::string geneticMapFile ;

std::string type ;

std::string region ;

std::map<int, double> geneticMap;

double af ;

}globalOpts;

using namespace std;

using namespace vcflib;

void printHelp(void){

cerr << endl << endl;

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

cerr << "INFO: description:" << endl;

cerr << " meltEHH provides the data to plot EHH curves. " << endl;

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

cerr << " 1. seqid " << endl;

cerr << " 2. position " << endl;

cerr << " 3. EHH " << endl;

cerr << " 4. ref or alt [0 == ref] " << endl;

cerr << "Usage:" << endl;

cerr << " meltEHH --target 0,1,2,3,4,5,6,7 --pos 10 --file my.phased.vcf \\" << endl;

cerr << " --region chr1:1-1000 > STDOUT 2> STDERR " << endl << endl;

cerr << "Params:" << endl;

cerr << " required: t,target <STRING> A zero base comma separated list of target" << endl;

cerr << " individuals corresponding to VCF columns " << endl;

cerr << " required: r,region <STRING> A tabix compliant genomic range " << endl;

cerr << " format: \"seqid:start-end\" or \"seqid\" " << endl;

cerr << " required: f,file <STRING> Proper formatted and phased VCF. " << endl;

cerr << " required: y,type <STRING> Genotype likelihood format: GT,PL,GL,GP " << endl;

cerr << " required: p,position <INT> Variant position to melt. " << endl;

cerr << " optional: a,af <DOUBLE> Alternative alleles with frequencies less " << endl;

cerr << " than [0.05] are skipped. " << endl;

cerr << endl;

printVersion();

exit(1);

}

bool gDist(int start, int end, double * gd){

if(globalOpts.geneticMap.find(start) == globalOpts.geneticMap.end()){

return false;

}

if(globalOpts.geneticMap.find(end) == globalOpts.geneticMap.end()){

return false;

}

*gd = abs(globalOpts.geneticMap[start] - globalOpts.geneticMap[end]);

return true;

}

void loadGeneticMap(int start, int end){

if(globalOpts.geneticMapFile.empty()){

std::cerr << "WARNING: No genetic map." << std::endl;

std::cerr << "WARNING: A constant genetic distance is being used: 0.001." << std::endl;

return;

}

ifstream featureFile (globalOpts.geneticMapFile.c_str());

string line;

100

101

int lastpos = 0;

102

double lastvalue = 0;

103

104

if(featureFile.is_open()){

105

106

while(getline(featureFile, line)){

107

108

vector<string> region = split(line, "\t");

109

110

if(region.front() != globalOpts.seqid){

111

std::cerr << "WARNING: seqid MisMatch: " << region.front() << " " << globalOpts.seqid << std::endl;

112

continue;

113

}

114

115

int pos = atoi(region[3].c_str()) ;

116

double cm = atof(region[2].c_str()) ;

117

118

if(lastpos == 0 && start > pos){

119

lastpos = pos;

120

continue;

121

}

122

123

int diff = abs(pos - lastpos);

124

double vdiff = abs(lastvalue - cm );

125

double chunk = vdiff/double(diff);

126

127

double running = lastvalue;

128

129

for(int i = lastpos; i < pos; i++){

130

globalOpts.geneticMap[i] = running;

131

running += chunk;

132

}

133

134

if(pos > end){

135

break;

136

}

137

138

139

lastpos = pos;

140

lastvalue = cm;

141

}

142

}

143

144

featureFile.close();

145

146

if(globalOpts.geneticMap.size() < 1){

147

std::cerr << "FATAL: Problem loading genetic map" << std::endl;

148

exit(1);

149

}

150

}

151

152

153

void clearHaplotypes(string **haplotypes, int ntarget){

154

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

155

haplotypes[i][0].clear();

156

haplotypes[i][1].clear();

157

}

158

}

159

160

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

161

162

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

163

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

164

165

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

166

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

167

}

168

}

169

170

void countHaps(int nhaps, map<string, int> & targetH,

171

string **haplotypes, int start, int end){

172

173

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

174

175

std::string h1 = haplotypes[i][0].substr(start, (end - start)) ;

176

std::string h2 = haplotypes[i][1].substr(start, (end - start)) ;

177

178

if(targetH.find(h1) == targetH.end()){

179

targetH[h1] = 1;

180

}

181

else{

182

targetH[h1]++;

183

}

184

if(targetH.find(h2) == targetH.end()){

185

targetH[h2] = 1;

186

}

187

else{

188

targetH[h2]++;

189

}

190

}

191

}

192

193

void computeNs(map<string, int> & targetH, int start,

194

int end, double * sumT, char ref, bool dir){

195

196

for( map<string, int>::iterator th = targetH.begin();

197

th != targetH.end(); th++){

198

199

if(th->second < 2){

200

continue;

201

}

202

203

204

// end is extending ; check first base

205

if(dir){

206

if( th->first[0] == ref){

207

208

// std::cerr << "count dat: " << th->first << " " << th->second << " " << ref << " " << dir << endl;

209

210

211

*sumT += r8_choose(th->second, 2);

212

}

213

}

214

215

// start is extending ; check last base

216

else{

217

218

int last = th->first.size() -1;

219

if( th->first[last] == ref ){

220

// std::cerr << "count dat:" << th->first << " " << th->second << " " << ref << " " << dir << endl;

221

222

223

*sumT += r8_choose(th->second, 2);

224

}

225

}

226

}

227

}

228

229

bool calcEhh(string **haplotypes, int start,

230

int end, char ref, int nhaps,

231

double * ehh, double div, bool dir){

232

233

double sum = 0 ;

234

map<string , int> refH;

235

236

countHaps(nhaps, refH, haplotypes, start, end);

237

computeNs(refH, start, end, &sum, ref, dir );

238

239

double internalEHH = sum / (r8_choose(div, 2));

240

241

if(internalEHH > 1){

242

std::cerr << "FATAL: internal error." << std::endl;

243

exit(1);

244

}

245

246

*ehh = internalEHH;

247

248

return true;

249

}

250

251

int integrate(string ** haplotypes,

252

vector<long int> & pos,

253

bool direction,

254

int maxl,

255

int snp,

256

char ref,

257

int nhaps,

258

double * iHH,

259

double denom ){

260

261

double ehh = 1;

262

263

int start = snp;

264

int end = snp;

265

266

// controls the substring madness

267

if(!direction){

268

start += 1;

269

end += 1;

270

}

271

272

while(ehh > 0.01){

273

if(direction){

274

end += 1;

275

}

276

else{

277

start -= 1;

278

}

279

if(start < 0){

280

return 1;

281

}

282

if(end > maxl){

283

return 1;

284

}

285

double ehhRT = 0;

286

if(!calcEhh(haplotypes,

287

start, end,

288

ref, nhaps,

289

&ehhRT, denom,

290

direction)){

291

return 1;

292

}

293

294

if(ehhRT <= 0.01){

295

return 0;

296

}

297

298

double delta_gDist = 0.001;

299

300

if(direction){

301

gDist(pos[end-1], pos[end], &delta_gDist);

302

}

303

else{

304

gDist(pos[start + 1], pos[start], &delta_gDist);

305

}

306

*iHH += ((ehh + ehhRT)/2)*delta_gDist;

307

308

if(direction){

309

std::cout << pos[end] << "\t" << ehh << "\t" << ref << "\t" << direction << std::endl;

310

}

311

else{

312

std::cout << pos[start] << "\t" << ehh << "\t" << ref << "\t" << direction << std::endl;

313

}

314

315

316

ehh = ehhRT;

317

318

}

319

320

return 10;

321

}

322

323

void calc(string ** haplotypes, int nhaps,

324

vector<double> & afs, vector<long int> & pos,

325

vector<int> & target, vector<int> & background, string seqid){

326

327

int maxl = haplotypes[0][0].length();

328

329

330

for(int snp = 0; snp < maxl; snp++){

331

332

if(pos[snp] != globalOpts.pos ){

333

continue;

334

}

335

336

337

double ihhR = 0;

338

double ihhA = 0;

339

340

map<string , int> refH;

341

342

countHaps(nhaps, refH, haplotypes, snp, snp+1);

343

344

345

double denomP1 = double(refH["0"]);

346

double denomP2 = double(refH["1"]);

347

348

int refFail = 0;

349

int altFail = 0;

350

351

std::cout << pos[snp] << "\t" << "1" << "\t" << "0" << "\t" << "0" << std::endl;

352

353

354

355

refFail += integrate(haplotypes, pos, true, maxl, snp, '0', nhaps, &ihhR, denomP1);

356

357

refFail += integrate(haplotypes, pos, false, maxl, snp, '0', nhaps, &ihhR, denomP1);

358

359

altFail += integrate(haplotypes, pos, true, maxl, snp, '1', nhaps, &ihhA, denomP2);

360

361

altFail += integrate(haplotypes, pos, false, maxl, snp, '1', nhaps, &ihhA, denomP2);

362

363

}

364

}

365

366

void loadPhased(string **haplotypes, genotype * pop, int ntarget){

367

368

int indIndex = 0;

369

370

for(vector<string>::iterator ind = pop->gts.begin(); ind != pop->gts.end(); ind++){

371

string g = (*ind);

372

vector< string > gs = split(g, "|");

373

haplotypes[indIndex][0].append(gs[0]);

374

haplotypes[indIndex][1].append(gs[1]);

375

indIndex += 1;

376

}

377

}

378

379

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

380

381

globalOpts.threads = 1 ;

382

globalOpts.af = 0.05;

383

384

// zero based index for the target and background indivudals

385

386

map<int, int> it, ib;

387

388

const struct option longopts[] =

389

{

390

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

391

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

392

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

393

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

394

{"region" , 1, 0, 'r'},

395

{"gen" , 1, 0, 'g'},

396

{"type" , 1, 0, 'y'},

397

{"threads" , 1, 0, 'x'},

398

{"af" , 1, 0, 'a'},

399

{"pos" , 1, 0, 'p'},

400

{0,0,0,0}

401

};

402

403

int findex;

404

int iarg=0;

405

406

while(iarg != -1)

407

{

408

iarg = getopt_long(argc, argv, "a:x:g:y:r:d:t:b:f:p:hv", longopts, &findex);

409

410

switch (iarg)

411

{

412

case 'p':

413

{

414

globalOpts.pos = atoi(optarg);

415

break;

416

}

417

418

case 'a':

419

{

420

globalOpts.af = atof(optarg);

421

break;

422

}

423

case 'x':

424

{

425

globalOpts.threads = atoi(optarg);

426

break;

427

}

428

case 'g':

429

{

430

globalOpts.geneticMapFile = optarg;

431

break;

432

}

433

case 'h':

434

{

435

printHelp();

436

break;

437

}

438

case 'v':

439

{

440

printVersion();

441

break;

442

}

443

case 'y':

444

{

445

globalOpts.type = optarg;

446

break;

447

}

448

case 't':

449

{

450

loadIndices(it, optarg);

451

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

452

cerr << "INFO: target ids: " << optarg << endl;

453

break;

454

}

455

case 'f':

456

{

457

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

458

globalOpts.filename = optarg;

459

break;

460

}

461

case 'r':

462

{

463

cerr << "INFO: set seqid region to : " << optarg << endl;

464

globalOpts.region = optarg;

465

break;

466

default:

467

break;

468

}

469

}

470

}

471

472

omp_set_num_threads(globalOpts.threads);

473

474

map<string, int> okayGenotypeLikelihoods;

475

okayGenotypeLikelihoods["PL"] = 1;

476

okayGenotypeLikelihoods["GL"] = 1;

477

okayGenotypeLikelihoods["GP"] = 1;

478

okayGenotypeLikelihoods["GT"] = 1;

479

480

481

// add an option for dumping

482

483

// for(std::map<int, double>::iterator gm = geneticMap.begin(); gm != geneticMap.end(); gm++){

484

// cerr << "pos: " << gm->first << " cm: " << gm->second << endl;

485

// }

486

487

if(globalOpts.type.empty()){

488

cerr << "FATAL: failed to specify genotype likelihood format : PL or GL" << endl;

489

printHelp();

490

exit(1);

491

}

492

if(okayGenotypeLikelihoods.find(globalOpts.type) == okayGenotypeLikelihoods.end()){

493

cerr << "FATAL: genotype likelihood is incorrectly formatted, only use: PL or GL" << endl;

494

printHelp();

495

exit(1);

496

}

497

498

if(globalOpts.filename.empty()){

499

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

500

printHelp();

501

exit(1);

502

}

503

504

if(it.size() < 2){

505

cerr << "FATAL: target option is required -- or -- less than two individuals in target\n";

506

printHelp();

507

exit(1);

508

}

509

510

// using vcflib; thanksErik

511

512

VariantCallFile variantFile;

513

514

variantFile.open(globalOpts.filename);

515

516

if(globalOpts.region.empty()){

517

cerr << "FATAL: region required" << endl;

518

exit(1);

519

}

520

if(! variantFile.setRegion(globalOpts.region)){

521

cerr <<"FATAL: unable to set region" << endl;

522

exit(1);

523

}

524

525

if (!variantFile.is_open()) {

526

exit(1);

527

}

528

529

Variant var( variantFile );

530

vector<int> target_h, background_h;

531

532

int index = 0;

533

int indexi = 0;

534

535

536

vector<string> samples = variantFile.sampleNames;

537

int nsamples = samples.size();

538

539

for(vector<string>::iterator samp = samples.begin(); samp != samples.end(); samp++){

540

541

string sampleName = (*samp);

542

543

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

544

target_h.push_back(indexi);

545

indexi++;

546

}

547

index++;

548

}

549

550

551

vector<long int> positions;

552

553

vector<double> afs;

554

555

string **haplotypes = new string*[target_h.size()];

556

for (int i = 0; i < target_h.size(); i++) {

557

haplotypes[i] = new string[2];

558

}

559

560

561

while (variantFile.getNextVariant(var)) {

562

563

globalOpts.seqid = var.sequenceName;

564

565

if(!var.isPhased()){

566

cerr << "FATAL: Found an unphased variant. All genotypes must be phased!" << endl;

567

exit(1);

568

}

569

570

if(var.alleles.size() > 2){

571

continue;

572

}

573

574

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

575

576

int sindex = 0;

577

578

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

579

580

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

581

582

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

583

target.push_back(sample);

584

}

585

sindex += 1;

586

}

587

588

genotype * populationTarget ;

589

590

if(globalOpts.type == "PL"){

591

populationTarget = new pl();

592

}

593

if(globalOpts.type == "GL"){

594

populationTarget = new gl();

595

}

596

if(globalOpts.type == "GP"){

597

populationTarget = new gp();

598

}

599

if(globalOpts.type == "GT"){

600

populationTarget = new gt();

601

}

602

603

populationTarget->loadPop(target, var.sequenceName, var.position);

604

605

if(populationTarget->af <= globalOpts.af

606

|| populationTarget->nref < 2

607

|| populationTarget->nalt < 2){

608

delete populationTarget;

609

continue;

610

}

611

positions.push_back(var.position);

612

afs.push_back(populationTarget->af);

613

loadPhased(haplotypes, populationTarget, populationTarget->gts.size());

614

615

populationTarget = NULL;

616

delete populationTarget;

617

}

618

619

if(!globalOpts.geneticMapFile.empty()){

620

cerr << "INFO: loading genetics map" << endl;

621

loadGeneticMap(positions.front(), positions.back());

622

cerr << "INFO: finished loading genetics map" << endl;

623

}

624

625

calc(haplotypes, target_h.size(), afs, positions,

626

target_h, background_h, globalOpts.seqid);

627

clearHaplotypes(haplotypes, target_h.size());

628

629

exit(0);

630

631

}

Older »