← Back to branch summary

~ubuntu-branches/ubuntu/raring/plink/raring

« back to all changes in this revision

Viewing changes to plink.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Andreas Tille
  • Date: 2009-10-23 13:35:02 UTC
  • mfrom: (1.1.2 upstream)
  • Revision ID: james.westby@ubuntu.com-20091023133502-rmnkrfi314mevnvt
Tags: 1.07-1
* New upstream version
* Removed debian/patches/20_plink-1.06-gcc4.4.patch because it
  is applied upstream
* Standards-Version: 3.8.3 (no changes needed)
* Debhelper 7
* Build-Depends: zlib1g-dev

Show diffs side-by-side

added added

removed removed

Lines of Context:
plink.cpp
38
38
#include "sets.h"
39
39
#include "stats.h"
40
40
#include "idhelp.h"
 
41
#include "zed.h"
41
42
 
42
43
using namespace std;
43
44
 
50
51
 
51
52
int main(int argc, char* argv[]) 
52
53
{
53
 
  
 
54
 
 
55
 
54
56
  
55
57
  /////////////////////////
56
58
  // Setup, display title
60
62
  
61
63
  set_new_handler(NoMem);
62
64
 
63
 
  PVERSION = "1.06";        // 4 chars
 
65
  PVERSION = "1.07";        // 4 chars
64
66
  PREL = " ";               // space or p (full, or prelease) 
65
 
  PDATE    = "24/Apr/2009"; // 11 chars
 
67
  PDATE    = "10/Aug/2009"; // 11 chars
66
68
 
67
69
 
68
70
  //////////////////
193
195
 
194
196
 
195
197
  /////////////////////////
196
 
  // Input files
 
198
  // File compression utility
 
199
 
 
200
  if ( par::compress_file )
 
201
    {
 
202
      fileCompress();
 
203
      shutdown();
 
204
    }
 
205
 
 
206
  if ( par::uncompress_file )
 
207
    {
 
208
      fileUncompress();
 
209
      shutdown();
 
210
    }
 
211
 
 
212
 
 
213
 
 
214
 
 
215
  //////////////////////////////////////////////////
 
216
  // Main Input files
197
217
 
198
218
  // Simulate or read in data:
199
219
  // Binary or ASCII format; transposed/long/generic
200
220
 
201
221
  if (par::dummy) P.dummyLoader();
202
222
  else if (par::greport) P.displayGeneReport();
 
223
  else if (par::annot_file) P.annotateFile();
 
224
  else if (par::meta_analysis) P.metaAnalysis();
203
225
  else if (par::rare_test_score_range) P.displayRareRange();
204
 
  else if (par::simul) P.simulateSNPs();
 
226
  else if (par::simul) 
 
227
    {
 
228
      if ( par::simul_qt )
 
229
        P.simulateSNPs_QT();
 
230
      else
 
231
        P.simulateSNPs();      
 
232
    }
205
233
  else if (par::cnv_list) P.setUpForCNVList();
206
234
  else if (par::read_bitfile) P.readBinData();
207
235
  else if (par::lfile_input) P.readDataLongFormat();
211
239
    par::load_gvar=true; 
212
240
    P.readGenericVariantData();
213
241
  }
 
242
  else if ( par::dosage_assoc ) 
 
243
    {
 
244
      P.readFamFile(par::famfile);
 
245
      if ( par::dosage_hasMap )
 
246
        {
 
247
          checkFileExists( par::mapfile );
 
248
          vector<bool> include;
 
249
          vector<int> include_pos(0);
 
250
          int nvar = 0;
 
251
          P.readMapFile(par::mapfile,
 
252
                        include,
 
253
                        include_pos,
 
254
                        nvar);
 
255
        }
 
256
    }
214
257
 
215
258
  // Set number of individuals
216
259
  P.n = P.sample.size();
319
362
 
320
363
 
321
364
 
 
365
 
322
366
  //////////////////////////////////////////////////////////
323
367
  // Output a specific set of SNPs (--extract or --exclude)
324
368
  
465
509
    }
466
510
  else if ( par::sol_family ) 
