BGI/SOAPsnp/prior.cc

   1 #include "soap_snp.h"
   2
   3 /**
   4  * DATE:
   5  * FUNCTION:     generate the prior matrix
   6  * PARAMETER:
   7  * RETURN :     1
   8  */
   9
  10 int Prob_matrix::prior_gen(Parameter * para) {
  11         char t_base, allele1, allele2;
  12         // Note, the above parameter should be changed to a more reasonable one
  13         for(t_base=0;t_base!=4;t_base++) {
  14                 for(allele1=0;allele1!=4;allele1++) {
  15                         for(allele2=allele1;allele2!=4;allele2++) {
  16                                 if(allele1 == t_base && allele2 == t_base) {
  17                                         // refHOM
  18                                         p_prior[t_base<<4|allele1<<2|allele2] = 1;
  19                                 }
  20                                 else if (allele1 == t_base || allele2 == t_base) {
  21                                         // refHET: 1 ref 1 alt
  22                                         p_prior[t_base<<4|allele1<<2|allele2] = para->het_novel_r;
  23                                 }
  24                                 else if (allele1 == allele2) {
  25                                         // altHOM
  26                                         p_prior[t_base<<4|allele1<<2|allele2] = para->althom_novel_r;
  27                                 }
  28                                 else {
  29                                         // altHET: 2 diff alt base
  30                                         p_prior[t_base<<4|allele1<<2|allele2] = para->het_novel_r * para->althom_novel_r;
  31                                 }
  32                                 if( para->transition_dominant && ((allele1^t_base) == 0x3 || (allele2^t_base) == 0x3)) {
  33                                         // transition
  34                                         p_prior[t_base<<4|allele1<<2|allele2] *= 4;
  35                                 }
  36                                 //std::cerr<<"ACTG"[t_base]<<"\t"<<"ACTG"[allele1]<<"ACTG"[allele2]<<"\t"<<p_prior[t_base<<4|allele1<<2|allele2]<<endl;
  37                         }
  38                 }
  39         }
  40         for(allele1=0;allele1!=4;allele1++) {
  41                 for(allele2=allele1;allele2!=4;allele2++) {
  42                         // Deal with N
  43                         p_prior[0x4<<4|allele1<<2|allele2] = (allele1==allele2? 1: (2*para->het_novel_r)) * 0.25 *0.25;
  44                         p_prior[0x5<<4|allele1<<2|allele2] = (allele1==allele2? 1: (2*para->het_novel_r)) * 0.25 *0.25;
  45                         p_prior[0x6<<4|allele1<<2|allele2] = (allele1==allele2? 1: (2*para->het_novel_r)) * 0.25 *0.25;
  46                         p_prior[0x7<<4|allele1<<2|allele2] = (allele1==allele2? 1: (2*para->het_novel_r)) * 0.25 *0.25;
  47                 }
  48         }
  49         return 1;
  50 }
  51
  52 int Call_win::snp_p_prior_gen(double * real_p_prior, Snp_info* snp, Parameter * para, char ref) {
  53         if (snp->is_indel()) {
  54                 return 0;
  55         }
  56         char base, allele1, allele2;
  57         int allele_count;
  58         allele_count = 0;
  59         for (base=0; base != 4; base ++) {
  60                 if(snp->get_freq(base)>0) {
  61                         // The base is found in dbSNP
  62                         allele_count += 1;
  63                 }
  64         }
  65         if(allele_count <= 1) {
  66                 // Should never occur
  67                 cerr<<"Previous Extract SNP error."<<endl;
  68                 exit(255);
  69                 return -1;
  70         }
  71         char t_base = (ref&0x3); // ????????????????
  72         for(allele1=0;allele1!=4;allele1++) {
  73                 for(allele2=allele1;allele2!=4;allele2++) {
  74                         if (! snp->is_hapmap()) {
  75                                 if(snp->get_freq(allele1)>0 && snp->get_freq(allele2)>0) {
  76                                         // Here the frequency is just a tag to indicate SNP alleles in non-HapMap sites
  77                                         if(allele1 == allele2 && allele1 == t_base) {
  78                                                 // refHOM
  79                                                 real_p_prior[allele1<<2|allele2] = 1;
  80                                         }
  81                                         else if (allele1 == t_base || allele2 == t_base) {
  82                                                 // refHET: 1 ref 1 alt
  83                                                 real_p_prior[allele1<<2|allele2] = snp->is_validated()?para->het_val_r:para->het_unval_r;
  84                                         }
  85                                         else if (allele1 == allele2) {
  86                                                 real_p_prior[allele1<<2|allele2] =  snp->is_validated()?para->althom_val_r:para->althom_unval_r;
  87                                         }
  88                                         else {
  89                                                 // altHET: 2 diff alt base
  90                                                 real_p_prior[allele1<<2|allele2] = snp->is_validated()?para->het_val_r:para->het_unval_r;
  91                                         }
  92                                 }
  93                         }
  94                         else {
  95                                 // Real HapMap Sites
  96                                 if(snp->get_freq(allele1)>0 && snp->get_freq(allele2)>0) {
  97                                         real_p_prior[allele1<<2|allele2] = (allele1==allele2?1:(2*para->het_val_r))*snp->get_freq(allele1)*snp->get_freq(allele2);
  98                                 }
  99                         }
 100                         //cerr<<"ACTG"[t_base]<<"\t"<<"ACTG"[allele1]<<"ACTG"[allele2]<<"\t"<<p_prior[t_base<<4|allele1<<2|allele2]<<endl;
 101                 }
 102         }
 103         return 1;
 104 }