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

« back to all changes in this revision

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

This program was created at: Tue Sep 8 21:05:23 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 <string>

#include <iostream>

#include <fstream>

#include <math.h>

#include <cmath>

#include <stdlib.h>

#include <time.h>

#include <stdio.h>

#include <unistd.h>

#include <algorithm>

#include "split.h"

#include "gpatInfo.hpp"

struct options{

std::string file;

double afDiff;

}globalOpts;

struct iHSdat{

std::string seqid;

std::string start;

std::string F1 ;

std::string F2 ;

double af ;

double ehhR ;

double ehhA ;

double iHS ;

double niHS ;

};

using namespace std;

static const char *optString = "hf:s:";

void printHelp(void){

cerr << endl << endl;

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

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

cerr << " normalizes iHS or XP-EHH scores " << endl;

cerr << "Output : normalize-iHS adds one additional column to input (normalized score)." << endl;

cerr << "INFO: usage: normalizeHS -s 0.01 -f input.txt " << endl;

cerr << endl;

cerr << "INFO: required: -f -- Output from iHS or XPEHH " << endl;

cerr << "INFO: optional: -s -- Max AF diff for window [0.01]" << endl;

cerr << endl;

printVersion();

}

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

int parseOpts(int argc, char** argv)

{

int opt = 0;

opt = getopt(argc, argv, optString);

100

while(opt != -1){

101

switch(opt){

102

case 's':

103

{

104

string op = optarg;

105

globalOpts.afDiff = atof(op.c_str());

106

break;

107

}

108

case 'h':

109

{

110

printHelp();

111

exit(1);

112

break;

113

}

114

115

case 'f':

116

{

117

globalOpts.file = optarg;

118

break;

119

}

120

case '?':

121

{

122

break;

123

}

124

}

125

opt = getopt( argc, argv, optString );

126

}

127

return 1;

128

}

129

130

bool sortAF(iHSdat * L, iHSdat * R){

131

if(L->af < R->af){

132

return true;

133

}

134

return false;

135

}

136

137

138

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

139

140

Function input : vector of doubles

141

142

Function does : calculates the var

143

144

Function returns: double

145

146

147

148

double var(vector<double> & data, double mu){

149

double variance = 0;

150

151

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

152

variance += pow((*it) - mu,2);

153

}

154

155

return variance / (data.size() - 1);

156

}

157

158

159

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

160

161

Function input : vector of doubles

162

163

Function does : computes the mean

164

165

Function returns: the mean

166

167

168

169

170

double windowAvg(std::vector<double> & rangeData){

171

172

long double n = 0;

173

long double s = 0;

174

175

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

176

s += *it;

177

n += 1;

178

}

179

180

181

return (s/n);

182

}

183

184

185

186

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

187

188

Function input : vector of iHS data

189

190

Function does : normalizes

191

192

Function returns: nothing

193

194

195

196

void normalize(std::vector<iHSdat *> & data, int * pos){

197

198

std::vector<double> windat;

199

200

int start = *pos;

201

int end = *pos;

202

203

while((abs(data[start]->af - data[end]->af ) < globalOpts.afDiff)

204

&& end < data.size() -1 ){

205

end += 1;

206

}

207

208

for(int i = start; i <= end; i++){

209

windat.push_back(data[i]->iHS);

210

}

211

212

double avg = windowAvg(windat);

213

double sd = sqrt(var(windat, avg));

214

215

std::cerr << "start: " << data[start]->af << " "

216

<< "end: " << data[end]->af << " "

217

<< "n iHS scores: " << windat.size() << " "

218

<< "mean: " << avg << " "

219

<< "sd: " << sd << std::endl;

220

221

for(int i = start; i <= end; i++){

222

data[i]->niHS = (data[i]->iHS - avg) / (sd);

223

}

224

225

*pos = end;

226

227

}

228

229

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

230

231

Comments:

232

233

234

int main( int argc, char** argv)

235

{

236

globalOpts.afDiff = 0.01;

237

int parse = parseOpts(argc, argv);

238

239

if(globalOpts.file.empty()){

240

std::cerr << "FATAL: no file" << std::endl;

241

exit(1);

242

}

243

244

std::vector<iHSdat *> data;

245

246

string line;

247

ifstream myfile (globalOpts.file);

248

if (myfile.is_open())

249

{

250

while ( getline (myfile,line) ){

251

vector<string> lineDat = split(line, '\t');

252

253

iHSdat * tp = new iHSdat;

254

tp->seqid = lineDat[0];

255

tp->start = lineDat[1];

256

tp->af = atof(lineDat[2].c_str());

257

tp->ehhR = atof(lineDat[3].c_str());

258

tp->ehhA = atof(lineDat[4].c_str());

259

tp->iHS = atof(lineDat[5].c_str());

260

tp->F1 = lineDat[6].c_str();

261

tp->F2 = lineDat[7].c_str();

262

tp->niHS = 0;

263

264

data.push_back(tp);

265

266

}

267

268

myfile.close();

269

}

270

else{

271

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

272

exit(1);

273

}

274

275

276

std::cerr << "INFO: sorting " << data.size() << " scores by AF" << std::endl;

277

278

sort(data.begin(), data.end(), sortAF);

279

280

std::cerr << "INFO: finished sorting" << std::endl;

281

282

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

283

normalize(data, &i);

284

}

285

286

287

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

288

std::cout << data[i]->seqid << "\t"

289

<< data[i]->start << "\t"

290

<< data[i]->af << "\t"

291

<< data[i]->ehhR << "\t"

292

<< data[i]->ehhA << "\t"

293

<< data[i]->iHS << "\t"

294

<< data[i]->niHS << "\t"

295

<< data[i]->F1 << "\t"

296

<< data[i]->F2 << std::endl;

297

298

}

299

300

301

return 0;

302

}

Older »