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

« back to all changes in this revision

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

This program was created at: Fri Apr 17 14:59:53 2015

This program was created by: Zev N. Kronenberg

Contact: zev.kronenber@gmail.com

Organization: Unviersity of Utah

School of Medicine

Salt Lake City, Utah

The MIT License (MIT)

Permission is hereby granted, free of charge, to any person obtaining a copy

of this software and associated documentation files (the "Software"), to deal

in the Software without restriction, including without limitation the rights

to use, copy, modify, merge, publish, distribute, sublicense, and/or sell

copies of the Software, and to permit persons to whom the Software is

furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in

all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR

IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,

FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE

AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER

LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,

OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN

THE SOFTWARE.

#include <fstream>

#include "split.h"

#include <vector>

#include <map>

#include <string>

#include <iostream>

#include <math.h>

#include <cmath>

#include <stdlib.h>

#include <time.h>

#include <stdio.h>

#include <unistd.h>

#include <stdlib.h>

#include "gpatInfo.hpp"

#include <omp.h>

// print lock

omp_lock_t lock;

struct options{

std::string file ;

std::string smoothed;

std::string format ;

double npermutation ;

double nsuc ;

int threads ;

int chrIndex ;

int nIndex ;

int valueIndex ;

}globalOpts;

struct score{

std::string seqid;

long int pos ;

double score ;

};

struct smoothedInputData{

std::string line;

double score ;

double n ;

double nPer ;

double nSuc ;

double ePv ;

};

static const char *optString = "x:y:u:f:n:s:";

using namespace std;

std::map<std::string, int> FORMATS;

//------------------------------- SUBROUTINE --------------------------------

Function input : NA

Function does : prints help

Function returns: NA

void printHelp()

100

{

101

102

cerr << endl << endl;

103

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

104

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

105

cerr << " permuteSmooth is a method for adding empirical p-values smoothed wcFst scores." << endl ;

106

cerr << endl;

107

cerr << "INFO: usage: permuteSmooth -s wcFst.smooth.txt -f wcFst.txt -n 5 -s 1 "<< endl;

108

cerr << endl;

109

cerr << "Required:" << endl;

110

cerr << " file: f -- argument: original wcFst data "<< endl;

111

cerr << " smoothed: s -- argument: smoothed wcFst data "<< endl;

112

cerr << " format: y -- argument: [swcFst, segwcFst] "<< endl;

113

cerr << "Optional:" << endl;

114

cerr << " number: n -- argument: the number of permutations to run for each value [1000]" << endl;

115

cerr << " success: u -- argument: stop permutations after \'s\' successes [1]" << endl;

116

cerr << " success: x -- argument: number of threads [1]" << endl;

117

cerr << endl;

118

cerr << "OUTPUT: permuteSmooth will append three additional columns:" << endl;

119

cerr << " 1. The number of successes " << endl;

120

cerr << " 2. The number of trials " << endl;

121

cerr << " 3. The empirical p-value " << endl;

122

cerr << endl;

123

printVersion();

124

125

}

126

127

128

//------------------------------- OPTIONS --------------------------------

129

int parseOpts(int argc, char** argv)

130

{

131

int opt = 0;

132

globalOpts.file = "NA";

133

134

globalOpts.nsuc = 1;

135

globalOpts.npermutation = 1000;

136

137

opt = getopt(argc, argv, optString);

138

while(opt != -1){

139

switch(opt){

140

case 'x':

141

{

142

globalOpts.threads = atoi(optarg);

143

}

144

case 'f':

145

{

146

globalOpts.file = optarg;

147

break;

148

}

149

case 'y':

150

{

151

globalOpts.format = optarg;

152

if(FORMATS.find(globalOpts.format) == FORMATS.end()){

153

std::cerr << "FATAL: format not supported: " << globalOpts.format;

154

std::cerr << endl;

155

printHelp();

156

exit(1);

157

}

158

if(globalOpts.format == "swcFst"){

159

globalOpts.chrIndex = 0;

160

globalOpts.nIndex = 3;

161

globalOpts.valueIndex = 4;

162

}

163

if(globalOpts.format == "segwcFst"){

164

globalOpts.chrIndex = 0;

165

globalOpts.valueIndex = 3;

166

globalOpts.nIndex = 5;

167

}

168

break;

169

}

170

case 'n':

171

{

172

globalOpts.npermutation = atof(((string)optarg).c_str());

173

cerr << "INFO: permuteGPAT++ will do N permutations: " << globalOpts.npermutation << endl;

174

break;

175

}

176

case 's':

177

{

178

globalOpts.smoothed = optarg;

179

cerr << "INFO: smoothed file: " << globalOpts.smoothed << endl;

180

break;

181

}

182

case 'u':

183

{

184

globalOpts.nsuc = atof(((string)optarg).c_str());

185

cerr << "INFO: permuteGPAT++ will stop permutations after N successes: " << globalOpts.nsuc << endl;

186

break;

187

}

188

case '?':

189

{

190

break;

191

}

192

}

193

194

opt = getopt( argc, argv, optString );

195

}

196

return 1;

197

}

198

199

200

//------------------------------- SUBROUTINE --------------------------------

201

202

Function input : data, vector to load

203

204

Function does : returns a contguous window

205

206

Function returns: bool

207

208

209

210

