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.

Sign up for GitHub

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

Jump to bottom

Bgen madness #319

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Bgen madness #319

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
6 changes: 6 additions & 0 deletions Makefile
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,12 @@ FLAGS2 = $(CPPFLAGS) $(FLAGS) $(LDFLAGS)
CFLAGS := $(FLAGS2) $(CFLAGS)
CXXFLAGS := $(FLAGS2) $(CXXFLAGS)

#for compiling with ZSTD which is used for .bgen file format
ifeq ($(WITH_ZSTD),1)
LIBS += -lzstd
CXXFLAGS += -D__ZSTD__
endif

CSRC = $(wildcard *.c)
CXXSRC = $(wildcard *.cpp)
OBJ = $(CSRC:.c=.o) $(CXXSRC:.cpp=.o)
Expand Down
112 changes: 78 additions & 34 deletions abcAsso.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,8 @@ void abcAsso::printArg(FILE *argFile){
fprintf(argFile,"\t-yCount\t\t%s\t(File containing count phenotypes)\n",yfile2);
fprintf(argFile,"\t-yQuant\t\t%s\t(File containing phenotypes)\n",yfile3);
fprintf(argFile,"\t-cov\t\t%s\t(File containing additional covariates)\n",covfile);
fprintf(argFile,"\t-sampleFile\t\t%s\t(.sample File containing phenotypes and covariates)\n",sampleFile);

fprintf(argFile,"\t-model\t%d\n",model);
fprintf(argFile,"\t1: Additive/Log-Additive (Default)\n");
fprintf(argFile,"\t2: Dominant\n");
Expand All @@ -62,7 +64,6 @@ void abcAsso::printArg(FILE *argFile){

void abcAsso::getOptions(argStruct *arguments){


doAsso=angsd::getArg("-doAsso",doAsso,arguments);
if(doAsso==3){
fprintf(stderr,"\t-> -doAsso 3 is deprecated from version 0.615 \n");
Expand Down Expand Up @@ -103,6 +104,8 @@ void abcAsso::getOptions(argStruct *arguments){
if(yfile3!=NULL)
isQuant=1;

sampleFile=angsd::getArg("-sampleFile",sampleFile,arguments);

if(doAsso==7){
fprintf(stderr,"Warning: Wald test is not properly tested and might give inflated test statistics!\n");
}
Expand All @@ -119,15 +122,20 @@ void abcAsso::getOptions(argStruct *arguments){

}

if(doAsso && yfile1==NULL && yfile2==NULL && yfile3==NULL){
if(doAsso && (yfile1!=NULL || yfile2!=NULL || yfile3!=NULL) && sampleFile!=NULL){
fprintf(stderr,"Error: you must provide either a phenotype file (-yBin or -yQuant) or sample file (-sampleFile) not both \n");
exit(0);
}

if(doAsso && yfile1==NULL && yfile2==NULL && yfile3==NULL && sampleFile==NULL){
fprintf(stderr,"Error: you must provide a phenotype file (-yBin or -yQuant) to perform association \n");
exit(0);
}
if(doAsso && arguments->inputtype==INPUT_BEAGLE&&doAsso==1){
}
if(doAsso && (arguments->inputtype==INPUT_BEAGLE&&doAsso==1 || arguments->inputtype==INPUT_BGEN&&doAsso==1)){
fprintf(stderr,"Error: Only doAsso=2 can be performed on posterior input\n");
exit(0);
}
if(doAsso && arguments->inputtype!=INPUT_BEAGLE&&(doAsso==2)&&doPost==0){
if(doAsso && arguments->inputtype!=INPUT_BEAGLE && arguments->inputtype!=INPUT_BGEN &&(doAsso==2)&&doPost==0){
fprintf(stderr,"Error: For doAsso=2 you must estimate the posterior probabilites for the genotypes (-doPost 1) \n");
exit(0);
}
Expand All @@ -154,6 +162,7 @@ abcAsso::abcAsso(const char *outfiles,argStruct *arguments,int inputtype){
yfile1=NULL;
yfile2=NULL;
yfile3=NULL;
sampleFile=NULL;
minHigh=10;
minCount=10;
dynCov=0;//not for users
Expand All @@ -172,8 +181,9 @@ abcAsso::abcAsso(const char *outfiles,argStruct *arguments,int inputtype){
if(!strcasecmp(arguments->argv[1],"-doAsso")){
printArg(stdout);
exit(0);
}else
} else{
return;
}
}

getOptions(arguments);
Expand All @@ -186,22 +196,43 @@ abcAsso::abcAsso(const char *outfiles,argStruct *arguments,int inputtype){
printArg(arguments->argumentFile);

//read phenotype - binary or count (poisson)
if(isBinary)
ymat = angsd::getMatrix(yfile1,1,100000);
else if(isCount)
ymat = angsd::getMatrix(yfile2,1,100000);
else
ymat = angsd::getMatrix(yfile3,0,100000);
if(sampleFile==NULL){

//read covariates
if(covfile!=NULL)
covmat = angsd::getMatrix(covfile,0,100000);
else{
covmat.x=0;
covmat.y=0;
covmat.matrix=NULL;
if(isBinary)
ymat = angsd::getMatrix(yfile1,1,100000);
else if(isCount)
ymat = angsd::getMatrix(yfile2,1,100000);
else
ymat = angsd::getMatrix(yfile3,0,100000);

//read covariates
if(covfile!=NULL)
covmat = angsd::getMatrix(covfile,0,100000);
else{
covmat.x=0;
covmat.y=0;
covmat.matrix=NULL;
}

} else{
dT = angsd::getSample(sampleFile,100000);

ymat.matrix=dT.matrix0;
ymat.x = dT.x0;
ymat.y = dT.y0;
isBinary = dT.isBinary;

if(dT.x1>0){;
covmat.matrix=dT.matrix1;
covmat.x = dT.x1;
covmat.y = dT.y1;
} else{
covmat.x=0;
covmat.y=0;
covmat.matrix=NULL;
}
}

if(covfile!=NULL&&(covmat.x!=ymat.x)){
fprintf(stderr,"The number of covariates (%d) does not match the number of phenotypes (%d)\n",covmat.x,ymat.x);
exit(0);
Expand All @@ -211,7 +242,7 @@ abcAsso::abcAsso(const char *outfiles,argStruct *arguments,int inputtype){
fprintf(stderr,"Error: Only doAsso=2 can be performed on quantitative traits\n");
exit(0);
}

// check cov and ymat
check_pars(covmat,ymat,isBinary);

Expand Down Expand Up @@ -329,9 +360,10 @@ abcAsso::~abcAsso(){


void abcAsso::clean(funkyPars *pars){

if(doAsso==0)
return;

assoStruct *assoc =(assoStruct*) pars->extras[index];

if(assoc->stat!=NULL)
Expand Down Expand Up @@ -369,8 +401,17 @@ void abcAsso::clean(funkyPars *pars){
for( int yi =0;yi<ymat.y;yi++)
delete[] assoc->keepInd[yi];
delete[] assoc->keepInd;


delete assoc;

//to delete if has not been deleted in other abc Classes
if(pars->post!=NULL){
for(int i=0;i<pars->numSites;i++)
delete [] pars->post[i];
delete [] pars->post;
pars->post =NULL;
}
}


Expand Down Expand Up @@ -409,13 +450,13 @@ void abcAsso::run(funkyPars *pars){

assoStruct *assoc = allocAssoStruct();
pars->extras[index] = assoc;

if(doAsso==2 or doAsso==4 or doAsso==5 or doAsso==6 or doAsso==7){
assoc->highWt=new int[pars->numSites];
assoc->highHe=new int[pars->numSites];
assoc->highHo=new int[pars->numSites];
}

if(doAsso==1){
frequencyAsso(pars,assoc);
} else if(doAsso==2){
Expand All @@ -440,8 +481,7 @@ void abcAsso::run(funkyPars *pars){
assoc->emIter=new int[pars->numSites];
emAssoWald(pars,assoc);
}



}

void abcAsso::check_pars(angsd::Matrix<double> &cov, angsd::Matrix<double> &phe, int isBinary){
Expand Down Expand Up @@ -1825,7 +1865,7 @@ double abcAsso::doEMasso(funkyPars *p,angsd::Matrix<double> *design,angsd::Matri

double* post = new double[keepInd*3];
double* y = new double[3*keepInd];

int count=0;

design->x=3*keepInd;
Expand Down Expand Up @@ -2097,14 +2137,14 @@ void abcAsso::emAsso(funkyPars *pars,assoStruct *assoc){
int **keepInd = new int*[ymat.y];
// we can also get coef of genotype
double **stat = new double*[ymat.y*2];

for(int yi=0;yi<ymat.y;yi++){
stat[yi] = new double[pars->numSites];
keepInd[yi]= new int[pars->numSites];
}

angsd::Matrix<double> designNull;
designNull.x=pars->nInd;
designNull.x=3*pars->nInd;
//covars + intercept
designNull.y=covmat.y+1;
designNull.matrix=new double*[3*pars->nInd];
Expand All @@ -2130,7 +2170,7 @@ void abcAsso::emAsso(funkyPars *pars,assoStruct *assoc){
for(int s=0;s<pars->numSites;s++){//loop overs sites
if(pars->keepSites[s]==0)
continue;

int *keepListAll = new int[pars->nInd];
for(int i=0 ; i<pars->nInd ;i++){
keepListAll[i]=1;
Expand All @@ -2152,6 +2192,7 @@ void abcAsso::emAsso(funkyPars *pars,assoStruct *assoc){
if(keepList[i]==1)
keepInd[yi][s]++;
}


double *y = new double[pars->nInd];
for(int i=0 ; i<pars->nInd ;i++)
Expand Down Expand Up @@ -2181,7 +2222,7 @@ void abcAsso::emAsso(funkyPars *pars,assoStruct *assoc){

delete [] start;

designNull.x=pars->nInd;
designNull.x=3*pars->nInd;
design.x=3*pars->nInd;
postAll.x=3*pars->nInd;

Expand Down Expand Up @@ -2520,7 +2561,7 @@ void abcAsso::hybridAsso(funkyPars *pars,assoStruct *assoc){
}

angsd::Matrix<double> designNull;
designNull.x=pars->nInd;
designNull.x=3*pars->nInd;
//covars + intercept
designNull.y=covmat.y+1;
designNull.matrix=new double*[3*pars->nInd];
Expand Down Expand Up @@ -2601,7 +2642,7 @@ void abcAsso::hybridAsso(funkyPars *pars,assoStruct *assoc){
delete [] start;
delete chisq1;

designNull.x=pars->nInd;
designNull.x=3*pars->nInd;
design.x=3*pars->nInd;
postAll.x=3*pars->nInd;

Expand Down Expand Up @@ -3459,10 +3500,10 @@ double abcAsso::poisScoreEnv(double *post,int numInds, double *y, double *ytilde
void abcAsso::printDoAsso(funkyPars *pars){
if(doPrint)
fprintf(stderr,"staring [%s]\t[%s]\n",__FILE__,__FUNCTION__);


freqStruct *freq = (freqStruct *) pars->extras[6];
assoStruct *assoc= (assoStruct *) pars->extras[index];

for(int yi=0;yi<ymat.y;yi++){
bufstr.l=0;
for(int s=0;s<pars->numSites;s++){
Expand All @@ -3483,6 +3524,9 @@ void abcAsso::printDoAsso(funkyPars *pars){
}
}
aio::bgzf_write(multiOutfile[yi],bufstr.s,bufstr.l);bufstr.l=0;

}


}

4 changes: 3 additions & 1 deletion abcAsso.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -58,7 +58,9 @@ class abcAsso:public abc{
char *yfile1;
char *yfile2;
char *yfile3;

char *sampleFile;

angsd::doubleTrouble<double> dT;
angsd::Matrix<double> ymat;
angsd::Matrix<double> covmat;
void check_pars(angsd::Matrix<double> &cov, angsd::Matrix<double> &phe, int isBinary);
Expand Down
19 changes: 9 additions & 10 deletions abcFreq.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -208,16 +208,16 @@ void abcFreq::getOptions(argStruct *arguments){

if(abs(doMaf)==0 &&doPost==0)
return;
if(inputtype!=INPUT_VCF_GP&&inputtype!=INPUT_BEAGLE&&inputtype!=0&&doMajorMinor==0&&inputtype!=INPUT_VCF_GL){
if(inputtype!=INPUT_VCF_GP&&inputtype!=INPUT_BEAGLE&&inputtype!=INPUT_BGEN&&inputtype!=0&&doMajorMinor==0&&inputtype!=INPUT_VCF_GL){
fprintf(stderr,"You must specify \'-doMajorMinor\' to infer major/minor \n");
exit(0);
}

if(inputtype==INPUT_BEAGLE&&abs(doMaf)!=4){
if((inputtype==INPUT_BEAGLE||inputtype==INPUT_BGEN)&&abs(doMaf)!=4){
fprintf(stderr,"Only \'-doMaf 4\' can not be performed on posterior input\n");
exit(0);
}
if(inputtype!=INPUT_BEAGLE&&inputtype!=INPUT_VCF_GP&&abs(doMaf)==4){
if(inputtype!=INPUT_BEAGLE&&inputtype!=INPUT_BGEN&&inputtype!=INPUT_VCF_GP&&abs(doMaf)==4){
fprintf(stderr,"\t \'-doMaf 4\' can only be performed on genotype probabilities provided by the user (-beagle).\n");
exit(0);
}
Expand All @@ -231,7 +231,7 @@ void abcFreq::getOptions(argStruct *arguments){
fprintf(stderr,"Error: -doMaf 4 cannot be used for a likelihood ratio test (doSNP) \n");
exit(0);
}
if(inputtype==INPUT_BEAGLE&&doPost){
if((inputtype==INPUT_BEAGLE||inputtype==INPUT_BGEN)&&doPost){
fprintf(stderr,"Error: Cannot estimate post (doPost) based on posterior probabilites\n");
exit(0);
}
Expand Down Expand Up @@ -470,16 +470,15 @@ void abcFreq::clean(funkyPars *pars) {
delete [] freq->lrt_EM_unknown;
delete [] freq->phat;
delete freq;



if(pars->post!=NULL){
for(int i=0;i<pars->numSites;i++)
delete [] pars->post[i];
for(int i=0;i<pars->numSites;i++){
delete [] pars->post[i];
}
delete [] pars->post;
pars->post =NULL;
}



}

freqStruct *allocFreqStruct(){
Expand Down