bin/loading_scripts/load_genotypes.pl

   1
   2
   3 use strict;
   4
   5 use Getopt::Std;
   6 use CXGN::DB::InsertDBH;
   7 use CXGN::Marker;
   8 use CXGN::Phenome::Genotype;
   9 use CXGN::Phenome::GenotypeRegion;
  10 use CXGN::Phenome::GenotypeExperiment;
  11 use CXGN::Phenome::Individual;
  12 use CXGN::Cview::Map::Tools;
  13 use CXGN::Marker::Tools;
  14
  15 our ($opt_D, $opt_H, $opt_s, $opt_m, $opt_b, $opt_p);
  16
  17 getopts('D:H:s:m:b:p:');
  18
  19 if (!$opt_p) {
  20     print STDERR " -p population_id is required. Stop.\n";
  21     exit();
  22 }
  23
  24 if (!$opt_s) {
  25     $opt_s = "206";
  26     print STDERR "no -s specified. Using default of 206.\n";
  27 }
  28 my $sp_person_id=$opt_s;
  29
  30 if (!$opt_m) {
  31    $opt_m = 17;
  32    print STDERR "Using a default of $opt_m for the reference map id.\n";
  33 }
  34 my $reference_map_id=$opt_m;
  35
  36 if (!$opt_b) {
  37     $opt_b = 100;
  38     print STDERR "Using a default background accession id of $opt_b.\n";
  39 }
  40 my $background_accession_id=$opt_b;
  41
  42 my $file = shift;
  43
  44 my $dbh = CXGN::DB::InsertDBH->new({ dbname=>$opt_D,
  45                                      dbhost=>$opt_H,
  46                                      dbuser=>"postgres"
  47                                      });
  48
  49 if (!$file) {
  50     print <<TEXT;
  51
  52 Usage: load_genotypes.pl -H host -D dbname [-s sp_person_id] file\n
  53
  54 TEXT
  55
  56     exit();
  57 }
  58
  59 my $map_version_id = CXGN::Cview::Map::Tools::find_current_version($dbh, $reference_map_id);
  60
  61 open(F, "<$file") || die "Can't open file '$file'\n";
  62
  63 #my $first_line = (<F>); # throw away, contains chr info
  64 my $header = (<F>);
  65 chomp($header);
  66
  67 my (undef, @markers) = split /\t/, $header;
  68
  69 print STDERR  "Markers @markers\n";
  70
  71
  72 eval {
  73     my $experiment = CXGN::Phenome::GenotypeExperiment->new($dbh);
  74
  75     $experiment->set_background_accession_id($background_accession_id);
  76     $experiment->set_experiment_name("QTL experiment");
  77     $experiment->set_reference_map_id($reference_map_id);
  78     $experiment->set_sp_person_id($sp_person_id);
  79     $experiment->set_preferred(1);
  80     my $experiment_id = $experiment->store();
  81
  82
  83     while (<F>) {
  84         chomp;
  85         my ($individual_name, @scores) = split /\t/;
  86         print STDERR "Processing individual $individual_name...\n";
  87         my (@individuals) = CXGN::Phenome::Individual->new_with_name($dbh, $individual_name);
  88
  89         if (@individuals>1) {
  90             print STDERR "Found several individuals with name $individual_name. Skipping.\n";
  91             next();
  92         }
  93         my $individual = $individuals[0];
  94
  95         if ($individual && $individual->get_population_id() != $opt_p) {
  96             print STDERR "Individual $individual_name is not from population $opt_p\n";
  97             exit();
  98         }
  99
 100
 101         if (!$individual) {
 102             print STDERR "Individual $individual_name not found. \n";
 103
 104             print STDERR "Can't continue without individual $individual_name\n";
 105             print STDERR "Insert individual? Give population id or type N:";
 106             my $response = <STDIN>;
 107             chomp($response);
 108             if ($response=~/n/i) { next(); }
 109
 110             else {
 111
 112                 print STDERR "Creating new individual...\n";
 113                 $individual=CXGN::Phenome::Individual->new($dbh);
 114                 $individual->set_name($individual_name);
 115                 $individual->set_population_id($opt_p); #!!!!!!!!!!!
 116                 my $individual_id = $individual->store();
 117
 118                 print STDERR "Inserted Individual $individual_name with id $individual_id\n";
 119             }
 120
 121         }
 122
 123
 124
 125         my $genotype = CXGN::Phenome::Genotype->new($dbh);
 126
 127         print STDERR "experiment_id is $experiment_id\n";
 128         $genotype->set_genotype_experiment_id($experiment_id);
 129         $genotype->set_individual_id($individual->get_individual_id());
 130
 131
 132         my $genotype_id = $genotype->store();
 133
 134         for(my $i=0; $i<@scores; $i++) {
 135             my $genotype_region = CXGN::Phenome::GenotypeRegion->new($dbh);
 136
 137             my $clean_name = $markers[$i];
 138             $clean_name=~s/\*//g;
 139             $clean_name=~s/^ *(.*) *$/$1/; #clean leading/trailing spaces
 140             $clean_name = CXGN::Marker::Tools::clean_marker_name($clean_name);
 141             my $marker = CXGN::Marker->new_with_name($dbh, $clean_name);
 142
 143             if (!$scores[$i] || ($scores[$i] =~/\-/)) {
 144                 print STDERR "No zygocity information ($scores[$i]). Skipping...\n";
 145                 next();
 146             }
 147
 148             if (!$marker) {
 149                 print STDERR "Marker $clean_name not found. Skipping...\n";
 150                 next();
 151             }
 152
 153             my $experiments = $marker->experiments();
 154             my $lg_id = 0;
 155             foreach my $e (@$experiments) {
 156                 if (defined($e->{location})) {
 157
 158                     if ($e->{location}->map_version_id() eq $map_version_id) {
 159                         #print STDERR "Defined location for marker $markers[$i] on map $map_version_id...\n";
 160                         $lg_id = $e->{location}->lg_id();
 161                     }
 162                 }
 163             }
 164
 165             if (!$lg_id) {
 166                 print STDERR "Marker $markers[$i] has no mapped locations on map version $map_version_id. Inserting undef for lg_id \n";
 167
 168                 $lg_id = undef;
 169             }
 170             $genotype_region->set_phenome_genotype_id($genotype_id);
 171             $genotype_region->set_marker_id_nn($marker->marker_id());
 172             $genotype_region->set_marker_id_ns($marker->marker_id());
 173             $genotype_region->set_marker_id_sn($marker->marker_id());
 174             $genotype_region->set_marker_id_ss($marker->marker_id());
 175             $genotype_region->set_lg_id($lg_id);
 176             print STDERR "Setting score: $scores[$i]\n";
 177             $genotype_region->set_mapmaker_zygocity_code($scores[$i]);
 178
 179             $genotype_region->set_type("map");
 180
 181             $genotype_region->store();
 182
 183         }
 184     }
 185
 186 };
 187
 188 if ($@) {
 189     $dbh->rollback();
 190     print STDERR "An error occurred: $@. ROLLED BACK CHANGES.\n";
 191 }
 192 else {
 193     print STDERR "All fine. Committing...\n";
 194     $dbh->commit();
 195 }
 196 print STDERR "Done.\n";