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

« back to all changes in this revision

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

using namespace std;

using namespace vcflib;

void printVersion(void){

cerr << "INFO: version 1.0.1 ; date: April 2014 ; author: Zev Kronenberg; email : zev.kronenberg@utah.edu " << endl;

exit(1);

}

void printHelp(void){

cerr << endl << endl;

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

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

cerr << " plotHaps provides the formatted output that can be used with \'bin/plotHaplotypes.R\'. " << endl;

cerr << "Output : haplotype matrix and positions" << endl << endl;

cerr << "INFO: plotHaps --target 0,1,2,3,4,5,6,7 --file my.phased.vcf.gz " << endl << endl;

cerr << "INFO: required: t,target -- argument: a zero base comma separated list of target individuals corrisponding to VCF column s " << endl;

cerr << "INFO: required: r,region -- argument: a tabix compliant genomic range : \"seqid:start-end\" or \"seqid\" " << endl;

cerr << "INFO: required: f,file -- argument: proper formatted phased VCF file " << endl;

cerr << "INFO: required: y,type -- argument: genotype likelihood format: PL,GP,GP " << endl;

cerr << endl;

printVersion();

exit(1);

}

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

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

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

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

}

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

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

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

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

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

}

void calc(string **haplotypes, int nhaps, vector<double> afs, vector<long int> pos, vector<int> & target, vector<int> & background, string seqid){

for(int snp = 0; snp < haplotypes[0][0].length(); snp++){

double ehhA = 1;

double ehhR = 1;

double iHSA = 1;

double iHSR = 1;

int start = snp;

int end = snp;

int core = snp;

while( ehhA > 0.05 && ehhR > 0.05 ) {

start -= 1;

end += 1;

if(start == -1){

break;

}

if(end == haplotypes[0][0].length() - 1){

break;

}

map<string , int> targetH;

double sumrT = 0;

double sumaT = 0;

double nrefT = 0;

double naltT = 0;

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

targetH[ haplotypes[i][0].substr(start, (end - start)) ]++;

targetH[ haplotypes[i][1].substr(start, (end - start)) ]++;

}

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

if( (*th).first.substr((end-start)/2, 1) == "1"){

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

100

naltT += th->second;

101

}

102

else{

103

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

104

nrefT += th->second;

105

}

106

}

107

108

ehhA = sumaT / (r8_choose(naltT, 2));

109

ehhR = sumrT / (r8_choose(nrefT, 2));

110

111

iHSA += ehhA;

112

iHSR += ehhR;

113

}

114

cout << seqid << "\t" << pos[snp] << "\t" << afs[snp] << "\t" << iHSA << "\t" << iHSR << "\t" << iHSA/iHSR << endl;

115

}

116

}

117

118

double EHH(string **haplotypes, int nhaps){

119

120

map<string , int> hapcounts;

121

122

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

123

hapcounts[ haplotypes[i][0] ]++;

124

hapcounts[ haplotypes[i][1] ]++;

125

}

126

127

double sum = 0;

128

double nh = 0;

129

130

for( map<string, int>::iterator it = hapcounts.begin(); it != hapcounts.end(); it++){

131

nh += it->second;

132

sum += r8_choose(it->second, 2);

133

}

134

135

double max = (sum / r8_choose(nh, 2));

136

137

return max;

138

139

}

140

141

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

142

143

int indIndex = 0;

144

145

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

146

string g = (*ind);

147

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

148

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

149

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

150

indIndex += 1;

151

}

152

}

153

154

void printHaplotypes(string **haps, vector<int> target, vector<long int> pos){

155

for(int snp = 0; snp < haps[0][1].length(); snp++){

156

cout << pos[snp] << "\t" ;

157

for(int ind = 0; ind < target.size(); ind++){

158

cout << haps[target[ind]][0].substr(snp , 1) << "\t";

159

cout << haps[target[ind]][1].substr(snp , 1) << "\t";

160

}

161

cout << endl;

162

}

163

}

