improve dominant/recessive model

ChangeLog

@@ -1,3 +1,7 @@
+2014-02-28  Xiaowei Zhan  <zhanxw@fantasia.csgstat.sph.umich.edu>
+
+	* Change: for dominant/recessive model, automatically generate covariance
+
 2014-02-27  Xiaowei Zhan  <zhanxw@fantasia.csgstat.sph.umich.edu>
 
 	* Change: for X chrom, also use male to get AF

README.md

@@ -5,7 +5,7 @@
 	- [Groupwise tests](#groupwise-tests)
 	- [Related individual tests](#related-individual-tests)
 	- [Meta-analysis tests](#meta-analysis-tests)
-		- [Dominant model and recessive model](#dominant-model-and-recessive-model)
+		- [Dominant models and recessive models](#dominant-models-and-recessive-models)
 - [Input files](#input-files)
 	- [Genotype file (VCF)](#genotype-file-vcf)
 	- [Phenotype file](#phenotype-file)
@@ -106,15 +106,13 @@ We support both unrelated individuals and related indivudlas (e.g. family data).
 
 The file `input.kinship` is calculated by `vcf2kinship` program, and usage to this program is described in [Related individual tests](#related-individual-tests).
 
-### Dominant model and recessive model
+### Dominant models and recessive models
 
 Dominant and recessive disease models are supported by appending "dominant" and/or "recessive" after "--meta" option. For example, use "--meta dominant,recessive" will 
-generate two files ".MetaDominant.assoc" and ".MetaRecessive.assoc". Details: In dominant models, genotypes 0/1/2 are coded as 0/1/1. In recessive models, genotypes 0/1/2 are 
+generate two sets of files. For dominant model, they are "prefix.MetaDominant.assoc" and "prefix.MetaDominantCov.assoc.gz"; for recessive model,
+they are "prefix.MetaRecessive.assoc" and "prefix.MetaRecessiveCov.assoc.gz". Internally, in dominant models, genotypes 0/1/2 are coded as 0/1/1; in recessive models, genotypes 0/1/2 are 
 coded as 0/0/1. Missing genotypes will be imputed to the mean.
 
-It is of intersts to meta-analyze the association reuslts of dominant/recessive model. That requires covariance matrix generated under dominant (or recessive) model.
-Use `--meta dominantCov` (or `--meta recessiveCov`) will generate `prefix.MetaDominantCov.assoc.gz` (or `prefix.MetaRecessiveCov.assoc.gz`).
-
 # Input files
 
 ## Genotype file (VCF)
@@ -252,11 +250,9 @@ beta distribution parameters for upweighting rare variants. Rvtests will output
 Type | Model(*)    |Traits(#) | Covariates | Related / unrelated | Description
 :--------------|:---------:|:------:|:----------:|:-------------------:|:-----------
 Score test          |  score      | Q  |     Y      |         R, U           | standard score tests
-Dominant coding score test       |  dominant   | Q  |     Y      |         R, U           | score tests, dominant disease model
-Recessive coding score test      |  recessive  | Q  |     Y      |         R, U           | score tests, recessive disease model
+Dominant model       |  dominant   | Q  |     Y      |         R, U           | score tests and covariance matrix under dominant disease model
+Recessive model      |  recessive  | Q  |     Y      |         R, U           | score tests and covariance matrix under recessive disease model
 Covariance          |  cov      | Q  |     Y      |         R, U           | covariance matrix
-Dominant coding covariance          |  dominantCov      | Q  |     Y      |         R, U           | covariance matrix using dominant coding
-Recessive coding covariance          |  recessiveCov      | Q  |     Y      |         R, U           | covariance matrix using recessive coding
 
 (*) Model columns list the regconized names in rvtests. For example, use `--meta score,cov` will generate score statistics and covariance matrix for meta-analysis.
 (#) In trait column, B and Q stand for (b)inary, (q)uantitiave trait.

base/Pedigree.cpp

@@ -79,11 +79,11 @@ int Pedigree::add(const std::string& family,
                   const std::string& mother)  {
 
   if (person != "0" && (person == father || person == mother) ) {
-    fprintf(stderr, "One and his father/mother are both [ %s ]\n", person.c_str());
+    fprintf(stderr, "ERROR: Individual cannot be his own father/mother [ %s ]\n", person.c_str());
     return -1;
   }
   if (father != "0" && mother != "0" && father == mother) {
-    fprintf(stderr, "Father and mother are both [ %s ]\n", father.c_str());
+    fprintf(stderr, "ERROR: Father and mother are the same individual [ %s ]\n", father.c_str());
     return -1;
   }
 
@@ -110,14 +110,14 @@ int Pedigree::add(const std::string& family,
   if (fatherId >= 0) {
     this->people[fatherId].addChild(pid);
     if (this->people[fatherId].setGender(MALE) < 0) {
-      fprintf(stderr, "Father but not male according to pedigree!\n");
+      fprintf(stderr, "ERROR: Father [ %s ] is not male according to pedigree!\n", father.c_str());
       errorOccured = true;
     };
   }
   if (motherId >= 0) {
     this->people[motherId].addChild(pid);
     if (this->people[motherId].setGender(FEMALE) < 0 ) {
-      fprintf(stderr, "Mother but not female according to pedigree!\n");
+      fprintf(stderr, "ERROR: Mother [ %s ] is not female according to pedigree!\n", mother.c_str());
       errorOccured = true;
     }
   }
@@ -133,6 +133,7 @@ int loadPedigree(const std::string& fn, zhanxw::Pedigree* ped) {
   std::vector<std::string> fd;
   LineReader lr(fn);
   int lineNo = 0;
+  bool errorOccured = false;
   while (lr.readLineBySep(&fd, " \t")) {
     ++lineNo;
     removeEmptyField(&fd);
@@ -145,10 +146,12 @@ int loadPedigree(const std::string& fn, zhanxw::Pedigree* ped) {
     }
     if (p.add(fd[0], fd[1], fd[2], fd[3]) < 0) {
       fprintf(stderr, "Encounter error when adding line %d.\n", lineNo);
+      errorOccured = true;
     }
     if (fd.size() >= 5 )
       p.addGender(fd[1], fd[4]);
   }
-
+
+  if (errorOccured) return -1;
   return 0;
 }

src/Main.cpp

@@ -36,7 +36,7 @@
 
 Logger* logger = NULL;
 
-#define VERSION "20140227"
+#define VERSION "20140228"
 
 void banner(FILE* fp) {
   const char* string =
@@ -603,7 +603,7 @@ int main(int argc, char** argv){
       ADD_STRING_PARAMETER(pl, modelBurden, "--burden", "cmc, zeggini, mb, exactCMC, rarecover, cmat, cmcWald")
       ADD_STRING_PARAMETER(pl, modelVT, "--vt", "cmc, zeggini, mb, skat")
       ADD_STRING_PARAMETER(pl, modelKernel, "--kernel", "SKAT, KBAC")
-      ADD_STRING_PARAMETER(pl, modelMeta, "--meta", "score, cov, dominant, recessive, dominantCov, recessiveCov")
+      ADD_STRING_PARAMETER(pl, modelMeta, "--meta", "score, cov, dominant, recessive")
 
       ADD_PARAMETER_GROUP(pl, "Family-based Models")
       ADD_STRING_PARAMETER(pl, kinship, "--kinship", "Specify a kinship file for autosomal analysis, use vcf2kinship to generate")
@@ -711,7 +711,7 @@ int main(int argc, char** argv){
   }
   if (FLAG_freqLower > 0) {
     vin.setSiteFreqMin(FLAG_freqLower);
-    logger->info("Set upper frequency limit to %f", FLAG_freqLower);
+    logger->info("Set lower frequency limit to %f", FLAG_freqLower);
   }
 
   // add filters. e.g. put in VCFInputFile is a good method
@@ -1098,20 +1098,18 @@ int main(int argc, char** argv){
         model.push_back( new MetaScoreTest() );
       } else if (modelName == "dominant") {
         model.push_back( new MetaDominantTest() );
+        parser.assign("windowSize", &windowSize, 1000000);
+        logger->info("Meta analysis uses window size %s to produce covariance statistics under dominant model", toStringWithComma(windowSize).c_str());
+        model.push_back( new MetaDominantCovTest(windowSize) );
       } else if (modelName == "recessive") {
         model.push_back( new MetaRecessiveTest() );
+        parser.assign("windowSize", &windowSize, 1000000);
+        logger->info("Meta analysis uses window size %s to produce covariance statistics under recessive model", toStringWithComma(windowSize).c_str());
+        model.push_back( new MetaRecessiveCovTest(windowSize) );
       } else if (modelName == "cov") {
         parser.assign("windowSize", &windowSize, 1000000);
         logger->info("Meta analysis uses window size %s to produce covariance statistics", toStringWithComma(windowSize).c_str());
         model.push_back( new MetaCovTest(windowSize) );
-      } else if (modelName == "dominantcov") {
-        parser.assign("windowSize", &windowSize, 1000000);
-        logger->info("Meta analysis uses window size %s to produce covariance statistics", toStringWithComma(windowSize).c_str());
-        model.push_back( new MetaDominantCovTest(windowSize) );
-      } else if (modelName == "recessivecov") {
-        parser.assign("windowSize", &windowSize, 1000000);
-        logger->info("Meta analysis uses window size %s to produce covariance statistics", toStringWithComma(windowSize).c_str());
-        model.push_back( new MetaRecessiveCovTest(windowSize) );
       }
 #if 0
       else if (modelName == "skew") {
@@ -1152,16 +1150,24 @@ int main(int argc, char** argv){
     std::string s = FLAG_outPrefix;
     s += ".";
     s += model[i]->getModelName();
-    s += ".assoc";
     if (model[i]->getModelName() == "MetaCov" ||
         model[i]->getModelName() == "MetaSkew" ||
-        model[i]->getModelName() == "MetaKurt" ||
-        model[i]->getModelName() == "MetaDominantCov" ||                
-        model[i]->getModelName() == "MetaRecessiveCov") {
-      s += ".gz";
+        model[i]->getModelName() == "MetaKurt") {
+      s += "assoc.gz";
       fOuts[i] = new FileWriter(s.c_str(), BGZIP);
       metaFileToIndex.push_back(s);
+    } else if (model[i]->getModelName() == "MetaDominant" ||                
+               model[i]->getModelName() == "MetaRecessive")  {
+      s += ".assoc";
+      fOuts[i] = new FileWriter(s.c_str());
+
+      ++i;
+      s.resize(s.size() - strlen(".assoc"));
+      s += "Cov.assoc.gz";
+      fOuts[i] = new FileWriter(s.c_str(), BGZIP);
+      metaFileToIndex.push_back(s);      
     } else {
+      s += ".assoc";
       fOuts[i] = new FileWriter(s.c_str());
     }
   }

vcfUtils/vcf2kinship.cpp

@@ -386,7 +386,7 @@ void usage(int argc, char** argv) {
 }
 
 #define PROGRAM "vcf2kinship"
-#define VERSION "20140227"
+#define VERSION "20140228"
 void welcome() {
 #ifdef NDEBUG
   fprintf(stdout, "Thank you for using %s (version %s)\n", PROGRAM, VERSION);
@@ -503,7 +503,7 @@ int main(int argc, char** argv){
   zhanxw::Pedigree ped;
   if (!FLAG_ped.empty()) {
     if (loadPedigree(FLAG_ped, &ped)) {
-      fprintf(stderr, "Failed to load pedigree file [ %s ]!", FLAG_ped.c_str());
+      fprintf(stderr, "Failed to load pedigree file [ %s ]!\n", FLAG_ped.c_str());
       exit(1);
     }
   }