Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[feat] Merge ibdgl method branch #2

Merged
merged 2 commits into from
Dec 12, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
54 changes: 46 additions & 8 deletions argStruct.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -44,6 +44,7 @@ argStruct *argStruct_get(int argc, char **argv) {
// -doDxy <int> : estimate Dxy
// -doPhylo <int> : do neighbor-joining
// -doDist <int> : estimate pairwise distance matrix
// -doIbd <int> : detect IBD segments

// TODO idea multilayer argument reading using LUTs
// e.g. `-doAMOVA` in lut1 -> "args->doAMOVA" in lut2 ->integer in lut3
Expand All @@ -66,7 +67,9 @@ argStruct *argStruct_get(int argc, char **argv) {

} else if (strcasecmp("-doDist", arv) == 0) {
args->doDist = atoi(val);
}
} else if (strcasecmp("-doIbd", arv) == 0) {
args->doIbd = atoi(val);
}

// ################################
// # INPUT FILES [--in-XXX] #
Expand Down Expand Up @@ -345,6 +348,37 @@ argStruct *argStruct_get(int argc, char **argv) {

else if (strcasecmp("--minInd", arv) == 0)
args->minInd = atoi(val);
else if (strcasecmp("--min-af", arv) == 0){
args->min_af = atof(val);
}


else if (strcasecmp("--ibdseq-errorprop", arv) == 0){
args->ibdseq_errorprop = atof(val);
}

else if (strcasecmp("--ibdseq-errormax", arv) == 0){
args->ibdseq_errormax = atof(val);
}

else if (strcasecmp("--ibdseq-minalleles", arv) == 0){
args->ibdseq_minalleles = atoi(val);
if(args->ibdseq_minalleles<2){
ERROR("--ibdseq-minalleles is set to %d. Allowed range of values is 2<value<N");
}
}
else if (strcasecmp("--ibdseq-ibdlod", arv) == 0){
args->ibdseq_ibdlod= atof(val);
if(args->ibdseq_ibdlod<0){
ERROR("--ibdseq-ibdlod must be a positive number.");
}
}
else if (strcasecmp("--ibdseq-ibdtrim", arv) == 0){
args->ibdseq_ibdtrim= atof(val);
if(args->ibdseq_ibdtrim<0){
ERROR("--ibdseq-ibdtrim must be a nonnegative number.");
}
}

else if (strcasecmp("--em-tole", arv) == 0)
args->tole = atof(val);
Expand Down Expand Up @@ -457,10 +491,12 @@ void argStruct::check_arg_dependencies() {
// 2: estimate distance matrix from genotypes
// 3: read distance matrix from file

if (0 == doDist) {
fprintf(stderr, "\n\n-doDist 0: Nothing to do; will exit.\n\n");
exit(0);
} else if (1 == doDist) {
//TODO
// if (0 == doDist) {
// fprintf(stderr, "\n\n-doDist 0: Nothing to do; will exit.\n\n");
// exit(0);
// } else
if (1 == doDist) {
LOG("-doDist 1, will estimate distance matrix from genotype likelihoods.");

IO::requireArgFile(in_vcf_fn, "--in-vcf", "-doDist 1");
Expand All @@ -479,9 +515,11 @@ void argStruct::check_arg_dependencies() {
if (0 != printDistanceMatrix) {
ERROR("-doDist 3 cannot be used with --printDistanceMatrix.\n");
}
} else {
ERROR("-doDist %d is not a valid option.", doDist);
}
}
//TODO
// } else {
// ERROR("-doDist %d is not a valid option.", doDist);
// }

//----------------------------------------------------------------------------------//
// -doAMOVA
Expand Down
19 changes: 19 additions & 0 deletions argStruct.h
Original file line number Diff line number Diff line change
Expand Up @@ -168,6 +168,7 @@ struct argStruct {
int doDxy = 0;
int doPhylo = 0;
int doDist = 0;
int doIbd = 0;

int printAmovaTable = 0;
int printDistanceMatrix = 0;
Expand All @@ -184,6 +185,24 @@ struct argStruct {
double tole = -1;
int maxEmIter = -1;


// ------------------------------------------------
// VCF filters
double min_af = 0.0;
// ------------------------------------------------



// ------------------------------------------------
// IBDSEQ method implementation
int ibdseq_minalleles=2;
double ibdseq_errormax=0.001;
double ibdseq_errorprop=0.25;
double ibdseq_ibdlod=3.0;
double ibdseq_ibdtrim=0.0;



int nThreads = 1;
int seed = -1;
int nBootstraps = 0;
Expand Down
Loading