bool getContiguousWindow(vector<score *> & data,

211

vector<double> & load,

212

int n, int * nfail){

213

int r = rand() % data.size();

214

215

if(r+n >= data.size()){

216

*nfail+=1;

217

return false;

218

}

219

220

if(data[r]->seqid != data[r+n]->seqid){

221

*nfail+=1;

222

return false;

223

}

224

225

for(int i = r; i < r+n; i++){

226

load.push_back(data[i]->score);

227

}

228

return true;

229

}

230

231

232

//------------------------------- SUBROUTINE --------------------------------

233

234

double mean(vector<double> & data){

235

236

double sum = 0;

237

238

for(vector<double>::iterator it = data.begin(); it != data.end(); it++){

239

sum += (*it);

240

}

241

return sum / data.size();

242

}

243

244

245

246

//------------------------------- SUBROUTINE --------------------------------

247

248

Function input : score, n, data

249

250

Function does : permutes the score. requires that the window is contiguous

251

252

Function returns: NA

253

254

255

256

bool permute(double s, int n, vector<score *> & data,

257

double * nRep, double *nSuc, double * ePv){

258

259

260

*ePv = 1 / globalOpts.npermutation;

261

262

while(*nSuc < globalOpts.nsuc && *nRep < globalOpts.npermutation ){

263

*nRep += 1;

264

265

266

std::vector<double> scores;

267

268

bool getWindow = false;

269

270

int nfail = 0;

271

272

while(!getWindow){

273

if(nfail > globalOpts.npermutation){

274

return false;

275

}

276

getWindow = getContiguousWindow(data, scores, n, &nfail);

277

}

278

279

double ns = mean(scores);

280

281

if(ns > s){

282

*nSuc += 1;

283

}

284

}

285

if(*nSuc > 0){

286

*ePv = *nSuc / *nRep;

287

}

288

289

return true;

290

}

291

292

//------------------------------- MAIN --------------------------------

293

294

Comments:

295

296

297

int main( int argc, char** argv)

298

{

299

srand (time(NULL));

300

301

FORMATS["swcFst"] = 1;

302

FORMATS["segwcFst"] = 1;

303

304

globalOpts.threads = 1;

305

306

int parse = parseOpts(argc, argv);

307

308

omp_set_num_threads(globalOpts.threads);

309

310

311

if(globalOpts.file.compare("NA") == 0){

312

cerr << "FATAL: no file was provided" << endl;

313

printHelp();

314

exit(1);

315

}

316

if(globalOpts.format.empty()){

317

std::cerr << "FATAL: no format specified." << std::endl;

318

exit(1);

319

}

320

321

vector< score *> data;

322

323

ifstream wcDat (globalOpts.file.c_str());

324

325

string line;

326

327

if(wcDat.is_open()){

328

329

while(getline(wcDat, line)){

330

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

331

// will change for other output

332

double value = atof(region[4].c_str());

333

334

if(globalOpts.format == "swcFst" || globalOpts.format == "segwcFst"){

335

336

if(region.size() != 5){

337

cerr << "FATAL: wrong number of columns in wcFst input" << endl;

338

exit(1);

339

}

340

if(value < 0 ){

341

value = 0;

342

}

343

}

344

345

score * sp;

346

sp = new score;

347

sp->seqid = region[0] ;

348

sp->pos = atoi(region[1].c_str());

349

sp->score = value ;

350

351

data.push_back(sp);

352

}

353

}

354

else{

355

cerr << "FATAL: coult not open file: " << globalOpts.file << endl;

356

exit(1);

357

}

358

359

360

wcDat.close();

361

line.clear();

362

363

cerr << "INFO: raw values to permute against: " << data.size() << endl;

364

365

ifstream smoothedFile (globalOpts.smoothed.c_str());

366

367

vector<smoothedInputData*> sData;

368

369

if(smoothedFile.is_open()){

370

while(getline(smoothedFile, line)){

371

372

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

373

smoothedInputData * sp = new smoothedInputData;

374

375

sp->line = line;

376

sp->score = atof(region[globalOpts.valueIndex].c_str());

377

sp->n = atoi(region[globalOpts.nIndex].c_str());

378

sp->nPer = 0;

379

sp->nSuc = 0;

380

sp->ePv = 0;

381

382

sData.push_back(sp);

383

}

384

}

385

smoothedFile.close();

386

387

cerr << "INFO: Number of smoothed windows to permute : " << sData.size() << endl;

388

389

#pragma omp parallel for schedule(dynamic, 20)

390

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

391

bool per = permute(sData[i]->score, sData[i]->n,

392

data, &sData[i]->nPer,

393

&sData[i]->nSuc, &sData[i]->ePv);

394

395

396

if(per){

397

omp_set_lock(&lock);

398

cout << sData[i]->line

399

<< "\t" << sData[i]->nSuc

400

<< "\t" << sData[i]->nPer

401

<< "\t" << sData[i]->ePv << endl;

402

omp_unset_lock(&lock);

403

}

404

else{

405

omp_set_lock(&lock);

406

cout << sData[i]->line

407

<< "\t" << "NA"

408

<< "\t" << "NA"

409

<< "\t" << "NA" << endl;

410

omp_unset_lock(&lock);

411

}

412

}

413

414

for(vector<score*>::iterator itz = data.begin();

415

itz != data.end(); itz++){

416

delete *itz;

417

}

418

for(vector<smoothedInputData*>::iterator itz = sData.begin();

419

itz != sData.end(); itz++){

420

delete *itz;

421

}

422

423

424

return 0;

425

}

Older »