Main.cpp

   1 ////////////////////////////////////////////////////////////////////// \r
   2 // Main.cpp \r
   3 //////////////////////////////////////////////////////////////////////////////\r
   4 //              COPYRIGHT NOTICE FOR GENOME CODE\r
   5 //\r
   6 // Copyright (C) 2006 - 2009, Liming Liang and Goncalo Abecasis,\r
   7 // All rights reserved.\r
   8 //\r
   9 //   Redistribution and use in source and binary forms, with or without\r
  10 //   modification, are permitted provided that the following conditions\r
  11 //   are met:\r
  12 //\r
  13 //     1. Redistributions of source code must retain the above copyright\r
  14 //        notice, this list of conditions and the following disclaimer.\r
  15 //\r
  16 //     2. Redistributions in binary form must reproduce the above copyright\r
  17 //        notice, this list of conditions and the following disclaimer in the\r
  18 //        documentation and/or other materials provided with the distribution.\r
  19 //\r
  20 //     3. The names of its contributors may not be used to endorse or promote\r
  21 //        products derived from this software without specific prior written\r
  22 //        permission.\r
  23 //\r
  24 //   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS\r
  25 //   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT\r
  26 //   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR\r
  27 //   A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR\r
  28 //   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,\r
  29 //   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,\r
  30 //   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR\r
  31 //   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF\r
  32 //   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING\r
  33 //   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS\r
  34 //   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.\r
  35 //\r
  36 ///////////////////////////////////////////////////////////////////////////////\r
  37 \r
  38 \r
  39 #include <iostream>\r
  40 #include <math.h>\r
  41 #include <vector> \r
  42 #include <time.h>\r
  43 #include "genome.h"\r
  44 #include "Random.h"\r
  45 #include <stdio.h> \r
  46 #include <iomanip>\r
  47     \r
  48 using namespace std;\r
  49 \r
  50 \r
  51 int main(int argc, char ** argv)\r
  52    {\r
  53         string popSize("10000"), rec("0.0001");\r
  54         \r
  55         int numPieces=100, pieceLen=10000, numIndepRegion=1, SNP=-1;\r
  56         int nSubPOP=2; \r
  57         vector<int> nSubSample(nSubPOP);\r
  58         for(int i=0;i<nSubPOP;i++)nSubSample[i]=10;\r
  59         \r
  60         double mut=1e-8, mig=2.5e-4;\r
  61         double maf=0.0,prop=1.0;\r
  62         \r
  63         vector< vector<bool> > data;\r
  64 \r
  65         long seed=-1;\r
  66         \r
  67         int drawtree=0;\r
  68 \r
  69         if(argc==1){\r
  70                 cout << endl  << "GENOME-0.2: Whole Genome Coalescent Simulator. (2009.2) Liming Liang, Goncalo Abecasis"  << endl << endl;\r
  71                 cout << "     Parameters and default values:" << endl;\r
  72                 cout << "     -pop     " << "number of subpopulations and size of subsamples [ " <<nSubPOP<<" "; \r
  73                                                                         for(int i=0;i<nSubPOP;i++)cout<<nSubSample[i]<<" ";cout<<"]"<< endl;\r
  74                 cout << "     -N       " << "effective size of each subpopulation or the filename for population profile [" <<popSize.c_str()<<"]"<< endl;\r
  75                 cout << "     -c       " << "number of independent regions (chromosomes) to simulate [" <<numIndepRegion<<"]"<< endl;\r
  76                 cout << "     -pieces  " << "number of fragments per independent region [" <<numPieces<<"]"<< endl;\r
  77                 cout << "     -len     " << "length in base of each fragment [" <<pieceLen<<"]"<< endl;\r
  78                 cout << "     -s       " << "fixed number of SNPs per independent region (chromosome), -1 = number of SNPs follows Poisson distribution [" <<SNP<<"]"<< endl;\r
  79                 cout << "     -rec     " << "recombination rate bewteen consecutive fragments"<<" or the filename for recombination rate distribution ["<<rec.c_str()<<"]"<< endl;\r
  80                 cout << "     -mut     " << "mutation rate per generation per base pair [" <<mut<<"]"<< endl;\r
  81                 cout << "     -mig     " << "migration rate per generation per individual [" <<mig<<"]"<< endl;\r
  82                 cout << "     -seed    " << "random seed, -1 = use time as the seed [" <<seed<<"]"<< endl;\r
  83                 cout << "     -tree    " << "1=draw the genealogy trees, 0=do not output [" <<drawtree<<"]"<< endl;\r
  84                 cout << "     -maf     " << "Output SNPs with minor allele frequency greater than [" <<maf<<"]"<< endl;\r
  85                 cout << "     -prop    " << "To keep this proportion of SNPs with MAF < the value of -maf parameter [" <<prop<<"]"<< endl;\r
  86                 exit(0);\r
  87         }\r
  88  \r
  89 \r
  90         int k=0;\r
  91     for(int i=1;i<argc;i++){\r
  92         if(strcmp(argv[i],"-pop")==0){ \r
  93                 nSubPOP=atoi(argv[++i]); k++;\r
  94                 nSubSample.clear();\r
  95                 nSubSample.resize(nSubPOP);\r
  96                 for(int iter=0;iter<nSubPOP;iter++)nSubSample[iter]=atoi(argv[++i]);\r
  97             }\r
  98             else if(strcmp(argv[i],"-N")==0){ popSize.assign(argv[++i]); k++;}\r
  99             else if(strcmp(argv[i],"-c")==0){ numIndepRegion=atoi(argv[++i]); k++;}\r
 100             else if(strcmp(argv[i],"-pieces")==0){ numPieces=atoi(argv[++i]); k++;}\r
 101             else if(strcmp(argv[i],"-len")==0){ pieceLen=atoi(argv[++i]); k++;}\r
 102             else if(strcmp(argv[i],"-s")==0){ SNP=atoi(argv[++i]); k++;}\r
 103             else if(strcmp(argv[i],"-rec")==0){ rec.assign(argv[++i]); k++;}\r
 104             else if(strcmp(argv[i],"-mut")==0){ mut=atof(argv[++i]); k++;}\r
 105             else if(strcmp(argv[i],"-mig")==0){ mig=atof(argv[++i]); k++;}\r
 106             else if(strcmp(argv[i],"-seed")==0){ seed=atol(argv[++i]); k++;}\r
 107             else if(strcmp(argv[i],"-tree")==0){ drawtree=atoi(argv[++i]); k++;}\r
 108             else if(strcmp(argv[i],"-maf")==0){ maf=atof(argv[++i]); k++;}\r
 109             else if(strcmp(argv[i],"-prop")==0){ prop=atof(argv[++i]); k++;}\r
 110     }\r
 111     \r
 112     if(k!=(argc-1)/2 && k!=(argc-1-nSubPOP)/2){\r
 113         cout << "Parameters do not match! " << argc<< " " << k << endl; \r
 114         exit(1);\r
 115     } \r
 116 \r
 117     cout << endl \r
 118             << "GENOME-0.2: Whole Genome Coalescent Simulator. (2009.2) Liming Liang, Goncalo Abecasis"  << endl << endl;\r
 119         cout << "     Parameters and effective values:" << endl;\r
 120         cout << "     -pop     " << "number of subpopulations and size of subsamples [ " <<nSubPOP<<" "; \r
 121                                                                         for(int i=0;i<nSubPOP;i++)cout<<nSubSample[i]<<" ";cout<<"]"<< endl;\r
 122                 cout << "     -N       " << "effective size of each subpopulation or the filename for population profile [" <<popSize.c_str()<<"]"<< endl;\r
 123                 cout << "     -c       " << "number of independent regions (chromosomes) to simulate [" <<numIndepRegion<<"]"<< endl;\r
 124                 cout << "     -pieces  " << "number of fragments per independent region [" <<numPieces<<"]"<< endl;\r
 125                 cout << "     -len     " << "length in base of each fragment [" <<pieceLen<<"]"<< endl;\r
 126                 cout << "     -s       " << "fixed number of SNPs per independent region (chromosome), -1 = number of SNPs follows Poisson distribution [" <<SNP<<"]"<< endl;\r
 127                 cout << "     -rec     " << "recombination rate bewteen consecutive fragments or the filename for recombination rate distribution ["<<rec.c_str()<<"]"<< endl;\r
 128                 cout << "     -mut     " << "mutation rate per generation per base pair [" <<mut<<"]"<< endl;\r
 129                 cout << "     -mig     " << "migration rate per generation per individual [" <<mig<<"]"<< endl;\r
 130                 cout << "     -seed    " << "random seed, -1 = use time as the seed [" <<seed<<"]"<< endl;\r
 131                 cout << "     -tree    " << "1=draw the genealogy trees, 0=do not output [" <<drawtree<<"]"<< endl;\r
 132                 cout << "     -maf     " << "Output SNPs with minor allele frequency greater than [" <<maf<<"]"<< endl;\r
 133                 cout << "     -prop    " << "To keep this proportion of SNPs with MAF < the value of -maf parameter [" <<prop<<"]"<< endl<<endl;\r
 134 \r
 135         if(atoi(popSize.c_str())>0 && atof(rec.c_str())>0){     \r
 136                 printf("     Scaled recombination rate = %.4e\n", 2 * atoi(popSize.c_str()) * (numPieces - 1) * atof(rec.c_str()));\r
 137                 printf("     Scaled migration rate = %.4e\n", 2 * atoi(popSize.c_str()) * mig);\r
 138                 printf("     Scaled mutation rate = %.4e\n\n", 2 * atoi(popSize.c_str()) * mut * numPieces * pieceLen);\r
 139         }\r
 140         \r
 141         genome(popSize, nSubPOP, nSubSample, numPieces, pieceLen, numIndepRegion, SNP, rec, mut, mig, data, seed,drawtree);\r
 142         \r
 143 // output format for haploview\r
 144 \r
 145 //      for(int i=0;i<data.size();i++){\r
 146 //              printf("FAMILY1 MEMBER%d ",i/2+1);\r
 147 //              for(int j=0;j<data[i].size();j++){\r
 148 //                      cout<<data[i][j]+1<<" ";\r
 149 //              }\r
 150 //              cout<<endl;\r
 151 //      }\r
 152 \r
 153 \r
 154 // 0 = ancestral state, 1 = mutated state\r
 155 \r
 156         cout<<"Total number of SNPs = "<<data[0].size()<<endl;\r
 157         cout<<"Samples:"<<endl;\r
 158         \r
 159         if(maf>0.0 && prop < 1.0){ // filter out the SNPs with MAF < maf \r
 160                 \r
 161                 cout<<"Rare SNPs filtered according to parameters."<<endl;\r
 162                 \r
 163                 int n=0;\r
 164                 for(int k=0;k<nSubPOP;k++)n+=nSubSample[k];\r
 165                 \r
 166                 int m=(int)(data[0].size());\r
 167                 vector<bool> keep(m,true);\r
 168                         \r
 169                 float p;\r
 170                 for(int i=0;i<m;i++){\r
 171                         p=0.0;\r
 172                         for(int j=0;j<n;j++){\r
 173                                 p += data[j][i];        \r
 174                         }\r
 175                         p /= n;\r
 176                         \r
 177                         if( (p<maf || (1-p)<maf) && globalRandom.Next()>prop)keep[i]=0;\r
 178                 }\r
 179                 \r
 180 \r
 181                 int index=0;\r
 182                 for(int k=0;k<nSubPOP;k++){\r
 183                         for(int i=0;i<nSubSample[k];i++){\r
 184                                 cout<<"POP"<<k+1<<": ";\r
 185                                 for(unsigned int j=0;j<data[index].size();j++){\r
 186                                         if(keep[j])cout<<data[index][j];\r
 187                                 }\r
 188                                 cout<<endl;\r
 189                                 index++;\r
 190                         }       \r
 191                 }\r
 192         \r
 193         }\r
 194         else{\r
 195                 int index=0;\r
 196                 for(int k=0;k<nSubPOP;k++){\r
 197                         for(int i=0;i<nSubSample[k];i++){\r
 198                                 cout<<"POP"<<k+1<<": ";\r
 199                                 for(unsigned int j=0;j<data[index].size();j++){\r
 200                                         cout<<data[index][j];\r
 201                                 }\r
 202                                 cout<<endl;\r
 203                                 index++;\r
 204                         }       \r
 205                 }\r
 206         }\r
 207 \r
 208   \r
 209 }\r
 210 \r
 211 \r