164

165

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

166

167

// set the random seed for MCMC

168

169

srand((unsigned)time(NULL));

170

171

// the filename

172

173

string filename = "NA";

174

175

// set region to scaffold

176

177

string region = "NA";

178

179

// using vcflib; thanks to Erik Garrison

180

181

VariantCallFile variantFile;

182

183

// zero based index for the target and background indivudals

184

185

map<int, int> it, ib;

186

187

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

188

189

// ancestral state is set to zero by default

190

191

192

int counts = 0;

193

194

// phased

195

196

int phased = 0;

197

198

string type = "NA";

199

200

const struct option longopts[] =

201

{

202

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

203

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

204

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

205

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

206

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

207

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

208

{0,0,0,0}

209

};

210

211

int findex;

212

int iarg=0;

213

214

while(iarg != -1)

215

{

216

iarg = getopt_long(argc, argv, "y:r:t:f:hv", longopts, &findex);

217

218

switch (iarg)

219

{

220

case 'h':

221

printHelp();

222

case 'v':

223

printVersion();

224

case 'y':

225

type = optarg;

226

break;

227

case 't':

228

loadIndices(it, optarg);

229

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

230

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

231

break;

232

case 'f':

233

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

234

filename = optarg;

235

break;

236

case 'r':

237

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

238

region = optarg;

239

break;

240

default:

241

break;

242

}

243

}

244

245

map<string, int> okayGenotypeLikelihoods;

246

okayGenotypeLikelihoods["PL"] = 1;

247

okayGenotypeLikelihoods["GL"] = 1;

248

okayGenotypeLikelihoods["GP"] = 1;

249

okayGenotypeLikelihoods["GT"] = 1;

250

251

252

if(type == "NA"){

253

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

254

printHelp();

255

return 1;

256

}

257

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

258

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

259

printHelp();

260

return 1;

261

}

262

263

if(filename == "NA"){

264

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

265

printHelp();

266

return(1);

267

}

268

269

variantFile.open(filename);

270

271

if (!variantFile.is_open()) {

272

cerr << "FATAL: could not open VCF for reading" << endl;

273

printHelp();

274

return 1;

275

}

276

277

if(region != "NA"){

278

if(! variantFile.setRegion(region)){

279

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

280

return 1;

281

}

282

}

283

else{

284

cerr << "FATAL: must specifiy a region" << endl;

285

printHelp();

286

return 1;

287

}

288

289

Variant var(variantFile);

290

291

vector<string> samples = variantFile.sampleNames;

292

int nsamples = samples.size();

293

294

vector<int> target_h, background_h;

295

296

int index = 0;

297

int indexi = 0;

298

299

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

300

301

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

302

target_h.push_back(indexi);

303

indexi++;

304

}

305

index++;

306

}

307

308

vector<long int> positions;

309

310

vector<double> afs;

311

312

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

313

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

314

haplotypes[i] = new string[2];

315

}

316

317

string currentSeqid = "NA";

318

319

320

while (variantFile.getNextVariant(var)) {

321

322

if(!var.isPhased()){

323

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

324

return(1);

325

}

326

327

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

328

continue;

329

}

330

331

332

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

333

334

int sindex = 0;

335

336

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

337

338

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

339

340

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

341

target.push_back(sample);

342

}

343

sindex += 1;

344

}

345

346

genotype * populationTarget ;

347

348

if(type == "PL"){

349

populationTarget = new pl();

350

}

351

if(type == "GL"){

352

populationTarget = new gl();

353

}

354

if(type == "GP"){

355

populationTarget = new gp();

356

}

357

if(type == "GT"){

358

populationTarget = new gt();

359

}

360

361

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

362

363

positions.push_back(var.position);

364

afs.push_back(populationTarget->af);

365

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

366

}

367

368

printHaplotypes( haplotypes, target_h, positions);

369

370

return 0;

371

}

Older »