467
511
    {
468
 
      P.printLOG("Setting to permute within family ID only\n");
 
512
      P.printLOG("Setting clusters based on family IDs\n");
469
513
      vector<string> famlist;
470
514
      P.kname.resize(0);
471
515
      for (int i=0; i<P.n; i++)
516
560
  if ( par::zero_cluster )
517
561
    P.zeroOnCluster();
518
562
  
519
 
  
 
563
    
 
564
  /////////////////////////////////
 
565
  // Fix reference allele? 
 
566
 
 
567
  if ( par::set_reference_allele ) 
 
568
    P.setReferenceAllele();
 
569
 
 
570
 
520
571
 
521
572
  //////////////////////////////////
522
573
  // Determine formats for printing
525
576
  
526
577
 
527
578
 
 
579
 
 
580
  //////////////////////////////////////////////////
 
581
  //                                              //
 
582
  // Process a dosage file                        //
 
583
  //                                              //
 
584
  //////////////////////////////////////////////////
 
585
 
 
586
  if ( par::dosage_assoc )
 
587
    {
 
588
      // Normal behavior is to load data, and perform 
 
589
      // analysis; if the hard-call option is specified, 
 
590
      // then this will generate a dataset, that we can
 
591
      // subsequent filter and save, etc, as usual, i.e.
 
592
      // in that case, do not halt
 
593
 
 
594
      P.processDosageFile();
 
595
      
 
596
      if ( ! par::dosage_hard_call )
 
597
        shutdown();
 
598
 
 
599
    }
 
600
 
 
601
 
 
602
 
 
603
 
528
604
  //////////////////////////////////////////////////
529
605
  //                                              //
530
606
  // Handle CNV segments separately               //
566
642
      P.processGVAR();
567
643
      shutdown();
568
644
    }
 
645
 
 
646
 
 
647
 
 
648
  //////////////////////////////////////////////////
 
649
  //                                              //
 
650
  // Misc. .genome grouper utility                //
 
651
  //                                              //
 
652
  //////////////////////////////////////////////////
569
653
  
 
654
  if ( par::genome_groups )    
 
655
    {
 
656
      P.groupGenome();
 
657
      shutdown();
 
658
    }
570
659
 
571
660
 
572
661
  //////////////////////////////////
576
665
//     par::missing_phenotype = "0";
577
666
  
578
667
 
 
668
 
579
669
  //////////////////////////////////////////////////
580
670
  //                                              //
581
671
  // Basic MAF, genotyping filters & HWE/ME       //
615
705
  // Re-report basic case/control counts, etc
616
706
 
617
707
  summaryBasics(P);
618
 
 
619
708
  
 
709
 
 
710
 
620
711
  //////////////////////////////////////////////////
621
712
  // Any null allele codes (monomorhpics)?
622
713
  
626
717
        P.locus[l]->allele1 = "0";
627
718
    }
628
719
 
629
 
  
630
 
  /////////////////////////////////
631
 
  // Fix reference allele? 
632
 
 
633
 
  if (par::set_reference_allele ) 
634
 
    P.setReferenceAllele();
635
 
 
636
720
 
637
721
  /////////////////////////////////////////
638
722
  // SET statistics?
646
730
  P.pS = & S;
647
731
 
648
732
  // Remove any SNPs not in a set
649
 
  // unless outputing the SET TABLE
 
733
  // unless using particular commands
 
734
  // (set-by-all epistasis, set-table)
650
735
  
651
736
  if ( par::read_set || par::make_set )
652
737
    {
653
 
      if ( ! par::set_table ) 
 
738
      if ( par::drop_sets ) 
654
739
        P.pS->dropNotSet(P);
655
740
    }
656
741
  
674
759
  if (par::MENDEL_test || 
675
760
      par::MENDEL_report || 
676
761
      par::TDT_test ||
677
 
      par::QTDT_test ||
 
762
      par::QTDT_test ||      
678
763
      par::make_founders &&
679
764
      !par::built_families)
680
765
    {
931
1016
      shutdown();
932
1017
    }
933
1018
 
 
1019
  if (par::recode_long) 
 
1020
    {
 
1021
      P.display_recoded_LONG();
 
1022
      if (par::clist)
 
1023
        P.write_covariates();
 
1024
      shutdown();
 
1025
    }
 
1026
 
 
1027
  if (par::recode_mutlist) 
 
1028
    {
 
1029
      P.display_recoded_MUTLIST();
 
1030
      if (par::clist)
 
1031
        P.write_covariates();
 
1032
      shutdown();
 
1033
    }
 
1034
  
934
1035
 
935
1036
  if (par::list_by_allele)
936
1037
    {