#include "iqlbestimator.h" IQLBEstimator::IQLBEstimator(vector & fam_list, vector & geno_trio_list, vector & missing_trio_list, map & hap_list, int size, int oneparent) :fam_list(fam_list), geno_trio_list(geno_trio_list), missing_trio_list(missing_trio_list), hap_freq_table(hap_list), loci_num(size), oneparent(oneparent){ sample_size=0; for(int i=0;iN_typed; init_freq_table(); compute_count(); if(oneparent==1) compute_hst(); geno_missing_trio_patten(); pair_missing_pattern(); flag_pd = 1; df = 0; pvalue = 0; test_perform = 0; } void IQLBEstimator::init_freq_table(){ int tot; if(hap_freq_table.size()==0){ tot = 1<::iterator iter=hap_freq_table.begin(); while(iter!=hap_freq_table.end()){ iter->second = 1.0/double(tot); iter++; } } geno_freq_table.clear(); map::iterator iter1=hap_freq_table.begin(); map::iterator iter2; while(iter1!=hap_freq_table.end()){ iter2 = iter1; while(iter2!=hap_freq_table.end()){ int h1 = iter1->first; int h2 = iter2->first; if(h1 == h2) geno_freq_table[pair(h1,h2)] = iter1->second*iter2->second; else geno_freq_table[pair(h1,h2)] = 2*iter1->second*iter2->second; iter2++; } iter1++; } } void IQLBEstimator::reset_freq_table(map & hap_freq){ hap_freq_table.clear(); hap_freq_table = hap_freq; geno_freq_table.clear(); map::iterator iter1=hap_freq_table.begin(); map::iterator iter2; while(iter1!=hap_freq_table.end()){ iter2 = iter1; while(iter2!=hap_freq_table.end()){ int h1 = iter1->first; int h2 = iter2->first; if(h1 == h2) geno_freq_table[pair(h1,h2)] = iter1->second*iter2->second; else geno_freq_table[pair(h1,h2)] = 2*iter1->second*iter2->second; iter2++; } iter1++; } } void IQLBEstimator::show_hap_freq(){ fprintf(outfp, "SNP haplotype frequency estimates in full sample\n"); fprintf(outfp, "Haplotype\t Frequency (choice b)\n"); map::iterator iter = hap_freq_table.begin(); while(iter!=hap_freq_table.end()){ if(iter->second<0.00005){ iter++; continue; } translate_hap(outfp,iter->first,loci_num); fprintf(outfp, " \t %.4f\n", iter->second); iter++; } fprintf(outfp, "\n"); } void IQLBEstimator::show_hap_freq(vector &a0v, vector &a1v){ fprintf(outfp, "SNP haplotype frequency estimates in full sample\n"); fprintf(outfp, "Haplotype\t Frequency (choice b)\n"); map::iterator iter = hap_freq_table.begin(); while(iter!=hap_freq_table.end()){ if(iter->second<0.00005){ iter++; continue; } translate_hap(outfp,iter->first,loci_num,a0v,a1v); if(loci_num==2) fprintf(outfp, "\t"); fprintf(outfp, " \t %.4f\n", iter->second); iter++; } fprintf(outfp, "\n"); } void IQLBEstimator::show_geno_freq(){ map, double>::iterator giter = geno_freq_table.begin(); while(giter!=geno_freq_table.end()){ if(giter->second<0.00005){ giter++; continue; } printf("(%2d,%2d) ", (giter->first).first,(giter->first).second); printf(" %.4f\n", giter->second); giter++; } printf("\n"); } void IQLBEstimator::compute_count(){ int row = hap_freq_table.size(); int col = geno_freq_table.size(); count = new int*[row-1]; for(int i=0;i::iterator iter=hap_freq_table.begin(); int loc = 0; while(iter!=hap_freq_table.end()){ hap_map[iter->first] = loc; loc++; iter++; } loc = 0; map, double>::iterator giter=geno_freq_table.begin(); while(giter!=geno_freq_table.end()){ int h1=(giter->first).first; int h2=(giter->first).second; int loc1 = hap_map[h1]; int loc2 = hap_map[h2]; if(loc1!=row-1) count[loc1][loc]++; if(loc2!=row-1) count[loc2][loc]++; loc++; giter++; } } void IQLBEstimator::compute_hst(){ int row = hap_freq_table.size(); int col = geno_freq_table.size(); hst = new double**[row-1]; for(int i=0;i, double>::iterator giter1; map, double>::iterator giter2; loc1 = 0; giter1=geno_freq_table.begin(); while(giter1!=geno_freq_table.end()){ int i=(giter1->first).first; int j=(giter1->first).second; loc2 = 0; giter2=geno_freq_table.begin(); while(giter2!=geno_freq_table.end()){ int k=(giter2->first).first; int l=(giter2->first).second; if(k==l){ // parent is homozygote if(i==k){ loc = hap_map[j]; if(loc!=row-1) hst[loc][loc1][loc2]+=1; else{ for(int s=0;styped_o; typed_par1 = geno_trio_list[i]->typed_par1; typed_par2 = geno_trio_list[i]->typed_par2; if(typed_par1==1 || typed_par2==1){ int *a_o=(geno_trio_list[i]->offspring)->get_gcode(); int *a_f=(geno_trio_list[i]->par1)->get_gcode(); int *a_m=(geno_trio_list[i]->par2)->get_gcode(); for(int j=0;jtyped_o || typed_par1!=missing_trio_list[j]->typed_par1 || typed_par2!=missing_trio_list[j]->typed_par2){ flag = 0; continue; } int *b_o=(missing_trio_list[j]->M_offspring)->get_gcode(); int *b_f=(missing_trio_list[j]->M_par1)->get_gcode(); int *b_m=(missing_trio_list[j]->M_par2)->get_gcode(); flag = compare_equ(a_o,b_o,loci_num); if(flag==0) continue; if((compare_equ(a_f,b_f,loci_num)==1 && compare_equ(a_m,b_m,loci_num)==1) || (compare_equ(a_f,b_m,loci_num)==1 && compare_equ(a_m,b_f,loci_num)==1)){ geno_missing_trio_map[i] = j; break; } } } else if(typed_par1==0 && typed_par2==0){ int *a_o=(geno_trio_list[i]->offspring)->get_gcode(); for(int j=0;jtyped_o || typed_par1!=missing_trio_list[j]->typed_par1 || typed_par2!=missing_trio_list[j]->typed_par2){ flag = 0; continue; } int *b_o=(missing_trio_list[j]->M_offspring)->get_gcode(); flag = compare_equ(a_o,b_o,loci_num); if(flag==0) continue; geno_missing_trio_map[i] = j; break; } } } /* if(geno_missing_trio_map.size()!=geno_trio_list.size()) { printf("geno_missing_trio_map assignment is wrong. Please check.\n"); exit(1); } */ } void IQLBEstimator::show_test(vector *cluster, vector &a0v, vector &a1v){ map trans; map::iterator iter = hap_freq_table.begin(); int count = 0; while(iter!=hap_freq_table.end()){ trans[count] = iter->first; iter++; count++; } if(df<(1< dis; for(int i=0;i::iterator diter = dis.begin(); while(diter!=dis.end()){ translate_hap(outfp,diter->first,loci_num,a0v,a1v); fprintf(outfp, "\t"); if(loci_num==2) fprintf(outfp, "\t"); fprintf(outfp, " \t %.6f\t %g\n", Btest[diter->second],Bpvalue[diter->second]); diter++; } fprintf(outfp, "\n"); trans.clear(); dis.clear(); } void IQLBEstimator::update_hap_freq(){ map::iterator iter = hap_freq_table.begin(); while(iter!=hap_freq_table.end()){ iter->second = 0; iter++; } map, double>::iterator giter = geno_freq_table.begin(); while(giter!= geno_freq_table.end()){ pair conf = giter->first; if(conf.first==conf.second) hap_freq_table[conf.first] += geno_freq_table[conf]; else{ hap_freq_table[conf.first] += .5*geno_freq_table[conf]; hap_freq_table[conf.second] += .5*geno_freq_table[conf]; } giter++; } } void IQLBEstimator::update_geno_freq(){ map, double>::iterator giter = geno_freq_table.begin(); while(giter!=geno_freq_table.end()){ giter->second = 0; giter++; } giter = geno_freq_table.begin(); while(giter!=geno_freq_table.end()){ pair conf = giter->first; int k1 = conf.first; int k2 = conf.second; if(k1==k2) giter->second = hap_freq_table[k1]*hap_freq_table[k2]; else giter->second = 2*hap_freq_table[k1]*hap_freq_table[k2]; giter++; } } IQLBEstimator::~IQLBEstimator(){ int row = hap_freq_table.size(); int col = geno_freq_table.size(); for(int k=0;k, double>::iterator iterk; map, double>::iterator iter1; map, double>::iterator iter2; for(int fam=0;famN_typed>0){ for(int i=0;iN;i++){ trio_cat_1 = fam_list[fam]->trio_cat[i]; mis_trio_cat_1 = fam_list[fam]->mis_trio_cat[i]; dim_1 = missing_trio_list[mis_trio_cat_1]->dim; if(dim_1 > 0 && geno_trio_list[trio_cat_1]->info==1){ ind_1 = fam_list[fam]->fam_member[i][0]; par_1 = 0; par_1+=missing_trio_list[mis_trio_cat_1]->typed_par1; par_1+=missing_trio_list[mis_trio_cat_1]->typed_par2; ind_1f = fam_list[fam]->fam_member[i][1]; ind_1m = fam_list[fam]->fam_member[i][2]; for(int j=i+1;jN;j++){ trio_cat_2 = fam_list[fam]->trio_cat[j]; mis_trio_cat_2 = fam_list[fam]->mis_trio_cat[j]; dim_2 = missing_trio_list[mis_trio_cat_2]->dim; if(dim_2 > 0 && geno_trio_list[trio_cat_2]->info==1){ ind_2 = fam_list[fam]->fam_member[j][0]; iterk = (fam_list[fam]->kinship_coeff).find(pair(ind_1,ind_2)); if(iterk == (fam_list[fam]->kinship_coeff).end()){ printf("No IBD coefficient between individual %d and individual %d from family %d. Please check...\n\n", ind_1,ind_2,fam_list[fam]->fam_id); exit(1); } if(iterk->second>1e-8){ // i is related to j par_2 = 0; par_2+=missing_trio_list[mis_trio_cat_2]->typed_par1; par_2+=missing_trio_list[mis_trio_cat_2]->typed_par2; ind_2f = fam_list[fam]->fam_member[j][1]; ind_2m = fam_list[fam]->fam_member[j][2]; kin = iterk->second; iter2 = (fam_list[fam]->IBD2_coeff).find(pair(ind_1,ind_2)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1,ind_2)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(par_1==0 && par_2==0){ if(fam_list[fam]->IBD1_coeff[pair(ind_1,ind_2)]>1e-8){ if(off_one.find(mis_trio_cat_1)==off_one.end()) off_one[mis_trio_cat_1] = 0; if(off_one.find(mis_trio_cat_2)==off_one.end()) off_one[mis_trio_cat_2] = 0; } if(fam_list[fam]->IBD2_coeff[pair(ind_1,ind_2)]>1e-8){ if(off_two.find(mis_trio_cat_1)==off_two.end()) off_two[mis_trio_cat_1] = 0; if(off_two.find(mis_trio_cat_2)==off_two.end()) off_two[mis_trio_cat_2] = 0; } } else if(par_1==0 && (par_2==1 || par_2==2)){ int order = fam_list[fam]->mis_par_order[j]; int ind_2par1 = ind_2f; int ind_2par2 = ind_2m; if(order==1){ ind_2par1 = ind_2m; ind_2par2 = ind_2f; } int diff = 0; int index = relation(fam,ind_1,ind_2,ind_2par1,ind_2par2,kin,IBD_2,IBD_1,diff); if(index==4){ if(off_two.find(mis_trio_cat_1)==off_two.end()) off_two[mis_trio_cat_1] = 0; if(par1_two.find(mis_trio_cat_2)==par1_two.end()) par1_two[mis_trio_cat_2] = 0; } else if(index==5){ if(off_two.find(mis_trio_cat_1)==off_two.end()) off_two[mis_trio_cat_1] = 0; if(par2_two.find(mis_trio_cat_2)==par2_two.end()) par2_two[mis_trio_cat_2] = 0; } else if(index==1){ if(fam_list[fam]->IBD1_coeff[pair(ind_1,ind_2par1)]>1e-8){ if(off_one.find(mis_trio_cat_1)==off_one.end()) off_one[mis_trio_cat_1] = 0; if(par1_one.find(mis_trio_cat_2)==par1_one.end()) par1_one[mis_trio_cat_2] = 0; } if(fam_list[fam]->IBD2_coeff[pair(ind_1,ind_2par1)]>1e-8){ if(off_two.find(mis_trio_cat_1)==off_two.end()) off_two[mis_trio_cat_1] = 0; if(par1_two.find(mis_trio_cat_2)==par1_two.end()) par1_two[mis_trio_cat_2] = 0; } } else if(index==2){ if(fam_list[fam]->IBD1_coeff[pair(ind_1,ind_2par2)]>1e-8){ if(off_one.find(mis_trio_cat_1)==off_one.end()) off_one[mis_trio_cat_1] = 0; if(par2_one.find(mis_trio_cat_2)==par2_one.end()) par2_one[mis_trio_cat_2] = 0; } if(fam_list[fam]->IBD2_coeff[pair(ind_1,ind_2par2)]>1e-8){ if(off_two.find(mis_trio_cat_1)==off_two.end()) off_two[mis_trio_cat_1] = 0; if(par2_two.find(mis_trio_cat_2)==par2_two.end()) par2_two[mis_trio_cat_2] = 0; } } else if(index==3){ if(fam_list[fam]->IBD1_coeff[pair(ind_1,ind_2)]>1e-8){ if(off_one.find(mis_trio_cat_1)==off_one.end()) off_one[mis_trio_cat_1] = 0; if(off_one.find(mis_trio_cat_2)==off_one.end()) off_one[mis_trio_cat_2] = 0; } if(fam_list[fam]->IBD2_coeff[pair(ind_1,ind_2)]>1e-8){ if(off_two.find(mis_trio_cat_1)==off_two.end()) off_two[mis_trio_cat_1] = 0; if(off_two.find(mis_trio_cat_2)==off_two.end()) off_two[mis_trio_cat_2] = 0; } } else if(index==6){ // can not be sibs if(diff==1) printf("please check relationships\n"); if(off_one.find(mis_trio_cat_1)==off_one.end()) off_one[mis_trio_cat_1] = 0; if(par1_one.find(mis_trio_cat_2)==par1_one.end()) par1_one[mis_trio_cat_2] = 0; if(par2_one.find(mis_trio_cat_2)==par2_one.end()) par2_one[mis_trio_cat_2] = 0; } } else if((par_1==1 || par_1==2) && par_2==0){ int order = fam_list[fam]->mis_par_order[i]; int ind_1par1 = ind_1f; int ind_1par2 = ind_1m; if(order==1){ ind_1par1 = ind_1m; ind_1par2 = ind_1f; } int diff = 0; int index = relation(fam,ind_2,ind_1,ind_1par1,ind_1par2,kin,IBD_2,IBD_1,diff); if(index==4){ if(off_two.find(mis_trio_cat_2)==off_two.end()) off_two[mis_trio_cat_2] = 0; if(par1_two.find(mis_trio_cat_1)==par1_two.end()) par1_two[mis_trio_cat_1] = 0; } else if(index==5){ if(off_two.find(mis_trio_cat_2)==off_two.end()) off_two[mis_trio_cat_2] = 0; if(par2_two.find(mis_trio_cat_1)==par2_two.end()) par2_two[mis_trio_cat_1] = 0; } else if(index==1){ if(fam_list[fam]->IBD1_coeff[pair(ind_1par1,ind_2)]>1e-8){ if(off_one.find(mis_trio_cat_2)==off_one.end()) off_one[mis_trio_cat_2] = 0; if(par1_one.find(mis_trio_cat_1)==par1_one.end()) par1_one[mis_trio_cat_1] = 0; } if(fam_list[fam]->IBD2_coeff[pair(ind_1par1,ind_2)]>1e-8){ if(off_two.find(mis_trio_cat_2)==off_two.end()) off_two[mis_trio_cat_2] = 0; if(par1_two.find(mis_trio_cat_1)==par1_two.end()) par1_two[mis_trio_cat_1] = 0; } } else if(index==2){ if(fam_list[fam]->IBD1_coeff[pair(ind_1par2,ind_2)]>1e-8){ if(off_one.find(mis_trio_cat_2)==off_one.end()) off_one[mis_trio_cat_2] = 0; if(par2_one.find(mis_trio_cat_1)==par2_one.end()) par2_one[mis_trio_cat_1] = 0; } if(fam_list[fam]->IBD2_coeff[pair(ind_1par2,ind_2)]>1e-8){ if(off_two.find(mis_trio_cat_2)==off_two.end()) off_two[mis_trio_cat_2] = 0; if(par2_two.find(mis_trio_cat_1)==par2_two.end()) par2_two[mis_trio_cat_1] = 0; } } else if(index==3){ if(fam_list[fam]->IBD1_coeff[pair(ind_1,ind_2)]>1e-8){ if(off_one.find(mis_trio_cat_1)==off_one.end()) off_one[mis_trio_cat_1] = 0; if(off_one.find(mis_trio_cat_2)==off_one.end()) off_one[mis_trio_cat_2] = 0; } if(fam_list[fam]->IBD2_coeff[pair(ind_1,ind_2)]>1e-8){ if(off_two.find(mis_trio_cat_1)==off_two.end()) off_two[mis_trio_cat_1] = 0; if(off_two.find(mis_trio_cat_2)==off_two.end()) off_two[mis_trio_cat_2] = 0; } } else if(index==6){ // can not be sibs if(diff==1) printf("please check relationships\n"); if(off_one.find(mis_trio_cat_2)==off_one.end()) off_one[mis_trio_cat_2] = 0; if(par1_one.find(mis_trio_cat_1)==par1_one.end()) par1_one[mis_trio_cat_1] = 0; if(par2_one.find(mis_trio_cat_1)==par2_one.end()) par2_one[mis_trio_cat_1] = 0; } } else if((par_1==1 || par_1==2) && (par_2==1 || par_2==2)){ int order_1 = fam_list[fam]->mis_par_order[i]; int ind_1par1 = ind_1f; int ind_1par2 = ind_1m; if(order_1==1){ ind_1par1 = ind_1m; ind_1par2 = ind_1f; } int order_2 = fam_list[fam]->mis_par_order[j]; int ind_2par1 = ind_2f; int ind_2par2 = ind_2m; if(order_2==1){ ind_2par1 = ind_2m; ind_2par2 = ind_2f; } int diff = 0; int index = relation(fam,ind_1,ind_1par1,ind_1par2,ind_2,ind_2par1,ind_2par2,kin,IBD_2,IBD_1,diff); if(index==13){ if(off_two.find(mis_trio_cat_1)==off_two.end()) off_two[mis_trio_cat_1] = 0; if(par1_two.find(mis_trio_cat_2)==par1_two.end()) par1_two[mis_trio_cat_2] = 0; } else if(index==14){ if(off_two.find(mis_trio_cat_1)==off_two.end()) off_two[mis_trio_cat_1] = 0; if(par2_two.find(mis_trio_cat_2)==par2_two.end()) par2_two[mis_trio_cat_2] = 0; } else if(index==15){ if(off_two.find(mis_trio_cat_2)==off_two.end()) off_two[mis_trio_cat_2] = 0; if(par1_two.find(mis_trio_cat_1)==par1_two.end()) par1_two[mis_trio_cat_1] = 0; } else if(index==16){ if(off_two.find(mis_trio_cat_2)==off_two.end()) off_two[mis_trio_cat_2] = 0; if(par2_two.find(mis_trio_cat_1)==par2_two.end()) par2_two[mis_trio_cat_1] = 0; } else if(index==17){ if(parents_two.find(mis_trio_cat_1)==parents_two.end()) parents_two[mis_trio_cat_1] = 0; if(parents_two.find(mis_trio_cat_2)==parents_two.end()) parents_two[mis_trio_cat_2] = 0; } else if(index==18){ if(ind_1par1==ind_2par1){ if(par1_two.find(mis_trio_cat_1)==par1_two.end()) par1_two[mis_trio_cat_1] = 0; if(par1_two.find(mis_trio_cat_2)==par1_two.end()) par1_two[mis_trio_cat_2] = 0; } else if(ind_1par1==ind_2par2){ if(par1_two.find(mis_trio_cat_1)==par1_two.end()) par1_two[mis_trio_cat_1] = 0; if(par2_two.find(mis_trio_cat_2)==par2_two.end()) par2_two[mis_trio_cat_2] = 0; } else if(ind_1par2==ind_2par1){ if(par2_two.find(mis_trio_cat_1)==par2_two.end()) par2_two[mis_trio_cat_1] = 0; if(par1_two.find(mis_trio_cat_2)==par1_two.end()) par1_two[mis_trio_cat_2] = 0; } else if(ind_1par2==ind_2par2){ if(par2_two.find(mis_trio_cat_1)==par2_two.end()) par2_two[mis_trio_cat_1] = 0; if(par2_two.find(mis_trio_cat_2)==par2_two.end()) par2_two[mis_trio_cat_2] = 0; } } else if(index==1){ if(fam_list[fam]->IBD1_coeff[pair(ind_1par1,ind_2par1)]>1e-8){ if(par1_one.find(mis_trio_cat_1)==par1_one.end()) par1_one[mis_trio_cat_1] = 0; if(par1_one.find(mis_trio_cat_2)==par1_one.end()) par1_one[mis_trio_cat_2] = 0; } if(fam_list[fam]->IBD2_coeff[pair(ind_1par1,ind_2par1)]>1e-8){ if(par1_two.find(mis_trio_cat_1)==par1_two.end()) par1_two[mis_trio_cat_1] = 0; if(par1_two.find(mis_trio_cat_2)==par1_two.end()) par1_two[mis_trio_cat_2] = 0; } } else if(index==2){ if(fam_list[fam]->IBD1_coeff[pair(ind_1par1,ind_2par2)]>1e-8){ if(par1_one.find(mis_trio_cat_1)==par1_one.end()) par1_one[mis_trio_cat_1] = 0; if(par2_one.find(mis_trio_cat_2)==par2_one.end()) par2_one[mis_trio_cat_2] = 0; } if(fam_list[fam]->IBD2_coeff[pair(ind_1par1,ind_2par2)]>1e-8){ if(par1_two.find(mis_trio_cat_1)==par1_two.end()) par1_two[mis_trio_cat_1] = 0; if(par2_two.find(mis_trio_cat_2)==par2_two.end()) par2_two[mis_trio_cat_2] = 0; } } else if(index==3){ if(fam_list[fam]->IBD1_coeff[pair(ind_1par2,ind_2par1)]>1e-8){ if(par2_one.find(mis_trio_cat_1)==par2_one.end()) par2_one[mis_trio_cat_1] = 0; if(par1_one.find(mis_trio_cat_2)==par1_one.end()) par1_one[mis_trio_cat_2] = 0; } if(fam_list[fam]->IBD2_coeff[pair(ind_1par2,ind_2par1)]>1e-8){ if(par2_two.find(mis_trio_cat_1)==par2_two.end()) par2_two[mis_trio_cat_1] = 0; if(par1_two.find(mis_trio_cat_2)==par1_two.end()) par1_two[mis_trio_cat_2] = 0; } } else if(index==4){ if(fam_list[fam]->IBD1_coeff[pair(ind_1par2,ind_2par2)]>1e-8){ if(par2_one.find(mis_trio_cat_1)==par2_one.end()) par2_one[mis_trio_cat_1] = 0; if(par2_one.find(mis_trio_cat_2)==par2_one.end()) par2_one[mis_trio_cat_2] = 0; } if(fam_list[fam]->IBD2_coeff[pair(ind_1par2,ind_2par2)]>1e-8){ if(par2_two.find(mis_trio_cat_1)==par2_two.end()) par2_two[mis_trio_cat_1] = 0; if(par2_two.find(mis_trio_cat_2)==par2_two.end()) par2_two[mis_trio_cat_2] = 0; } } else if(index==5){ if(fam_list[fam]->IBD1_coeff[pair(ind_1par1,ind_2)]>1e-8){ if(par1_one.find(mis_trio_cat_1)==par1_one.end()) par1_one[mis_trio_cat_1] = 0; if(off_one.find(mis_trio_cat_2)==off_one.end()) off_one[mis_trio_cat_2] = 0; } if(fam_list[fam]->IBD2_coeff[pair(ind_1par1,ind_2)]>1e-8){ if(par1_two.find(mis_trio_cat_1)==par1_two.end()) par1_two[mis_trio_cat_1] = 0; if(off_two.find(mis_trio_cat_2)==off_two.end()) off_two[mis_trio_cat_2] = 0; } } else if(index==6){ if(fam_list[fam]->IBD1_coeff[pair(ind_1par2,ind_2)]>1e-8){ if(par2_one.find(mis_trio_cat_1)==par2_one.end()) par2_one[mis_trio_cat_1] = 0; if(off_one.find(mis_trio_cat_2)==off_one.end()) off_one[mis_trio_cat_2] = 0; } if(fam_list[fam]->IBD2_coeff[pair(ind_1par2,ind_2)]>1e-8){ if(par2_two.find(mis_trio_cat_1)==par2_two.end()) par2_two[mis_trio_cat_1] = 0; if(off_two.find(mis_trio_cat_2)==off_two.end()) off_two[mis_trio_cat_2] = 0; } } else if(index==7){ if(fam_list[fam]->IBD1_coeff[pair(ind_1,ind_2par1)]>1e-8){ if(off_one.find(mis_trio_cat_1)==off_one.end()) off_one[mis_trio_cat_1] = 0; if(par1_one.find(mis_trio_cat_2)==par1_one.end()) par1_one[mis_trio_cat_2] = 0; } if(fam_list[fam]->IBD2_coeff[pair(ind_1,ind_2par1)]>1e-8){ if(off_two.find(mis_trio_cat_1)==off_two.end()) off_two[mis_trio_cat_1] = 0; if(par1_two.find(mis_trio_cat_2)==par1_two.end()) par1_two[mis_trio_cat_2] = 0; } } else if(index==8){ if(fam_list[fam]->IBD1_coeff[pair(ind_1,ind_2par2)]>1e-8){ if(off_one.find(mis_trio_cat_1)==off_one.end()) off_one[mis_trio_cat_1] = 0; if(par2_one.find(mis_trio_cat_2)==par2_one.end()) par2_one[mis_trio_cat_2] = 0; } if(fam_list[fam]->IBD2_coeff[pair(ind_1,ind_2par2)]>1e-8){ if(off_two.find(mis_trio_cat_1)==off_two.end()) off_two[mis_trio_cat_1] = 0; if(par2_two.find(mis_trio_cat_2)==par2_two.end()) par2_two[mis_trio_cat_2] = 0; } } else if(index==9){ if(diff==1){ // j and if are sibs if(par1_two.find(mis_trio_cat_1)==par1_two.end()) par1_two[mis_trio_cat_1] = 0; if(parents_one.find(mis_trio_cat_2)==parents_one.end()) parents_one[mis_trio_cat_2] = 0; } else{ if(par1_one.find(mis_trio_cat_1)==par1_one.end()) par1_one[mis_trio_cat_1] = 0; if(par1_one.find(mis_trio_cat_2)==par1_one.end()) par1_one[mis_trio_cat_2] = 0; if(par2_one.find(mis_trio_cat_2)==par2_one.end()) par2_one[mis_trio_cat_2] = 0; } } else if(index==10){ if(diff==1){ // j and im are sibs if(par2_two.find(mis_trio_cat_1)==par2_two.end()) par2_two[mis_trio_cat_1] = 0; if(parents_one.find(mis_trio_cat_2)==parents_one.end()) parents_one[mis_trio_cat_2] = 0; } else{ if(par2_one.find(mis_trio_cat_1)==par2_one.end()) par2_one[mis_trio_cat_1] = 0; if(par1_one.find(mis_trio_cat_2)==par1_one.end()) par1_one[mis_trio_cat_2] = 0; if(par2_one.find(mis_trio_cat_2)==par2_one.end()) par2_one[mis_trio_cat_2] = 0; } } else if(index==11){ if(diff==1){ // i and jf are sibs if(par1_two.find(mis_trio_cat_2)==par1_two.end()) par1_two[mis_trio_cat_2] = 0; if(parents_one.find(mis_trio_cat_1)==parents_one.end()) parents_one[mis_trio_cat_1] = 0; } else{ if(par1_one.find(mis_trio_cat_2)==par1_one.end()) par1_one[mis_trio_cat_2] = 0; if(par1_one.find(mis_trio_cat_1)==par1_one.end()) par1_one[mis_trio_cat_1] = 0; if(par2_one.find(mis_trio_cat_1)==par2_one.end()) par2_one[mis_trio_cat_1] = 0; } } else if(index==12){ if(diff==1){ // i and jm are sibs if(par2_two.find(mis_trio_cat_2)==par2_two.end()) par2_two[mis_trio_cat_2] = 0; if(parents_one.find(mis_trio_cat_1)==parents_one.end()) parents_one[mis_trio_cat_1] = 0; } else{ if(par2_one.find(mis_trio_cat_2)==par2_one.end()) par2_one[mis_trio_cat_2] = 0; if(par1_one.find(mis_trio_cat_1)==par1_one.end()) par1_one[mis_trio_cat_1] = 0; if(par2_one.find(mis_trio_cat_1)==par2_one.end()) par2_one[mis_trio_cat_1] = 0; } } } } // i, j related } // j is typed } // complete j } // i is typed } // complete i } } } void IQLBEstimator::allocate(){ map, double>::iterator giter=geno_freq_table.begin(); int loc = 0; while(giter!=geno_freq_table.end()){ hap_pair_map[giter->first] = loc; loc++; giter++; } int h_size = hap_freq_table.size(); int col = geno_freq_table.size(); rst = new double*[h_size-1]; for(int i=0;i0){ result_par1_one = new double **[par1_one_size]; map::iterator pi=par1_one.begin(); for(int i=0;isecond = i; int mis = pi->first; int dim_1 = missing_trio_list[mis]->dim; result_par1_one[i] = new double*[dim_1]; for(int j=0;j0){ result_par1_two = new double **[par1_two_size]; map::iterator pi=par1_two.begin(); for(int i=0;isecond = i; int mis = pi->first; int dim_1 = missing_trio_list[mis]->dim; result_par1_two[i] = new double*[dim_1]; for(int j=0;j0){ result_par2_one = new double **[par2_one_size]; map::iterator pi=par2_one.begin(); for(int i=0;isecond = i; int mis = pi->first; int dim_1 = missing_trio_list[mis]->dim; result_par2_one[i] = new double*[dim_1]; for(int j=0;j0){ result_par2_two = new double **[par2_two_size]; map::iterator pi=par2_two.begin(); for(int i=0;isecond = i; int mis = pi->first; int dim_1 = missing_trio_list[mis]->dim; result_par2_two[i] = new double*[dim_1]; for(int j=0;j0){ result_off_one = new double **[off_one_size]; map::iterator pi=off_one.begin(); for(int i=0;isecond = i; int mis = pi->first; int dim_1 = missing_trio_list[mis]->dim; result_off_one[i] = new double*[dim_1]; for(int j=0;j0){ result_off_two = new double **[off_two_size]; map::iterator pi=off_two.begin(); for(int i=0;isecond = i; int mis = pi->first; int dim_1 = missing_trio_list[mis]->dim; result_off_two[i] = new double*[dim_1]; for(int j=0;j0){ result_parents_one = new double ***[parents_one_size]; map::iterator pi=parents_one.begin(); for(int i=0;isecond = i; int mis = pi->first; int dim_1 = missing_trio_list[mis]->dim; result_parents_one[i] = new double**[dim_1]; for(int j=0;j0){ result_parents_two = new double ***[parents_two_size]; map::iterator pi=parents_two.begin(); for(int i=0;isecond = i; int mis = pi->first; int dim_1 = missing_trio_list[mis]->dim; result_parents_two[i] = new double**[dim_1]; for(int j=0;jN_typed>0){ int trio_cat; int mis_trio_cat; int dim=0; for(int i=0;iN;i++){ trio_cat = fam_list[fam]->trio_cat[i]; if(geno_trio_list[trio_cat]->info==1){ mis_trio_cat = fam_list[fam]->mis_trio_cat[i]; dim+=missing_trio_list[mis_trio_cat]->dim; } } if(dim>tot_dim) tot_dim = dim; } } cholent = new double*[tot_dim]; for(int i=0;i0){ map::iterator pi=par1_one.begin(); for(int i=0;ifirst]->dim; for(int j=0;j0){ map::iterator pi=par1_two.begin(); for(int i=0;ifirst]->dim; for(int j=0;j0){ map::iterator pi=par2_one.begin(); for(int i=0;ifirst]->dim; for(int j=0;j0){ map::iterator pi=par2_two.begin(); for(int i=0;ifirst]->dim; for(int j=0;j0){ map::iterator pi=off_one.begin(); for(int i=0;ifirst]->dim; for(int j=0;j0){ map::iterator pi=off_two.begin(); for(int i=0;ifirst]->dim; for(int j=0;j0){ map::iterator pi=parents_one.begin(); for(int i=0;ifirst]->dim; for(int j=0;j0){ map::iterator pi=parents_two.begin(); for(int i=0;ifirst]->dim; for(int j=0;j::iterator iter=hap_freq_table.begin(); while(iter!=hap_freq_table.end()){ freq[loc] = iter->second; loc++; iter++; } loc = 0; map, double>::iterator giter=geno_freq_table.begin(); while(giter!=geno_freq_table.end()){ genofreq[loc] = giter->second; loc++; giter++; } int h_last=0; iter = hap_freq_table.end(); iter--; h_last=iter->first; giter=geno_freq_table.begin(); while(giter!=geno_freq_table.end()){ pair conf = giter->first; loc = 0; iter = hap_freq_table.begin(); while(locfirst; if(conf.first == conf.second && conf.first == h) rst[loc][hap_pair_map[conf]] = 2*hap_freq_table[h]; else if(conf.first == conf.second && conf.first == h_last) rst[loc][hap_pair_map[conf]] = -2*hap_freq_table[h_last]; else if(conf.first < conf.second && conf.first == h && conf.second != h_last) rst[loc][hap_pair_map[conf]] = 2*hap_freq_table[conf.second]; else if(conf.first < conf.second && conf.second == h) rst[loc][hap_pair_map[conf]] = 2*hap_freq_table[conf.first]; else if(conf.first == h && conf.second == h_last) rst[loc][hap_pair_map[conf]] = 2*(hap_freq_table[h_last]-hap_freq_table[h]); else if(conf.first != h && conf.first != h_last && conf.second == h_last) rst[loc][hap_pair_map[conf]] = -2*hap_freq_table[conf.first]; loc++; iter++; } giter++; } for(int i=0;iget_num()>0 && geno_trio_list[i]->info==1 && missing_trio_list[geno_missing_trio_map[i]]->dim>0) geno_trio_list[i]->expected_hap_num(hap_freq_table,freq,geno_freq_table,hap_pair_map,count,missing_trio_list[geno_missing_trio_map[i]]->sites); } for(int i=0;iget_num()>0 && missing_trio_list[i]->dim>0){ missing_trio_list[i]->set_weight(hap_freq_table,geno_freq_table,hap_pair_map,count); missing_trio_list[i]->compute_centmatrix(hap_freq_table,freq,geno_freq_table); missing_trio_list[i]->compute_derivative(hap_freq_table,geno_freq_table,hap_pair_map,count,rst,hst); missing_trio_list[i]->expected_hap_num(hap_freq_table,freq,geno_freq_table,hap_pair_map,count); } } int par1_one_size = par1_one.size(); if(par1_one_size>0){ map::iterator pi=par1_one.begin(); for(int i=0;isecond = i; int mis = pi->first; int dim_1 = missing_trio_list[mis]->dim; for(int j=0;jcond_par1_one(result_par1_one[i],hap_freq_table,geno_freq_table,hap_map); pi++; } } int par1_two_size = par1_two.size(); if(par1_two_size>0){ map::iterator pi=par1_two.begin(); for(int i=0;isecond = i; int mis = pi->first; int dim_1 = missing_trio_list[mis]->dim; for(int j=0;jcond_par1_two(result_par1_two[i],hap_freq_table,geno_freq_table,hap_pair_map); pi++; } } int par2_one_size = par2_one.size(); if(par2_one_size>0){ map::iterator pi=par2_one.begin(); for(int i=0;isecond = i; int mis = pi->first; int dim_1 = missing_trio_list[mis]->dim; for(int j=0;jcond_par2_one(result_par2_one[i],hap_freq_table,geno_freq_table,hap_map); pi++; } } int par2_two_size = par2_two.size(); if(par2_two_size>0){ map::iterator pi=par2_two.begin(); for(int i=0;isecond = i; int mis = pi->first; int dim_1 = missing_trio_list[mis]->dim; for(int j=0;jcond_par2_two(result_par2_two[i],hap_freq_table,geno_freq_table,hap_pair_map); pi++; } } int off_one_size = off_one.size(); if(off_one_size>0){ map::iterator pi=off_one.begin(); for(int i=0;isecond = i; int mis = pi->first; int dim_1 = missing_trio_list[mis]->dim; for(int j=0;jcond_off_one(result_off_one[i],hap_freq_table,geno_freq_table,hap_map); pi++; } } int off_two_size = off_two.size(); if(off_two_size>0){ map::iterator pi=off_two.begin(); for(int i=0;isecond = i; int mis = pi->first; int dim_1 = missing_trio_list[mis]->dim; for(int j=0;jcond_off_two(result_off_two[i],hap_freq_table,geno_freq_table,hap_pair_map); pi++; } } int parents_one_size = parents_one.size(); if(parents_one_size>0){ map::iterator pi=parents_one.begin(); for(int i=0;isecond = i; int mis = pi->first; int dim_1 = missing_trio_list[mis]->dim; for(int j=0;jcond_parents_one(result_parents_one[i],hap_freq_table,geno_freq_table,hap_map); pi++; } } int parents_two_size = parents_two.size(); if(parents_two_size>0){ map::iterator pi=parents_two.begin(); for(int i=0;isecond = i; int mis = pi->first; int dim_1 = missing_trio_list[mis]->dim; for(int j=0;jcond_parents_two(result_parents_two[i],hap_freq_table,geno_freq_table,hap_pair_map); pi++; } } for(int i=0;iN_typed>0){ int trio_cat; int mis_trio_cat; int dim=0; for(int i=0;iN;i++){ trio_cat = fam_list[fam]->trio_cat[i]; if(geno_trio_list[trio_cat]->info==1){ mis_trio_cat = fam_list[fam]->mis_trio_cat[i]; dim+=missing_trio_list[mis_trio_cat]->dim; } } if(dim<=0) continue; for(int i=0;i, double>::iterator iterk; map, double>::iterator iter1; map, double>::iterator iter2; double value; flag_pd = 1; for(int i=0;iN;i++){ trio_cat_1 = fam_list[fam]->trio_cat[i]; mis_trio_cat_1 = fam_list[fam]->mis_trio_cat[i]; dim_1 = missing_trio_list[mis_trio_cat_1]->dim; if(dim_1 > 0 && geno_trio_list[trio_cat_1]->info==1){ nb2 = nb1; ind_1 = fam_list[fam]->fam_member[i][0]; ind_2 = ind_1; iterk = (fam_list[fam]->kinship_coeff).find(pair(ind_1,ind_2)); if(iterk == (fam_list[fam]->kinship_coeff).end()){ printf("No IBD coefficient for individual %d from family %d. Please check...\n\n", ind_1,fam_list[fam]->fam_id); exit(1); } for(int j=0;jcentmatrix[j][k]; nb2+=dim_1; par_1 = 0; par_1+=missing_trio_list[mis_trio_cat_1]->typed_par1; par_1+=missing_trio_list[mis_trio_cat_1]->typed_par2; ind_1f = fam_list[fam]->fam_member[i][1]; ind_1m = fam_list[fam]->fam_member[i][2]; for(int j=i+1;jN;j++){ trio_cat_2 = fam_list[fam]->trio_cat[j]; mis_trio_cat_2 = fam_list[fam]->mis_trio_cat[j]; dim_2 = missing_trio_list[mis_trio_cat_2]->dim; if(dim_2 > 0 && geno_trio_list[trio_cat_2]->info==1){ ind_2 = fam_list[fam]->fam_member[j][0]; iterk = (fam_list[fam]->kinship_coeff).find(pair(ind_1,ind_2)); if(iterk == (fam_list[fam]->kinship_coeff).end()){ printf("No IBD coefficient between individual %d and individual %d from family %d. Please check...\n\n", ind_1,ind_2,fam_list[fam]->fam_id); exit(1); } if(iterk->second>1e-8){ // i is related to j par_2 = 0; par_2+=missing_trio_list[mis_trio_cat_2]->typed_par1; par_2+=missing_trio_list[mis_trio_cat_2]->typed_par2; ind_2f = fam_list[fam]->fam_member[j][1]; ind_2m = fam_list[fam]->fam_member[j][2]; kin = iterk->second; iter2 = (fam_list[fam]->IBD2_coeff).find(pair(ind_1,ind_2)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1,ind_2)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(par_1==0 && par_2==0){ if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;kmis_par_order[j]; int ind_2par1 = ind_2f; int ind_2par2 = ind_2m; if(order==1){ ind_2par1 = ind_2m; ind_2par2 = ind_2f; } int diff = 0; int index = relation(fam,ind_1,ind_2,ind_2par1,ind_2par2,kin,IBD_2,IBD_1,diff); if(index==4){ for(int k=0;kIBD2_coeff).find(pair(ind_1,ind_2par1)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1,ind_2par1)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;kIBD2_coeff).find(pair(ind_1,ind_2par2)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1,ind_2par2)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;k1e-8){ for(int k=0;k1e-8){ for(int k=0;kIBD1_coeff).find(pair(ind_1,ind_2par1)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;kIBD1_coeff).find(pair(ind_1,ind_2par2)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;kmis_par_order[i]; int ind_1par1 = ind_1f; int ind_1par2 = ind_1m; if(order==1){ ind_1par1 = ind_1m; ind_1par2 = ind_1f; } int diff = 0; int index = relation(fam,ind_2,ind_1,ind_1par1,ind_1par2,kin,IBD_2,IBD_1,diff); if(index==4){ for(int k=0;kIBD2_coeff).find(pair(ind_1par1,ind_2)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1par1,ind_2)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;kIBD2_coeff).find(pair(ind_1par2,ind_2)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1par2,ind_2)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;k1e-8){ for(int k=0;k1e-8){ for(int k=0;kIBD1_coeff).find(pair(ind_1par1,ind_2)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;kIBD1_coeff).find(pair(ind_1par2,ind_2)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;kmis_par_order[i]; int ind_1par1 = ind_1f; int ind_1par2 = ind_1m; if(order_1==1){ ind_1par1 = ind_1m; ind_1par2 = ind_1f; } int order_2 = fam_list[fam]->mis_par_order[j]; int ind_2par1 = ind_2f; int ind_2par2 = ind_2m; if(order_2==1){ ind_2par1 = ind_2m; ind_2par2 = ind_2f; } int diff = 0; int index = relation(fam,ind_1,ind_1par1,ind_1par2,ind_2,ind_2par1,ind_2par2,kin,IBD_2,IBD_1,diff); if(index==13){ for(int k=0;kIBD2_coeff).find(pair(ind_1par1,ind_2par1)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1par1,ind_2par1)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;kIBD2_coeff).find(pair(ind_1par1,ind_2par2)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1par1,ind_2par2)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;kIBD2_coeff).find(pair(ind_1par2,ind_2par1)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1par2,ind_2par1)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;kIBD2_coeff).find(pair(ind_1par2,ind_2par2)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1par2,ind_2par2)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;kIBD2_coeff).find(pair(ind_1par1,ind_2)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1par1,ind_2)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;kIBD2_coeff).find(pair(ind_1par2,ind_2)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1par2,ind_2)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;kIBD2_coeff).find(pair(ind_1,ind_2par1)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1,ind_2par1)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;kIBD2_coeff).find(pair(ind_1,ind_2par2)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1,ind_2par2)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;kfirst).first; int h2 = (giter->first).second; if(h1==h2){ for(int k=0;kIBD1_coeff).find(pair(ind_1par1,ind_2par1)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;kIBD1_coeff).find(pair(ind_1par1,ind_2par2)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;kfirst).first; int h2 = (giter->first).second; if(h1==h2){ for(int k=0;kIBD1_coeff).find(pair(ind_1par2,ind_2par1)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;kIBD1_coeff).find(pair(ind_1par2,ind_2par2)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;kfirst).first; int h2 = (giter->first).second; if(h1==h2){ for(int k=0;kIBD1_coeff).find(pair(ind_1par1,ind_2par1)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;kIBD1_coeff).find(pair(ind_1par2,ind_2par1)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;kfirst).first; int h2 = (giter->first).second; if(h1==h2){ for(int k=0;kIBD1_coeff).find(pair(ind_1par1,ind_2par2)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;kIBD1_coeff).find(pair(ind_1par2,ind_2par2)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k hap_num = geno_trio_list[trio_cat_1]->get_hap_num(); for(int j=0;jderivative[j][k]; cholent[nb1+j][h_size-1] = hap_num[j]; } nb1+=dim_1; } // i is typed } // complete i if(cholesky(covMatrix,dim,cholent,h_size,chol,cholaug,0) != 1){ printf("cholesky decomposition of the cov matrix failed for family %d. \n", fam_list[fam]->fam_id); flag_pd = 0; exit(1); } if(flag_pd==1){ double **s = self_colmultiply(cholaug,dim,1,h_size-1,1,h_size-1); double **t = self_colmultiply(cholaug,dim,1,h_size-1,h_size,h_size); for(int i=0;isecond+=s[i][0]; n+=iter->second; iter++; } iter->second = 1-n; clear_matrix(s,h_size-1,1); } else{ double var; double &determ = var; double **v = inverse(denominator,h_size-1,determ); double **w = new double*[h_size-1]; for(int i=0;isecond+=s[i][0]; n+=iter->second; iter++; } iter->second = 1-n; clear_matrix(s,h_size-1,1); clear_matrix(v,h_size-1,h_size-1); clear_matrix(w,h_size-1,1); } clear_matrix(cholentmain,h_size-1,h_size); clear_matrix(cholmain,h_size-1,h_size-1); clear_matrix(cholaugmain,h_size-1,h_size); update_geno_freq(); } void IQLBEstimator::Score_test_new(){ int h_size = hap_freq_table.size(); int col = geno_freq_table.size(); for(int i=0;i::iterator iter=hap_freq_table.begin(); while(iter!=hap_freq_table.end()){ freq[loc] = iter->second; loc++; iter++; } loc = 0; map, double>::iterator giter=geno_freq_table.begin(); while(giter!=geno_freq_table.end()){ genofreq[loc] = giter->second; loc++; giter++; } int h_last=0; iter = hap_freq_table.end(); iter--; h_last=iter->first; giter=geno_freq_table.begin(); while(giter!=geno_freq_table.end()){ pair conf = giter->first; loc = 0; iter = hap_freq_table.begin(); while(locfirst; if(conf.first == conf.second && conf.first == h) rst[loc][hap_pair_map[conf]] = 2*hap_freq_table[h]; else if(conf.first == conf.second && conf.first == h_last) rst[loc][hap_pair_map[conf]] = -2*hap_freq_table[h_last]; else if(conf.first < conf.second && conf.first == h && conf.second != h_last) rst[loc][hap_pair_map[conf]] = 2*hap_freq_table[conf.second]; else if(conf.first < conf.second && conf.second == h) rst[loc][hap_pair_map[conf]] = 2*hap_freq_table[conf.first]; else if(conf.first == h && conf.second == h_last) rst[loc][hap_pair_map[conf]] = 2*(hap_freq_table[h_last]-hap_freq_table[h]); else if(conf.first != h && conf.first != h_last && conf.second == h_last) rst[loc][hap_pair_map[conf]] = -2*hap_freq_table[conf.first]; loc++; iter++; } giter++; } double ***drst; if(oneparent==1){ drst = new double**[h_size-1]; for(int i=0;i, double>::iterator giter1; map, double>::iterator giter2; loc1 = 0; giter1=geno_freq_table.begin(); while(giter1!=geno_freq_table.end()){ int i=(giter1->first).first; int j=(giter1->first).second; value1 = 1.0/hap_freq_table[i]; value2 = 1.0/hap_freq_table[j]; value3 = 1.0/(hap_freq_table[i]+hap_freq_table[j]); loc2 = 0; giter2=geno_freq_table.begin(); while(giter2!=geno_freq_table.end()){ int k=(giter2->first).first; int l=(giter2->first).second; if(k==l){ // parent is homozygote if(i==k){ loc = hap_map[j]; if(loc!=h_size-1) drst[loc][loc1][loc2]+=value2; else{ for(int s=0;sget_num()>0 && geno_trio_list[i]->info==1 && missing_trio_list[geno_missing_trio_map[i]]->dim>0) geno_trio_list[i]->expected_hap_num(hap_freq_table,freq,geno_freq_table,hap_pair_map,count,missing_trio_list[geno_missing_trio_map[i]]->sites); } for(int i=0;iget_num()>0 && missing_trio_list[i]->dim>0){ missing_trio_list[i]->set_weight(hap_freq_table,geno_freq_table,hap_pair_map,count); missing_trio_list[i]->compute_centmatrix(hap_freq_table,freq,geno_freq_table); missing_trio_list[i]->compute_derivative(hap_freq_table,geno_freq_table,hap_pair_map,count,rst,hst); missing_trio_list[i]->expected_hap_num(hap_freq_table,freq,geno_freq_table,hap_pair_map,count); if(missing_trio_list[i]->typed_par1==1 || missing_trio_list[i]->typed_par2==1) missing_trio_list[i]->compute_deriv1to3(hap_freq_table,geno_freq_table,hap_pair_map,count,rst,drst); } } if(oneparent==1){ for(int i=0;i0){ map::iterator pi=par1_one.begin(); for(int i=0;isecond = i; int mis = pi->first; int dim_1 = missing_trio_list[mis]->dim; for(int j=0;jcond_par1_one(result_par1_one[i],hap_freq_table,geno_freq_table,hap_map); pi++; } } int par1_two_size = par1_two.size(); if(par1_two_size>0){ map::iterator pi=par1_two.begin(); for(int i=0;isecond = i; int mis = pi->first; int dim_1 = missing_trio_list[mis]->dim; for(int j=0;jcond_par1_two(result_par1_two[i],hap_freq_table,geno_freq_table,hap_pair_map); pi++; } } int par2_one_size = par2_one.size(); if(par2_one_size>0){ map::iterator pi=par2_one.begin(); for(int i=0;isecond = i; int mis = pi->first; int dim_1 = missing_trio_list[mis]->dim; for(int j=0;jcond_par2_one(result_par2_one[i],hap_freq_table,geno_freq_table,hap_map); pi++; } } int par2_two_size = par2_two.size(); if(par2_two_size>0){ map::iterator pi=par2_two.begin(); for(int i=0;isecond = i; int mis = pi->first; int dim_1 = missing_trio_list[mis]->dim; for(int j=0;jcond_par2_two(result_par2_two[i],hap_freq_table,geno_freq_table,hap_pair_map); pi++; } } int off_one_size = off_one.size(); if(off_one_size>0){ map::iterator pi=off_one.begin(); for(int i=0;isecond = i; int mis = pi->first; int dim_1 = missing_trio_list[mis]->dim; for(int j=0;jcond_off_one(result_off_one[i],hap_freq_table,geno_freq_table,hap_map); pi++; } } int off_two_size = off_two.size(); if(off_two_size>0){ map::iterator pi=off_two.begin(); for(int i=0;isecond = i; int mis = pi->first; int dim_1 = missing_trio_list[mis]->dim; for(int j=0;jcond_off_two(result_off_two[i],hap_freq_table,geno_freq_table,hap_pair_map); pi++; } } int parents_one_size = parents_one.size(); if(parents_one_size>0){ map::iterator pi=parents_one.begin(); for(int i=0;isecond = i; int mis = pi->first; int dim_1 = missing_trio_list[mis]->dim; for(int j=0;jcond_parents_one(result_parents_one[i],hap_freq_table,geno_freq_table,hap_map); pi++; } } int parents_two_size = parents_two.size(); if(parents_two_size>0){ map::iterator pi=parents_two.begin(); for(int i=0;isecond = i; int mis = pi->first; int dim_1 = missing_trio_list[mis]->dim; for(int j=0;jcond_parents_two(result_parents_two[i],hap_freq_table,geno_freq_table,hap_pair_map); pi++; } } for(int i=0;iN_typed>0){ int trio_cat; int mis_trio_cat; int dim=0; for(int i=0;iN;i++){ trio_cat = fam_list[fam]->trio_cat[i]; if(geno_trio_list[trio_cat]->info==1){ mis_trio_cat = fam_list[fam]->mis_trio_cat[i]; dim+=missing_trio_list[mis_trio_cat]->dim; } } if(dim<=0) continue; for(int i=0;i, double>::iterator iterk; map, double>::iterator iter1; map, double>::iterator iter2; double value; double test_weight; double test_trio_weight1, test_trio_weight2, test_trio_weight3; flag_pd = 1; for(int i=0;iN;i++){ trio_cat_1 = fam_list[fam]->trio_cat[i]; mis_trio_cat_1 = fam_list[fam]->mis_trio_cat[i]; dim_1 = missing_trio_list[mis_trio_cat_1]->dim; if(dim_1 > 0 && geno_trio_list[trio_cat_1]->info==1){ nb2 = nb1; ind_1 = fam_list[fam]->fam_member[i][0]; ind_2 = ind_1; iterk = (fam_list[fam]->kinship_coeff).find(pair(ind_1,ind_2)); if(iterk == (fam_list[fam]->kinship_coeff).end()){ printf("No IBD coefficient for individual %d from family %d. Please check...\n\n", ind_1,fam_list[fam]->fam_id); exit(1); } for(int j=0;jcentmatrix[j][k]; nb2+=dim_1; par_1 = 0; par_1+=missing_trio_list[mis_trio_cat_1]->typed_par1; par_1+=missing_trio_list[mis_trio_cat_1]->typed_par2; ind_1f = fam_list[fam]->fam_member[i][1]; ind_1m = fam_list[fam]->fam_member[i][2]; for(int j=i+1;jN;j++){ trio_cat_2 = fam_list[fam]->trio_cat[j]; mis_trio_cat_2 = fam_list[fam]->mis_trio_cat[j]; dim_2 = missing_trio_list[mis_trio_cat_2]->dim; if(dim_2 > 0 && geno_trio_list[trio_cat_2]->info==1){ ind_2 = fam_list[fam]->fam_member[j][0]; iterk = (fam_list[fam]->kinship_coeff).find(pair(ind_1,ind_2)); if(iterk == (fam_list[fam]->kinship_coeff).end()){ printf("No IBD coefficient between individual %d and individual %d from family %d. Please check...\n\n", ind_1,ind_2,fam_list[fam]->fam_id); exit(1); } if(iterk->second>1e-8){ // i is related to j par_2 = 0; par_2+=missing_trio_list[mis_trio_cat_2]->typed_par1; par_2+=missing_trio_list[mis_trio_cat_2]->typed_par2; ind_2f = fam_list[fam]->fam_member[j][1]; ind_2m = fam_list[fam]->fam_member[j][2]; kin = iterk->second; iter2 = (fam_list[fam]->IBD2_coeff).find(pair(ind_1,ind_2)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1,ind_2)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(par_1==0 && par_2==0){ if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;kmis_par_order[j]; int ind_2par1 = ind_2f; int ind_2par2 = ind_2m; if(order==1){ ind_2par1 = ind_2m; ind_2par2 = ind_2f; } int diff = 0; int index = relation(fam,ind_1,ind_2,ind_2par1,ind_2par2,kin,IBD_2,IBD_1,diff); if(index==4){ for(int k=0;kIBD2_coeff).find(pair(ind_1,ind_2par1)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1,ind_2par1)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;kIBD2_coeff).find(pair(ind_1,ind_2par2)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1,ind_2par2)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;k1e-8){ for(int k=0;k1e-8){ for(int k=0;kIBD1_coeff).find(pair(ind_1,ind_2par1)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;kIBD1_coeff).find(pair(ind_1,ind_2par2)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;kmis_par_order[i]; int ind_1par1 = ind_1f; int ind_1par2 = ind_1m; if(order==1){ ind_1par1 = ind_1m; ind_1par2 = ind_1f; } int diff = 0; int index = relation(fam,ind_2,ind_1,ind_1par1,ind_1par2,kin,IBD_2,IBD_1,diff); if(index==4){ for(int k=0;kIBD2_coeff).find(pair(ind_1par1,ind_2)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1par1,ind_2)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;kIBD2_coeff).find(pair(ind_1par2,ind_2)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1par2,ind_2)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;k1e-8){ for(int k=0;k1e-8){ for(int k=0;kIBD1_coeff).find(pair(ind_1par1,ind_2)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;kIBD1_coeff).find(pair(ind_1par2,ind_2)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;kmis_par_order[i]; int ind_1par1 = ind_1f; int ind_1par2 = ind_1m; if(order_1==1){ ind_1par1 = ind_1m; ind_1par2 = ind_1f; } int order_2 = fam_list[fam]->mis_par_order[j]; int ind_2par1 = ind_2f; int ind_2par2 = ind_2m; if(order_2==1){ ind_2par1 = ind_2m; ind_2par2 = ind_2f; } int diff = 0; int index = relation(fam,ind_1,ind_1par1,ind_1par2,ind_2,ind_2par1,ind_2par2,kin,IBD_2,IBD_1,diff); if(index==13){ for(int k=0;kIBD2_coeff).find(pair(ind_1par1,ind_2par1)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1par1,ind_2par1)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;kIBD2_coeff).find(pair(ind_1par1,ind_2par2)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1par1,ind_2par2)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;kIBD2_coeff).find(pair(ind_1par2,ind_2par1)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1par2,ind_2par1)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;kIBD2_coeff).find(pair(ind_1par2,ind_2par2)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1par2,ind_2par2)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;kIBD2_coeff).find(pair(ind_1par1,ind_2)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1par1,ind_2)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;kIBD2_coeff).find(pair(ind_1par2,ind_2)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1par2,ind_2)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;kIBD2_coeff).find(pair(ind_1,ind_2par1)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1,ind_2par1)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;kIBD2_coeff).find(pair(ind_1,ind_2par2)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1,ind_2par2)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;kfirst).first; int h2 = (giter->first).second; if(h1==h2){ for(int k=0;kIBD1_coeff).find(pair(ind_1par1,ind_2par1)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;kIBD1_coeff).find(pair(ind_1par1,ind_2par2)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;kfirst).first; int h2 = (giter->first).second; if(h1==h2){ for(int k=0;kIBD1_coeff).find(pair(ind_1par2,ind_2par1)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;kIBD1_coeff).find(pair(ind_1par2,ind_2par2)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;kfirst).first; int h2 = (giter->first).second; if(h1==h2){ for(int k=0;kIBD1_coeff).find(pair(ind_1par1,ind_2par1)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;kIBD1_coeff).find(pair(ind_1par2,ind_2par1)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;kfirst).first; int h2 = (giter->first).second; if(h1==h2){ for(int k=0;kIBD1_coeff).find(pair(ind_1par1,ind_2par2)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;kIBD1_coeff).find(pair(ind_1par2,ind_2par2)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k hap_num = geno_trio_list[trio_cat_1]->get_hap_num(); test_weight = fam_list[fam]->test_weight[i]; if(missing_trio_list[mis_trio_cat_1]->typed_par1==1 || missing_trio_list[mis_trio_cat_1]->typed_par2==1){ test_trio_weight1 = fam_list[fam]->test_trio_weight1[i]; test_trio_weight2 = fam_list[fam]->test_trio_weight2[i]; test_trio_weight3 = fam_list[fam]->test_trio_weight3[i]; for(int j=0;jderivative[j][k]; cholent[nb1+j][h_size+k] = test_trio_weight1*(missing_trio_list[mis_trio_cat_1]->deriv1[j][k])+test_trio_weight2*(missing_trio_list[mis_trio_cat_1]->deriv2[j][k])+test_trio_weight3*(missing_trio_list[mis_trio_cat_1]->deriv3[j][k]); } cholent[nb1+j][h_size+missing_trio_list[mis_trio_cat_1]->sites[j]]+=test_weight; cholent[nb1+j][h_size-1] = hap_num[j]; } } else{ for(int j=0;jderivative[j][k]; cholent[nb1+j][h_size+k] = test_weight*cholent[nb1+j][k]; } cholent[nb1+j][h_size-1] = hap_num[j]; } } nb1+=dim_1; } // i is typed } // complete i if(cholesky(covMatrix,dim,cholent,2*h_size-1,chol,cholaug,0) != 1){ printf("cholesky decomposition of the cov matrix failed for family %d. \n", fam_list[fam]->fam_id); flag_pd = 0; exit(1); } if(flag_pd==1){ double **s = self_colmultiply(cholaug,dim,1,h_size-1,1,h_size-1); double **t = self_colmultiply(cholaug,dim,h_size+1,2*h_size-1,h_size+1,2*h_size-1); double **u = self_colmultiply(cholaug,dim,1,h_size-1,h_size+1,2*h_size-1); double **v = self_colmultiply(cholaug,dim,h_size+1,2*h_size-1,h_size,h_size); for(int i=0;i *cluster){ col_1df = cluster->size()-1; int col_clus = col_1df; if((*cluster)[col_clus][0]==0) col_clus--; df = col_clus; if(df<=0){ fprintf(outfp, "The df is %d, please check!\n\n", df); return -1; } int h_size = hap_freq_table.size(); int col = geno_freq_table.size(); for(int i=0;i::iterator iter=hap_freq_table.begin(); while(iter!=hap_freq_table.end()){ freq[loc] = iter->second; loc++; iter++; } loc = 0; map, double>::iterator giter=geno_freq_table.begin(); while(giter!=geno_freq_table.end()){ genofreq[loc] = giter->second; loc++; giter++; } // create a matrix to cluster haplotypes double **mul_clus = new double*[h_size-1]; for(int i=0;ifirst; giter=geno_freq_table.begin(); while(giter!=geno_freq_table.end()){ pair conf = giter->first; loc = 0; iter = hap_freq_table.begin(); while(locfirst; if(conf.first == conf.second && conf.first == h) rst[loc][hap_pair_map[conf]] = 2*hap_freq_table[h]; else if(conf.first == conf.second && conf.first == h_last) rst[loc][hap_pair_map[conf]] = -2*hap_freq_table[h_last]; else if(conf.first < conf.second && conf.first == h && conf.second != h_last) rst[loc][hap_pair_map[conf]] = 2*hap_freq_table[conf.second]; else if(conf.first < conf.second && conf.second == h) rst[loc][hap_pair_map[conf]] = 2*hap_freq_table[conf.first]; else if(conf.first == h && conf.second == h_last) rst[loc][hap_pair_map[conf]] = 2*(hap_freq_table[h_last]-hap_freq_table[h]); else if(conf.first != h && conf.first != h_last && conf.second == h_last) rst[loc][hap_pair_map[conf]] = -2*hap_freq_table[conf.first]; loc++; iter++; } giter++; } double ***drst; if(oneparent==1){ drst = new double**[h_size-1]; for(int i=0;i, double>::iterator giter1; map, double>::iterator giter2; loc1 = 0; giter1=geno_freq_table.begin(); while(giter1!=geno_freq_table.end()){ int i=(giter1->first).first; int j=(giter1->first).second; value1 = 1.0/hap_freq_table[i]; value2 = 1.0/hap_freq_table[j]; value3 = 1.0/(hap_freq_table[i]+hap_freq_table[j]); loc2 = 0; giter2=geno_freq_table.begin(); while(giter2!=geno_freq_table.end()){ int k=(giter2->first).first; int l=(giter2->first).second; if(k==l){ // parent is homozygote if(i==k){ loc = hap_map[j]; if(loc!=h_size-1) drst[loc][loc1][loc2]+=value2; else{ for(int s=0;sget_num()>0 && geno_trio_list[i]->info==1 && missing_trio_list[geno_missing_trio_map[i]]->dim>0) geno_trio_list[i]->expected_hap_num(hap_freq_table,freq,geno_freq_table,hap_pair_map,count,missing_trio_list[geno_missing_trio_map[i]]->sites); } for(int i=0;iget_num()>0 && missing_trio_list[i]->dim>0){ missing_trio_list[i]->set_weight(hap_freq_table,geno_freq_table,hap_pair_map,count); missing_trio_list[i]->compute_centmatrix(hap_freq_table,freq,geno_freq_table); missing_trio_list[i]->compute_derivative(hap_freq_table,geno_freq_table,hap_pair_map,count,rst,hst); missing_trio_list[i]->compute_derivative_cluster(hap_freq_table,mul_clus,col_clus); missing_trio_list[i]->compute_derivative_1df(hap_freq_table,mul_1df,col_1df); missing_trio_list[i]->expected_hap_num(hap_freq_table,freq,geno_freq_table,hap_pair_map,count); if(missing_trio_list[i]->typed_par1==1 || missing_trio_list[i]->typed_par2==1){ missing_trio_list[i]->compute_deriv1to3(hap_freq_table,geno_freq_table,hap_pair_map,count,rst,drst); missing_trio_list[i]->compute_deriv1to3_cluster(hap_freq_table,mul_clus,col_clus); missing_trio_list[i]->compute_deriv1to3_1df(hap_freq_table,mul_1df,col_1df); } } } if(oneparent==1){ for(int i=0;i0){ map::iterator pi=par1_one.begin(); for(int i=0;isecond = i; int mis = pi->first; int dim_1 = missing_trio_list[mis]->dim; for(int j=0;jcond_par1_one(result_par1_one[i],hap_freq_table,geno_freq_table,hap_map); pi++; } } int par1_two_size = par1_two.size(); if(par1_two_size>0){ map::iterator pi=par1_two.begin(); for(int i=0;isecond = i; int mis = pi->first; int dim_1 = missing_trio_list[mis]->dim; for(int j=0;jcond_par1_two(result_par1_two[i],hap_freq_table,geno_freq_table,hap_pair_map); pi++; } } int par2_one_size = par2_one.size(); if(par2_one_size>0){ map::iterator pi=par2_one.begin(); for(int i=0;isecond = i; int mis = pi->first; int dim_1 = missing_trio_list[mis]->dim; for(int j=0;jcond_par2_one(result_par2_one[i],hap_freq_table,geno_freq_table,hap_map); pi++; } } int par2_two_size = par2_two.size(); if(par2_two_size>0){ map::iterator pi=par2_two.begin(); for(int i=0;isecond = i; int mis = pi->first; int dim_1 = missing_trio_list[mis]->dim; for(int j=0;jcond_par2_two(result_par2_two[i],hap_freq_table,geno_freq_table,hap_pair_map); pi++; } } int off_one_size = off_one.size(); if(off_one_size>0){ map::iterator pi=off_one.begin(); for(int i=0;isecond = i; int mis = pi->first; int dim_1 = missing_trio_list[mis]->dim; for(int j=0;jcond_off_one(result_off_one[i],hap_freq_table,geno_freq_table,hap_map); pi++; } } int off_two_size = off_two.size(); if(off_two_size>0){ map::iterator pi=off_two.begin(); for(int i=0;isecond = i; int mis = pi->first; int dim_1 = missing_trio_list[mis]->dim; for(int j=0;jcond_off_two(result_off_two[i],hap_freq_table,geno_freq_table,hap_pair_map); pi++; } } int parents_one_size = parents_one.size(); if(parents_one_size>0){ map::iterator pi=parents_one.begin(); for(int i=0;isecond = i; int mis = pi->first; int dim_1 = missing_trio_list[mis]->dim; for(int j=0;jcond_parents_one(result_parents_one[i],hap_freq_table,geno_freq_table,hap_map); pi++; } } int parents_two_size = parents_two.size(); if(parents_two_size>0){ map::iterator pi=parents_two.begin(); for(int i=0;isecond = i; int mis = pi->first; int dim_1 = missing_trio_list[mis]->dim; for(int j=0;jcond_parents_two(result_parents_two[i],hap_freq_table,geno_freq_table,hap_pair_map); pi++; } } for(int i=0;iN_typed>0){ int trio_cat; int mis_trio_cat; int dim=0; for(int i=0;iN;i++){ trio_cat = fam_list[fam]->trio_cat[i]; if(geno_trio_list[trio_cat]->info==1){ mis_trio_cat = fam_list[fam]->mis_trio_cat[i]; dim+=missing_trio_list[mis_trio_cat]->dim; } } if(dim<=0) continue; for(int i=0;i, double>::iterator iterk; map, double>::iterator iter1; map, double>::iterator iter2; double value; double test_weight; double test_trio_weight1, test_trio_weight2, test_trio_weight3; flag_pd = 1; for(int i=0;iN;i++){ trio_cat_1 = fam_list[fam]->trio_cat[i]; mis_trio_cat_1 = fam_list[fam]->mis_trio_cat[i]; dim_1 = missing_trio_list[mis_trio_cat_1]->dim; if(dim_1 > 0 && geno_trio_list[trio_cat_1]->info==1){ nb2 = nb1; ind_1 = fam_list[fam]->fam_member[i][0]; ind_2 = ind_1; iterk = (fam_list[fam]->kinship_coeff).find(pair(ind_1,ind_2)); if(iterk == (fam_list[fam]->kinship_coeff).end()){ printf("No IBD coefficient for individual %d from family %d. Please check...\n\n", ind_1,fam_list[fam]->fam_id); exit(1); } for(int j=0;jcentmatrix[j][k]; nb2+=dim_1; par_1 = 0; par_1+=missing_trio_list[mis_trio_cat_1]->typed_par1; par_1+=missing_trio_list[mis_trio_cat_1]->typed_par2; ind_1f = fam_list[fam]->fam_member[i][1]; ind_1m = fam_list[fam]->fam_member[i][2]; for(int j=i+1;jN;j++){ trio_cat_2 = fam_list[fam]->trio_cat[j]; mis_trio_cat_2 = fam_list[fam]->mis_trio_cat[j]; dim_2 = missing_trio_list[mis_trio_cat_2]->dim; if(dim_2 > 0 && geno_trio_list[trio_cat_2]->info==1){ ind_2 = fam_list[fam]->fam_member[j][0]; iterk = (fam_list[fam]->kinship_coeff).find(pair(ind_1,ind_2)); if(iterk == (fam_list[fam]->kinship_coeff).end()){ printf("No IBD coefficient between individual %d and individual %d from family %d. Please check...\n\n", ind_1,ind_2,fam_list[fam]->fam_id); exit(1); } if(iterk->second>1e-8){ // i is related to j par_2 = 0; par_2+=missing_trio_list[mis_trio_cat_2]->typed_par1; par_2+=missing_trio_list[mis_trio_cat_2]->typed_par2; ind_2f = fam_list[fam]->fam_member[j][1]; ind_2m = fam_list[fam]->fam_member[j][2]; kin = iterk->second; iter2 = (fam_list[fam]->IBD2_coeff).find(pair(ind_1,ind_2)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1,ind_2)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(par_1==0 && par_2==0){ if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;kmis_par_order[j]; int ind_2par1 = ind_2f; int ind_2par2 = ind_2m; if(order==1){ ind_2par1 = ind_2m; ind_2par2 = ind_2f; } int diff = 0; int index = relation(fam,ind_1,ind_2,ind_2par1,ind_2par2,kin,IBD_2,IBD_1,diff); if(index==4){ for(int k=0;kIBD2_coeff).find(pair(ind_1,ind_2par1)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1,ind_2par1)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;kIBD2_coeff).find(pair(ind_1,ind_2par2)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1,ind_2par2)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;k1e-8){ for(int k=0;k1e-8){ for(int k=0;kIBD1_coeff).find(pair(ind_1,ind_2par1)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;kIBD1_coeff).find(pair(ind_1,ind_2par2)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;kmis_par_order[i]; int ind_1par1 = ind_1f; int ind_1par2 = ind_1m; if(order==1){ ind_1par1 = ind_1m; ind_1par2 = ind_1f; } int diff = 0; int index = relation(fam,ind_2,ind_1,ind_1par1,ind_1par2,kin,IBD_2,IBD_1,diff); if(index==4){ for(int k=0;kIBD2_coeff).find(pair(ind_1par1,ind_2)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1par1,ind_2)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;kIBD2_coeff).find(pair(ind_1par2,ind_2)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1par2,ind_2)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;k1e-8){ for(int k=0;k1e-8){ for(int k=0;kIBD1_coeff).find(pair(ind_1par1,ind_2)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;kIBD1_coeff).find(pair(ind_1par2,ind_2)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;kmis_par_order[i]; int ind_1par1 = ind_1f; int ind_1par2 = ind_1m; if(order_1==1){ ind_1par1 = ind_1m; ind_1par2 = ind_1f; } int order_2 = fam_list[fam]->mis_par_order[j]; int ind_2par1 = ind_2f; int ind_2par2 = ind_2m; if(order_2==1){ ind_2par1 = ind_2m; ind_2par2 = ind_2f; } int diff = 0; int index = relation(fam,ind_1,ind_1par1,ind_1par2,ind_2,ind_2par1,ind_2par2,kin,IBD_2,IBD_1,diff); if(index==13){ for(int k=0;kIBD2_coeff).find(pair(ind_1par1,ind_2par1)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1par1,ind_2par1)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;kIBD2_coeff).find(pair(ind_1par1,ind_2par2)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1par1,ind_2par2)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;kIBD2_coeff).find(pair(ind_1par2,ind_2par1)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1par2,ind_2par1)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;kIBD2_coeff).find(pair(ind_1par2,ind_2par2)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1par2,ind_2par2)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;kIBD2_coeff).find(pair(ind_1par1,ind_2)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1par1,ind_2)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;kIBD2_coeff).find(pair(ind_1par2,ind_2)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1par2,ind_2)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;kIBD2_coeff).find(pair(ind_1,ind_2par1)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1,ind_2par1)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;kIBD2_coeff).find(pair(ind_1,ind_2par2)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1,ind_2par2)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;kfirst).first; int h2 = (giter->first).second; if(h1==h2){ for(int k=0;kIBD1_coeff).find(pair(ind_1par1,ind_2par1)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;kIBD1_coeff).find(pair(ind_1par1,ind_2par2)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;kfirst).first; int h2 = (giter->first).second; if(h1==h2){ for(int k=0;kIBD1_coeff).find(pair(ind_1par2,ind_2par1)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;kIBD1_coeff).find(pair(ind_1par2,ind_2par2)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;kfirst).first; int h2 = (giter->first).second; if(h1==h2){ for(int k=0;kIBD1_coeff).find(pair(ind_1par1,ind_2par1)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;kIBD1_coeff).find(pair(ind_1par2,ind_2par1)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;kfirst).first; int h2 = (giter->first).second; if(h1==h2){ for(int k=0;kIBD1_coeff).find(pair(ind_1par1,ind_2par2)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;kIBD1_coeff).find(pair(ind_1par2,ind_2par2)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k hap_num = geno_trio_list[trio_cat_1]->get_hap_num(); test_weight = fam_list[fam]->test_weight[i]; if(missing_trio_list[mis_trio_cat_1]->typed_par1==1 || missing_trio_list[mis_trio_cat_1]->typed_par2==1){ test_trio_weight1 = fam_list[fam]->test_trio_weight1[i]; test_trio_weight2 = fam_list[fam]->test_trio_weight2[i]; test_trio_weight3 = fam_list[fam]->test_trio_weight3[i]; for(int j=0;jderivative[j][k]; cholent[nb1+j][h_size-1] = hap_num[j]; } for(int j=0;jderiv1_clus[j][k])+test_trio_weight2*(missing_trio_list[mis_trio_cat_1]->deriv2_clus[j][k])+test_trio_weight3*(missing_trio_list[mis_trio_cat_1]->deriv3_clus[j][k]); cholent[nb1+j][h_size+k]+=test_weight*mul_clus[missing_trio_list[mis_trio_cat_1]->sites[j]][k]; } for(int k=0;kderiv1_1df[j][k])+test_trio_weight2*(missing_trio_list[mis_trio_cat_1]->deriv2_1df[j][k])+test_trio_weight3*(missing_trio_list[mis_trio_cat_1]->deriv3_1df[j][k]); cholent[nb1+j][h_size+col_clus+k]+=test_weight*mul_1df[missing_trio_list[mis_trio_cat_1]->sites[j]][k]; } } } else{ for(int j=0;jderivative[j][k]; cholent[nb1+j][h_size-1] = hap_num[j]; } for(int j=0;jderivative_clus[j][k]; for(int k=0;kderivative_1df[j][k]; } } nb1+=dim_1; } // i is typed } // complete i if(cholesky(covMatrix,dim,cholent,h_size+col_clus+col_1df,chol,cholaug,0) != 1){ printf("cholesky decomposition of the cov matrix failed for family %d. \n", fam_list[fam]->fam_id); flag_pd = 0; exit(1); } if(flag_pd==1){ double **s = self_colmultiply(cholaug,dim,1,h_size-1,1,h_size-1); double **t = self_colmultiply(cholaug,dim,h_size+1,h_size+col_clus,h_size+1,h_size+col_clus); double **u = self_colmultiply(cholaug,dim,1,h_size-1,h_size+1,h_size+col_clus); double **v = self_colmultiply(cholaug,dim,h_size+1,h_size+col_clus,h_size,h_size); for(int i=0;i, double>::iterator iter; iter = (fam_list[fam]->kinship_coeff).find(pair(ind_1,ind_2f)); if(iter == (fam_list[fam]->kinship_coeff).end()) { printf("No kinship coefficient between individual %d and individual %d from family %d. Please check...\n\n", ind_1,ind_2f,fam_list[fam]->fam_id); exit(1); } value_f = iter->second; iter = (fam_list[fam]->kinship_coeff).find(pair(ind_1,ind_2m)); if(iter == (fam_list[fam]->kinship_coeff).end()) { printf("No kinship coefficient between individual %d and individual %d from family %d. Please check...\n\n", ind_1,ind_2m,fam_list[fam]->fam_id); exit(1); } value_m = iter->second; if(value_m<1e-8) // i is unrelated to jm return 1; else if(value_f<1e-8) // i is unrelated to jf return 2; else if(value_o>value_f && value_o>value_m) // j is an ancestor of i return 3; else if(IBD_2<1e-8) // jf and jm are great...great grandparents of i { double gen = log(2.0/IBD_1)/log(2.0); diff = (gen-floor(gen)>=0.5) ? int(ceil(gen)) : int(floor(gen)); return 6; } return 0; } } int IQLBEstimator::relation(int fam, int ind_1, int ind_1f, int ind_1m, int ind_2, int ind_2f, int ind_2m, double value_o, double IBD_2, double IBD_1, int &diff) { // assume ind_1 != ind_2 if(ind_1==ind_2f) return 13; else if(ind_1==ind_2m) return 14; else if(ind_2==ind_1f) return 15; else if(ind_2==ind_1m) return 16; else if(ind_1f==ind_2f && ind_1m==ind_2m) // sib return 17; else if((ind_1f==ind_2f || ind_1m==ind_2m || ind_1f==ind_2m || ind_1m==ind_2f) && IBD_2<1e-8 && fabs(IBD_1-0.5)<1e-8) // half-sib return 18; else if(fabs(IBD_2-0.0625)<1e-8 && fabs(IBD_1-0.375)<1e-8) // double first cousin return 19; else if(fabs(IBD_2-0.125)<1e-8 && fabs(IBD_1-0.5)<1e-8) // half-sib plus first cousin return 20; else { double value_of, value_om, value_fo, value_mo, value_ff, value_fm, value_mf, value_mm; map, double>::iterator iter; iter = (fam_list[fam]->kinship_coeff).find(pair(ind_1,ind_2f)); if(iter == (fam_list[fam]->kinship_coeff).end()) { printf("No kinship coefficient between individual %d and individual %d from family %d. Please check...\n\n", ind_1,ind_2f,fam_list[fam]->fam_id); exit(1); } value_of = iter->second; iter = (fam_list[fam]->kinship_coeff).find(pair(ind_1,ind_2m)); if(iter == (fam_list[fam]->kinship_coeff).end()) { printf("No kinship coefficient between individual %d and individual %d from family %d. Please check...\n\n", ind_1,ind_2m,fam_list[fam]->fam_id); exit(1); } value_om = iter->second; iter = (fam_list[fam]->kinship_coeff).find(pair(ind_2,ind_1f)); if(iter == (fam_list[fam]->kinship_coeff).end()) { printf("No kinship coefficient between individual %d and individual %d from family %d. Please check...\n\n", ind_2,ind_1f,fam_list[fam]->fam_id); exit(1); } value_fo = iter->second; iter = (fam_list[fam]->kinship_coeff).find(pair(ind_2,ind_1m)); if(iter == (fam_list[fam]->kinship_coeff).end()) { printf("No kinship coefficient between individual %d and individual %d from family %d. Please check...\n\n", ind_2,ind_1m,fam_list[fam]->fam_id); exit(1); } value_mo = iter->second; iter = (fam_list[fam]->kinship_coeff).find(pair(ind_1f,ind_2f)); if(iter == (fam_list[fam]->kinship_coeff).end()) { printf("No kinship coefficient between individual %d and individual %d from family %d. Please check...\n\n", ind_1f,ind_2f,fam_list[fam]->fam_id); exit(1); } value_ff = iter->second; iter = (fam_list[fam]->kinship_coeff).find(pair(ind_1f,ind_2m)); if(iter == (fam_list[fam]->kinship_coeff).end()) { printf("No kinship coefficient between individual %d and individual %d from family %d. Please check...\n\n", ind_1f,ind_2m,fam_list[fam]->fam_id); exit(1); } value_fm = iter->second; iter = (fam_list[fam]->kinship_coeff).find(pair(ind_1m,ind_2f)); if(iter == (fam_list[fam]->kinship_coeff).end()) { printf("No kinship coefficient between individual %d and individual %d from family %d. Please check...\n\n", ind_1m,ind_2f,fam_list[fam]->fam_id); exit(1); } value_mf = iter->second; iter = (fam_list[fam]->kinship_coeff).find(pair(ind_1m,ind_2m)); if(iter == (fam_list[fam]->kinship_coeff).end()) { printf("No kinship coefficient between individual %d and individual %d from family %d. Please check...\n\n", ind_1m,ind_2m,fam_list[fam]->fam_id); exit(1); } value_mm = iter->second; if(value_om<1e-8 && value_mo<1e-8) // if is related to jf return 1; else if(value_of<1e-8 && value_mo<1e-8) // if is related to jm return 2; else if(value_om<1e-8 && value_fo<1e-8) // im is related to jf return 3; else if(value_of<1e-8 && value_fo<1e-8) // im is related to jm return 4; else if(value_mo<1e-8 && value_fo>value_ff && value_fo>value_fm) // j is an ancestor of if return 5; else if(value_fo<1e-8 && value_mo>value_mf && value_mo>value_mm) // j is an ancestor of im return 6; else if(value_om<1e-8 && value_of>value_ff && value_of>value_mf) // i is an ancestor of jf return 7; else if(value_of<1e-8 && value_om>value_fm && value_om>value_mm) // i is an ancestor of jm return 8; else if(value_mo<1e-8 && IBD_2<1e-8) // jf and jm are great...great grandparents of i, if { double gen = log(2.0/IBD_1)/log(2.0)-1; diff = (gen-floor(gen)>=0.5) ? int(ceil(gen)) : int(floor(gen)); return 9; } else if(value_fo<1e-8 && IBD_2<1e-8) // jf and jm are great...great grandparents of i, im { double gen = log(2.0/IBD_1)/log(2.0)-1; diff = (gen-floor(gen)>=0.5) ? int(ceil(gen)) : int(floor(gen)); return 10; } else if(value_om<1e-8 && IBD_2<1e-8) // if and im are great...great grandparents of j, jf { double gen = log(2.0/IBD_1)/log(2.0)-1; diff = (gen-floor(gen)>=0.5) ? int(ceil(gen)) : int(floor(gen)); return 11; } else if(value_of<1e-8 && IBD_2<1e-8) // if and im are great...great grandparents of j, jm { double gen = log(2.0/IBD_1)/log(2.0)-1; diff = (gen-floor(gen)>=0.5) ? int(ceil(gen)) : int(floor(gen)); return 12; } return 0; } } void IQLBEstimator::estimate_subgroup(){ full_hap_freq_table = hap_freq_table; // affected // unaffected // unknown aff=0; unaff=0; unknown=0; int trio_cat_1; int mis_trio_cat_1; int dim_1; for(int fam=0;famN_typed>0){ for(int i=0;iN;i++){ trio_cat_1 = fam_list[fam]->trio_cat[i]; mis_trio_cat_1 = fam_list[fam]->mis_trio_cat[i]; dim_1 = missing_trio_list[mis_trio_cat_1]->dim; if(dim_1 > 0 && geno_trio_list[trio_cat_1]->info==1){ if(fam_list[fam]->fam_member[i][4]==0) unknown++; else if(fam_list[fam]->fam_member[i][4]==1) unaff++; else if(fam_list[fam]->fam_member[i][4]==2) aff++; } } } } if(aff>=15){ reset_freq_table(full_hap_freq_table); int count=1; int flag; while(1){ map of = hap_freq_table; NR_Iteration_subgroup(2); map nf = hap_freq_table; double diff = 0; flag = 1; map::iterator iter=hap_freq_table.begin(); while(iter!=hap_freq_table.end()){ diff += (nf[iter->first]-of[iter->first])*(nf[iter->first]-of[iter->first]); if(nf[iter->first]<0.0000001){ flag = 0; break; } iter++; } if(diff<1e-5 || count>5 || flag==0) break; count++; } aff_hap_freq_table = hap_freq_table; } if(unaff>=15){ reset_freq_table(full_hap_freq_table); int count=1; int flag; while(1){ map of = hap_freq_table; NR_Iteration_subgroup(1); map nf = hap_freq_table; double diff = 0; flag = 1; map::iterator iter=hap_freq_table.begin(); while(iter!=hap_freq_table.end()){ diff += (nf[iter->first]-of[iter->first])*(nf[iter->first]-of[iter->first]); if(nf[iter->first]<0.0000001){ flag = 0; break; } iter++; } if(diff<1e-5 || count>5 || flag==0) break; count++; } unaff_hap_freq_table = hap_freq_table; } if(unknown>=15){ reset_freq_table(full_hap_freq_table); int count=1; int flag; while(1){ map of = hap_freq_table; NR_Iteration_subgroup(0); map nf = hap_freq_table; double diff = 0; flag = 1; map::iterator iter=hap_freq_table.begin(); while(iter!=hap_freq_table.end()){ diff += (nf[iter->first]-of[iter->first])*(nf[iter->first]-of[iter->first]); if(nf[iter->first]<0.0000001){ flag = 0; break; } iter++; } if(diff<1e-5 || count>5 || flag==0) break; count++; } unknown_hap_freq_table = hap_freq_table; } reset_freq_table(full_hap_freq_table); } void IQLBEstimator::NR_Iteration_subgroup(int group){ int h_size = hap_freq_table.size(); int col = geno_freq_table.size(); for(int i=0;i::iterator iter=hap_freq_table.begin(); while(iter!=hap_freq_table.end()){ freq[loc] = iter->second; loc++; iter++; } loc = 0; map, double>::iterator giter=geno_freq_table.begin(); while(giter!=geno_freq_table.end()){ genofreq[loc] = giter->second; loc++; giter++; } int h_last=0; iter = hap_freq_table.end(); iter--; h_last=iter->first; giter=geno_freq_table.begin(); while(giter!=geno_freq_table.end()){ pair conf = giter->first; loc = 0; iter = hap_freq_table.begin(); while(locfirst; if(conf.first == conf.second && conf.first == h) rst[loc][hap_pair_map[conf]] = 2*hap_freq_table[h]; else if(conf.first == conf.second && conf.first == h_last) rst[loc][hap_pair_map[conf]] = -2*hap_freq_table[h_last]; else if(conf.first < conf.second && conf.first == h && conf.second != h_last) rst[loc][hap_pair_map[conf]] = 2*hap_freq_table[conf.second]; else if(conf.first < conf.second && conf.second == h) rst[loc][hap_pair_map[conf]] = 2*hap_freq_table[conf.first]; else if(conf.first == h && conf.second == h_last) rst[loc][hap_pair_map[conf]] = 2*(hap_freq_table[h_last]-hap_freq_table[h]); else if(conf.first != h && conf.first != h_last && conf.second == h_last) rst[loc][hap_pair_map[conf]] = -2*hap_freq_table[conf.first]; loc++; iter++; } giter++; } for(int i=0;iget_num()>0 && geno_trio_list[i]->info==1 && missing_trio_list[geno_missing_trio_map[i]]->dim>0) geno_trio_list[i]->expected_hap_num(hap_freq_table,freq,geno_freq_table,hap_pair_map,count,missing_trio_list[geno_missing_trio_map[i]]->sites); } for(int i=0;iget_num()>0 && missing_trio_list[i]->dim>0){ missing_trio_list[i]->set_weight(hap_freq_table,geno_freq_table,hap_pair_map,count); missing_trio_list[i]->compute_centmatrix(hap_freq_table,freq,geno_freq_table); missing_trio_list[i]->compute_derivative(hap_freq_table,geno_freq_table,hap_pair_map,count,rst,hst); missing_trio_list[i]->expected_hap_num(hap_freq_table,freq,geno_freq_table,hap_pair_map,count); } } int par1_one_size = par1_one.size(); if(par1_one_size>0){ map::iterator pi=par1_one.begin(); for(int i=0;isecond = i; int mis = pi->first; int dim_1 = missing_trio_list[mis]->dim; for(int j=0;jcond_par1_one(result_par1_one[i],hap_freq_table,geno_freq_table,hap_map); pi++; } } int par1_two_size = par1_two.size(); if(par1_two_size>0){ map::iterator pi=par1_two.begin(); for(int i=0;isecond = i; int mis = pi->first; int dim_1 = missing_trio_list[mis]->dim; for(int j=0;jcond_par1_two(result_par1_two[i],hap_freq_table,geno_freq_table,hap_pair_map); pi++; } } int par2_one_size = par2_one.size(); if(par2_one_size>0){ map::iterator pi=par2_one.begin(); for(int i=0;isecond = i; int mis = pi->first; int dim_1 = missing_trio_list[mis]->dim; for(int j=0;jcond_par2_one(result_par2_one[i],hap_freq_table,geno_freq_table,hap_map); pi++; } } int par2_two_size = par2_two.size(); if(par2_two_size>0){ map::iterator pi=par2_two.begin(); for(int i=0;isecond = i; int mis = pi->first; int dim_1 = missing_trio_list[mis]->dim; for(int j=0;jcond_par2_two(result_par2_two[i],hap_freq_table,geno_freq_table,hap_pair_map); pi++; } } int off_one_size = off_one.size(); if(off_one_size>0){ map::iterator pi=off_one.begin(); for(int i=0;isecond = i; int mis = pi->first; int dim_1 = missing_trio_list[mis]->dim; for(int j=0;jcond_off_one(result_off_one[i],hap_freq_table,geno_freq_table,hap_map); pi++; } } int off_two_size = off_two.size(); if(off_two_size>0){ map::iterator pi=off_two.begin(); for(int i=0;isecond = i; int mis = pi->first; int dim_1 = missing_trio_list[mis]->dim; for(int j=0;jcond_off_two(result_off_two[i],hap_freq_table,geno_freq_table,hap_pair_map); pi++; } } int parents_one_size = parents_one.size(); if(parents_one_size>0){ map::iterator pi=parents_one.begin(); for(int i=0;isecond = i; int mis = pi->first; int dim_1 = missing_trio_list[mis]->dim; for(int j=0;jcond_parents_one(result_parents_one[i],hap_freq_table,geno_freq_table,hap_map); pi++; } } int parents_two_size = parents_two.size(); if(parents_two_size>0){ map::iterator pi=parents_two.begin(); for(int i=0;isecond = i; int mis = pi->first; int dim_1 = missing_trio_list[mis]->dim; for(int j=0;jcond_parents_two(result_parents_two[i],hap_freq_table,geno_freq_table,hap_pair_map); pi++; } } for(int i=0;iN_typed>0){ int trio_cat; int mis_trio_cat; int dim=0; for(int i=0;iN;i++){ trio_cat = fam_list[fam]->trio_cat[i]; if(geno_trio_list[trio_cat]->info==1 && fam_list[fam]->fam_member[i][4]==group){ mis_trio_cat = fam_list[fam]->mis_trio_cat[i]; dim+=missing_trio_list[mis_trio_cat]->dim; } } if(dim<=0) continue; for(int i=0;i, double>::iterator iterk; map, double>::iterator iter1; map, double>::iterator iter2; double value; flag_pd = 1; for(int i=0;iN;i++){ trio_cat_1 = fam_list[fam]->trio_cat[i]; mis_trio_cat_1 = fam_list[fam]->mis_trio_cat[i]; dim_1 = missing_trio_list[mis_trio_cat_1]->dim; if(dim_1 > 0 && geno_trio_list[trio_cat_1]->info==1 && fam_list[fam]->fam_member[i][4]==group){ nb2 = nb1; ind_1 = fam_list[fam]->fam_member[i][0]; ind_2 = ind_1; iterk = (fam_list[fam]->kinship_coeff).find(pair(ind_1,ind_2)); if(iterk == (fam_list[fam]->kinship_coeff).end()){ printf("No IBD coefficient for individual %d from family %d. Please check...\n\n", ind_1,fam_list[fam]->fam_id); exit(1); } for(int j=0;jcentmatrix[j][k]; nb2+=dim_1; par_1 = 0; par_1+=missing_trio_list[mis_trio_cat_1]->typed_par1; par_1+=missing_trio_list[mis_trio_cat_1]->typed_par2; ind_1f = fam_list[fam]->fam_member[i][1]; ind_1m = fam_list[fam]->fam_member[i][2]; for(int j=i+1;jN;j++){ trio_cat_2 = fam_list[fam]->trio_cat[j]; mis_trio_cat_2 = fam_list[fam]->mis_trio_cat[j]; dim_2 = missing_trio_list[mis_trio_cat_2]->dim; if(dim_2 > 0 && geno_trio_list[trio_cat_2]->info==1 && fam_list[fam]->fam_member[j][4]==group){ ind_2 = fam_list[fam]->fam_member[j][0]; iterk = (fam_list[fam]->kinship_coeff).find(pair(ind_1,ind_2)); if(iterk == (fam_list[fam]->kinship_coeff).end()){ printf("No IBD coefficient between individual %d and individual %d from family %d. Please check...\n\n", ind_1,ind_2,fam_list[fam]->fam_id); exit(1); } if(iterk->second>1e-8){ // i is related to j par_2 = 0; par_2+=missing_trio_list[mis_trio_cat_2]->typed_par1; par_2+=missing_trio_list[mis_trio_cat_2]->typed_par2; ind_2f = fam_list[fam]->fam_member[j][1]; ind_2m = fam_list[fam]->fam_member[j][2]; kin = iterk->second; iter2 = (fam_list[fam]->IBD2_coeff).find(pair(ind_1,ind_2)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1,ind_2)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(par_1==0 && par_2==0){ if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;kmis_par_order[j]; int ind_2par1 = ind_2f; int ind_2par2 = ind_2m; if(order==1){ ind_2par1 = ind_2m; ind_2par2 = ind_2f; } int diff = 0; int index = relation(fam,ind_1,ind_2,ind_2par1,ind_2par2,kin,IBD_2,IBD_1,diff); if(index==4){ for(int k=0;kIBD2_coeff).find(pair(ind_1,ind_2par1)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1,ind_2par1)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;kIBD2_coeff).find(pair(ind_1,ind_2par2)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1,ind_2par2)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;k1e-8){ for(int k=0;k1e-8){ for(int k=0;kIBD1_coeff).find(pair(ind_1,ind_2par1)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;kIBD1_coeff).find(pair(ind_1,ind_2par2)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;kmis_par_order[i]; int ind_1par1 = ind_1f; int ind_1par2 = ind_1m; if(order==1){ ind_1par1 = ind_1m; ind_1par2 = ind_1f; } int diff = 0; int index = relation(fam,ind_2,ind_1,ind_1par1,ind_1par2,kin,IBD_2,IBD_1,diff); if(index==4){ for(int k=0;kIBD2_coeff).find(pair(ind_1par1,ind_2)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1par1,ind_2)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;kIBD2_coeff).find(pair(ind_1par2,ind_2)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1par2,ind_2)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;k1e-8){ for(int k=0;k1e-8){ for(int k=0;kIBD1_coeff).find(pair(ind_1par1,ind_2)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;kIBD1_coeff).find(pair(ind_1par2,ind_2)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;kmis_par_order[i]; int ind_1par1 = ind_1f; int ind_1par2 = ind_1m; if(order_1==1){ ind_1par1 = ind_1m; ind_1par2 = ind_1f; } int order_2 = fam_list[fam]->mis_par_order[j]; int ind_2par1 = ind_2f; int ind_2par2 = ind_2m; if(order_2==1){ ind_2par1 = ind_2m; ind_2par2 = ind_2f; } int diff = 0; int index = relation(fam,ind_1,ind_1par1,ind_1par2,ind_2,ind_2par1,ind_2par2,kin,IBD_2,IBD_1,diff); if(index==13){ for(int k=0;kIBD2_coeff).find(pair(ind_1par1,ind_2par1)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1par1,ind_2par1)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;kIBD2_coeff).find(pair(ind_1par1,ind_2par2)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1par1,ind_2par2)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;kIBD2_coeff).find(pair(ind_1par2,ind_2par1)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1par2,ind_2par1)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;kIBD2_coeff).find(pair(ind_1par2,ind_2par2)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1par2,ind_2par2)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;kIBD2_coeff).find(pair(ind_1par1,ind_2)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1par1,ind_2)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;kIBD2_coeff).find(pair(ind_1par2,ind_2)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1par2,ind_2)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;kIBD2_coeff).find(pair(ind_1,ind_2par1)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1,ind_2par1)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;kIBD2_coeff).find(pair(ind_1,ind_2par2)); iter1 = (fam_list[fam]->IBD1_coeff).find(pair(ind_1,ind_2par2)); IBD_2 = iter2->second; IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k1e-8){ for(int k=0;kfirst).first; int h2 = (giter->first).second; if(h1==h2){ for(int k=0;kIBD1_coeff).find(pair(ind_1par1,ind_2par1)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;kIBD1_coeff).find(pair(ind_1par1,ind_2par2)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;kfirst).first; int h2 = (giter->first).second; if(h1==h2){ for(int k=0;kIBD1_coeff).find(pair(ind_1par2,ind_2par1)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;kIBD1_coeff).find(pair(ind_1par2,ind_2par2)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;kfirst).first; int h2 = (giter->first).second; if(h1==h2){ for(int k=0;kIBD1_coeff).find(pair(ind_1par1,ind_2par1)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;kIBD1_coeff).find(pair(ind_1par2,ind_2par1)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;kfirst).first; int h2 = (giter->first).second; if(h1==h2){ for(int k=0;kIBD1_coeff).find(pair(ind_1par1,ind_2par2)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;kIBD1_coeff).find(pair(ind_1par2,ind_2par2)); IBD_1 = iter1->second; if(IBD_1>1e-8){ for(int k=0;k hap_num = geno_trio_list[trio_cat_1]->get_hap_num(); for(int j=0;jderivative[j][k]; cholent[nb1+j][h_size-1] = hap_num[j]; } nb1+=dim_1; } // i is typed } // complete i if(cholesky(covMatrix,dim,cholent,h_size,chol,cholaug,0) != 1){ printf("cholesky decomposition of the cov matrix failed for family %d. \n", fam_list[fam]->fam_id); flag_pd = 0; exit(1); } if(flag_pd==1){ double **s = self_colmultiply(cholaug,dim,1,h_size-1,1,h_size-1); double **t = self_colmultiply(cholaug,dim,1,h_size-1,h_size,h_size); for(int i=0;isecond+=s[i][0]; n+=iter->second; iter++; } iter->second = 1-n; clear_matrix(s,h_size-1,1); } else{ double var; double &determ = var; double **v = inverse(denominator,h_size-1,determ); double **w = new double*[h_size-1]; for(int i=0;isecond+=s[i][0]; n+=iter->second; iter++; } iter->second = 1-n; clear_matrix(s,h_size-1,1); clear_matrix(v,h_size-1,h_size-1); clear_matrix(w,h_size-1,1); } clear_matrix(cholentmain,h_size-1,h_size); clear_matrix(cholmain,h_size-1,h_size-1); clear_matrix(cholaugmain,h_size-1,h_size); update_geno_freq(); } int compare_equ(int *a, int *b, int loci_num) { for(int k=0;k &a0v, vector &a1v){ fprintf(outfp, "Haplotype frequency estimates by phenotype category (IQL_b estimator)\n"); fprintf(outfp, "Haplotype\tOverall\t"); if(aff>=15) fprintf(outfp, "\tAffected"); if(unaff>=15) fprintf(outfp, "\tUnaffected"); if(unknown>=15) fprintf(outfp, "\tUnknown"); fprintf(outfp, "\n"); double sum=0, sum_aff=0, sum_unaff=0, sum_unknown=0; int ind=0, ind_aff=0, ind_unaff=0, ind_unknown=0; map::iterator iter = hap_freq_table.begin(); while(iter!=hap_freq_table.end()){ if(full_hap_freq_table[iter->first]>=0) sum+=full_hap_freq_table[iter->first]; else ind=1; if(aff>=15){ if(aff_hap_freq_table[iter->first]>=0) sum_aff+=aff_hap_freq_table[iter->first]; else ind_aff=1; } if(unaff>=15){ if(unaff_hap_freq_table[iter->first]>=0) sum_unaff+=unaff_hap_freq_table[iter->first]; else ind_unaff=1; } if(unknown>=15){ if(unknown_hap_freq_table[iter->first]>=0) sum_unknown+=unknown_hap_freq_table[iter->first]; else ind_unknown=1; } iter++; } iter = hap_freq_table.begin(); while(iter!=hap_freq_table.end()){ translate_hap(outfp,iter->first,loci_num,a0v,a1v); if(loci_num==2) fprintf(outfp, "\t"); if(full_hap_freq_table[iter->first]>=0){ if(ind==0) fprintf(outfp, "\t%.4f\t", full_hap_freq_table[iter->first]); else fprintf(outfp, "\t%.4f\t", full_hap_freq_table[iter->first]/sum); } else fprintf(outfp, "\t0.0000\t"); if(aff>=15){ if(aff_hap_freq_table[iter->first]>=0){ if(ind_aff==0) fprintf(outfp, "\t%.4f\t", aff_hap_freq_table[iter->first]); else fprintf(outfp, "\t%.4f\t", aff_hap_freq_table[iter->first]/sum_aff); } else fprintf(outfp, "\t0.0000\t"); } if(unaff>=15){ if(unaff_hap_freq_table[iter->first]>=0){ if(ind_unaff==0) fprintf(outfp, "\t%.4f\t", unaff_hap_freq_table[iter->first]); else fprintf(outfp, "\t%.4f\t", unaff_hap_freq_table[iter->first]/sum_unaff); } else fprintf(outfp, "\t0.0000\t"); } if(unknown>=15){ if(unknown_hap_freq_table[iter->first]>=0){ if(ind_unknown==0) fprintf(outfp, "\t%.4f", unknown_hap_freq_table[iter->first]); else fprintf(outfp, "\t%.4f", unknown_hap_freq_table[iter->first]/sum_unknown); } else fprintf(outfp, "\t0.0000"); } fprintf(outfp, "\n"); iter++; } fprintf(outfp, "\n"); } void IQLBEstimator::NR_Iteration_EW(){ int h_size = hap_freq_table.size(); int col = geno_freq_table.size(); for(int i=0;i::iterator iter=hap_freq_table.begin(); while(iter!=hap_freq_table.end()){ freq[loc] = iter->second; loc++; iter++; } loc = 0; map, double>::iterator giter=geno_freq_table.begin(); while(giter!=geno_freq_table.end()){ genofreq[loc] = giter->second; loc++; giter++; } int h_last=0; iter = hap_freq_table.end(); iter--; h_last=iter->first; giter=geno_freq_table.begin(); while(giter!=geno_freq_table.end()){ pair conf = giter->first; loc = 0; iter = hap_freq_table.begin(); while(locfirst; if(conf.first == conf.second && conf.first == h) rst[loc][hap_pair_map[conf]] = 2*hap_freq_table[h]; else if(conf.first == conf.second && conf.first == h_last) rst[loc][hap_pair_map[conf]] = -2*hap_freq_table[h_last]; else if(conf.first < conf.second && conf.first == h && conf.second != h_last) rst[loc][hap_pair_map[conf]] = 2*hap_freq_table[conf.second]; else if(conf.first < conf.second && conf.second == h) rst[loc][hap_pair_map[conf]] = 2*hap_freq_table[conf.first]; else if(conf.first == h && conf.second == h_last) rst[loc][hap_pair_map[conf]] = 2*(hap_freq_table[h_last]-hap_freq_table[h]); else if(conf.first != h && conf.first != h_last && conf.second == h_last) rst[loc][hap_pair_map[conf]] = -2*hap_freq_table[conf.first]; loc++; iter++; } giter++; } for(int i=0;iget_num()>0 && geno_trio_list[i]->info==1 && missing_trio_list[geno_missing_trio_map[i]]->dim>0) geno_trio_list[i]->expected_hap_num(hap_freq_table,freq,geno_freq_table,hap_pair_map,count,missing_trio_list[geno_missing_trio_map[i]]->sites); } for(int i=0;iget_num()>0 && missing_trio_list[i]->dim>0){ missing_trio_list[i]->set_weight(hap_freq_table,geno_freq_table,hap_pair_map,count); missing_trio_list[i]->compute_centmatrix(hap_freq_table,freq,geno_freq_table); missing_trio_list[i]->compute_derivative(hap_freq_table,geno_freq_table,hap_pair_map,count,rst,hst); missing_trio_list[i]->expected_hap_num(hap_freq_table,freq,geno_freq_table,hap_pair_map,count); } } for(int i=0;iN_typed>0){ int trio_cat; int mis_trio_cat; int dim=0; for(int i=0;iN;i++){ trio_cat = fam_list[fam]->trio_cat[i]; if(geno_trio_list[trio_cat]->info==1){ mis_trio_cat = fam_list[fam]->mis_trio_cat[i]; dim+=missing_trio_list[mis_trio_cat]->dim; } } if(dim<=0) continue; for(int i=0;i, double>::iterator iterk; flag_pd = 1; for(int i=0;iN;i++){ trio_cat_1 = fam_list[fam]->trio_cat[i]; mis_trio_cat_1 = fam_list[fam]->mis_trio_cat[i]; dim_1 = missing_trio_list[mis_trio_cat_1]->dim; if(dim_1 > 0 && geno_trio_list[trio_cat_1]->info==1){ ind_1 = fam_list[fam]->fam_member[i][0]; iterk = (fam_list[fam]->kinship_coeff).find(pair(ind_1,ind_1)); if(iterk == (fam_list[fam]->kinship_coeff).end()){ printf("No IBD coefficient for individual %d from family %d. Please check...\n\n", ind_1,fam_list[fam]->fam_id); exit(1); } for(int j=0;jcentmatrix[j][k]; const vector hap_num = geno_trio_list[trio_cat_1]->get_hap_num(); for(int j=0;jderivative[j][k]; cholent[nb1+j][h_size-1] = hap_num[j]; } nb1+=dim_1; } // i is typed } // complete i if(cholesky(covMatrix,dim,cholent,h_size,chol,cholaug,0) != 1){ printf("cholesky decomposition of the cov matrix failed for family %d. \n", fam_list[fam]->fam_id); flag_pd = 0; exit(1); } if(flag_pd==1){ double **s = self_colmultiply(cholaug,dim,1,h_size-1,1,h_size-1); double **t = self_colmultiply(cholaug,dim,1,h_size-1,h_size,h_size); for(int i=0;isecond+=s[i][0]; n+=iter->second; iter++; } iter->second = 1-n; clear_matrix(s,h_size-1,1); } else{ double var; double &determ = var; double **v = inverse(denominator,h_size-1,determ); double **w = new double*[h_size-1]; for(int i=0;isecond+=s[i][0]; n+=iter->second; iter++; } iter->second = 1-n; clear_matrix(s,h_size-1,1); clear_matrix(v,h_size-1,h_size-1); clear_matrix(w,h_size-1,1); } clear_matrix(cholentmain,h_size-1,h_size); clear_matrix(cholmain,h_size-1,h_size-1); clear_matrix(cholaugmain,h_size-1,h_size); update_geno_freq(); } void IQLBEstimator::NR_Iteration_subgroup_EW(int group){ int h_size = hap_freq_table.size(); int col = geno_freq_table.size(); for(int i=0;i::iterator iter=hap_freq_table.begin(); while(iter!=hap_freq_table.end()){ freq[loc] = iter->second; loc++; iter++; } loc = 0; map, double>::iterator giter=geno_freq_table.begin(); while(giter!=geno_freq_table.end()){ genofreq[loc] = giter->second; loc++; giter++; } int h_last=0; iter = hap_freq_table.end(); iter--; h_last=iter->first; giter=geno_freq_table.begin(); while(giter!=geno_freq_table.end()){ pair conf = giter->first; loc = 0; iter = hap_freq_table.begin(); while(locfirst; if(conf.first == conf.second && conf.first == h) rst[loc][hap_pair_map[conf]] = 2*hap_freq_table[h]; else if(conf.first == conf.second && conf.first == h_last) rst[loc][hap_pair_map[conf]] = -2*hap_freq_table[h_last]; else if(conf.first < conf.second && conf.first == h && conf.second != h_last) rst[loc][hap_pair_map[conf]] = 2*hap_freq_table[conf.second]; else if(conf.first < conf.second && conf.second == h) rst[loc][hap_pair_map[conf]] = 2*hap_freq_table[conf.first]; else if(conf.first == h && conf.second == h_last) rst[loc][hap_pair_map[conf]] = 2*(hap_freq_table[h_last]-hap_freq_table[h]); else if(conf.first != h && conf.first != h_last && conf.second == h_last) rst[loc][hap_pair_map[conf]] = -2*hap_freq_table[conf.first]; loc++; iter++; } giter++; } for(int i=0;iget_num()>0 && geno_trio_list[i]->info==1 && missing_trio_list[geno_missing_trio_map[i]]->dim>0) geno_trio_list[i]->expected_hap_num(hap_freq_table,freq,geno_freq_table,hap_pair_map,count,missing_trio_list[geno_missing_trio_map[i]]->sites); } for(int i=0;iget_num()>0 && missing_trio_list[i]->dim>0){ missing_trio_list[i]->set_weight(hap_freq_table,geno_freq_table,hap_pair_map,count); missing_trio_list[i]->compute_centmatrix(hap_freq_table,freq,geno_freq_table); missing_trio_list[i]->compute_derivative(hap_freq_table,geno_freq_table,hap_pair_map,count,rst,hst); missing_trio_list[i]->expected_hap_num(hap_freq_table,freq,geno_freq_table,hap_pair_map,count); } } for(int i=0;iN_typed>0){ int trio_cat; int mis_trio_cat; int dim=0; for(int i=0;iN;i++){ trio_cat = fam_list[fam]->trio_cat[i]; if(geno_trio_list[trio_cat]->info==1 && fam_list[fam]->fam_member[i][4]==group){ mis_trio_cat = fam_list[fam]->mis_trio_cat[i]; dim+=missing_trio_list[mis_trio_cat]->dim; } } if(dim<=0) continue; for(int i=0;i, double>::iterator iterk; flag_pd = 1; for(int i=0;iN;i++){ trio_cat_1 = fam_list[fam]->trio_cat[i]; mis_trio_cat_1 = fam_list[fam]->mis_trio_cat[i]; dim_1 = missing_trio_list[mis_trio_cat_1]->dim; if(dim_1 > 0 && geno_trio_list[trio_cat_1]->info==1 && fam_list[fam]->fam_member[i][4]==group){ ind_1 = fam_list[fam]->fam_member[i][0]; iterk = (fam_list[fam]->kinship_coeff).find(pair(ind_1,ind_1)); if(iterk == (fam_list[fam]->kinship_coeff).end()){ printf("No IBD coefficient for individual %d from family %d. Please check...\n\n", ind_1,fam_list[fam]->fam_id); exit(1); } for(int j=0;jcentmatrix[j][k]; const vector hap_num = geno_trio_list[trio_cat_1]->get_hap_num(); for(int j=0;jderivative[j][k]; cholent[nb1+j][h_size-1] = hap_num[j]; } nb1+=dim_1; } // i is typed } // complete i if(cholesky(covMatrix,dim,cholent,h_size,chol,cholaug,0) != 1){ printf("cholesky decomposition of the cov matrix failed for family %d. \n", fam_list[fam]->fam_id); flag_pd = 0; exit(1); } if(flag_pd==1){ double **s = self_colmultiply(cholaug,dim,1,h_size-1,1,h_size-1); double **t = self_colmultiply(cholaug,dim,1,h_size-1,h_size,h_size); for(int i=0;isecond+=s[i][0]; n+=iter->second; iter++; } iter->second = 1-n; clear_matrix(s,h_size-1,1); } else{ double var; double &determ = var; double **v = inverse(denominator,h_size-1,determ); double **w = new double*[h_size-1]; for(int i=0;isecond+=s[i][0]; n+=iter->second; iter++; } iter->second = 1-n; clear_matrix(s,h_size-1,1); clear_matrix(v,h_size-1,h_size-1); clear_matrix(w,h_size-1,1); } clear_matrix(cholentmain,h_size-1,h_size); clear_matrix(cholmain,h_size-1,h_size-1); clear_matrix(cholaugmain,h_size-1,h_size); update_geno_freq(); } void IQLBEstimator::estimate_subgroup_EW(){ reset_freq_table(full_hap_freq_table); int count=1; int flag; while(1){ map of = hap_freq_table; NR_Iteration_EW(); map nf = hap_freq_table; double diff = 0; flag = 1; map::iterator iter=hap_freq_table.begin(); while(iter!=hap_freq_table.end()){ diff += (nf[iter->first]-of[iter->first])*(nf[iter->first]-of[iter->first]); if(nf[iter->first]<0.0000001){ flag = 0; break; } iter++; } if(diff<1e-5 || count>5 || flag==0) break; count++; } full_hap_freq_table_EW = hap_freq_table; if(aff>=15){ reset_freq_table(full_hap_freq_table); int count=1; int flag; while(1){ map of = hap_freq_table; NR_Iteration_subgroup_EW(2); map nf = hap_freq_table; double diff = 0; flag = 1; map::iterator iter=hap_freq_table.begin(); while(iter!=hap_freq_table.end()){ diff += (nf[iter->first]-of[iter->first])*(nf[iter->first]-of[iter->first]); if(nf[iter->first]<0.0000001){ flag = 0; break; } iter++; } if(diff<1e-5 || count>5 || flag==0) break; count++; } aff_hap_freq_table_EW = hap_freq_table; } if(unaff>=15){ reset_freq_table(full_hap_freq_table); int count=1; int flag; while(1){ map of = hap_freq_table; NR_Iteration_subgroup_EW(1); map nf = hap_freq_table; double diff = 0; flag = 1; map::iterator iter=hap_freq_table.begin(); while(iter!=hap_freq_table.end()){ diff += (nf[iter->first]-of[iter->first])*(nf[iter->first]-of[iter->first]); if(nf[iter->first]<0.0000001){ flag = 0; break; } iter++; } if(diff<1e-5 || count>5 || flag==0) break; count++; } unaff_hap_freq_table_EW = hap_freq_table; } if(unknown>=15){ reset_freq_table(full_hap_freq_table); int count=1; int flag; while(1){ map of = hap_freq_table; NR_Iteration_subgroup_EW(0); map nf = hap_freq_table; double diff = 0; flag = 1; map::iterator iter=hap_freq_table.begin(); while(iter!=hap_freq_table.end()){ diff += (nf[iter->first]-of[iter->first])*(nf[iter->first]-of[iter->first]); if(nf[iter->first]<0.0000001){ flag = 0; break; } iter++; } if(diff<1e-5 || count>5 || flag==0) break; count++; } unknown_hap_freq_table_EW = hap_freq_table; } reset_freq_table(full_hap_freq_table); } void IQLBEstimator::print_subgroup_EW(vector &a0v, vector &a1v){ fprintf(outfp, "Haplotype frequency estimates by phenotype category (EW_b estimator)\n"); fprintf(outfp, "Haplotype\tOverall\t"); if(aff>=15) fprintf(outfp, "\tAffected"); if(unaff>=15) fprintf(outfp, "\tUnaffected"); if(unknown>=15) fprintf(outfp, "\tUnknown"); fprintf(outfp, "\n"); double sum=0, sum_aff=0, sum_unaff=0, sum_unknown=0; int ind=0, ind_aff=0, ind_unaff=0, ind_unknown=0; map::iterator iter = hap_freq_table.begin(); while(iter!=hap_freq_table.end()){ if(full_hap_freq_table_EW[iter->first]>=0) sum+=full_hap_freq_table_EW[iter->first]; else ind=1; if(aff>=15){ if(aff_hap_freq_table_EW[iter->first]>=0) sum_aff+=aff_hap_freq_table_EW[iter->first]; else ind_aff=1; } if(unaff>=15){ if(unaff_hap_freq_table_EW[iter->first]>=0) sum_unaff+=unaff_hap_freq_table_EW[iter->first]; else ind_unaff=1; } if(unknown>=15){ if(unknown_hap_freq_table_EW[iter->first]>=0) sum_unknown+=unknown_hap_freq_table_EW[iter->first]; else ind_unknown=1; } iter++; } iter = hap_freq_table.begin(); while(iter!=hap_freq_table.end()){ translate_hap(outfp,iter->first,loci_num,a0v,a1v); if(loci_num==2) fprintf(outfp, "\t"); if(full_hap_freq_table_EW[iter->first]>=0){ if(ind==0) fprintf(outfp, "\t%.4f\t", full_hap_freq_table_EW[iter->first]); else fprintf(outfp, "\t%.4f\t", full_hap_freq_table_EW[iter->first]/sum); } else fprintf(outfp, "\t0.0000\t"); if(aff>=15){ if(aff_hap_freq_table_EW[iter->first]>=0){ if(ind_aff==0) fprintf(outfp, "\t%.4f\t", aff_hap_freq_table_EW[iter->first]); else fprintf(outfp, "\t%.4f\t", aff_hap_freq_table_EW[iter->first]/sum_aff); } else fprintf(outfp, "\t0.0000\t"); } if(unaff>=15){ if(unaff_hap_freq_table_EW[iter->first]>=0){ if(ind_unaff==0) fprintf(outfp, "\t%.4f\t", unaff_hap_freq_table_EW[iter->first]); else fprintf(outfp, "\t%.4f\t", unaff_hap_freq_table_EW[iter->first]/sum_unaff); } else fprintf(outfp, "\t0.0000\t"); } if(unknown>=15){ if(unknown_hap_freq_table_EW[iter->first]>=0){ if(ind_unknown==0) fprintf(outfp, "\t%.4f", unknown_hap_freq_table_EW[iter->first]); else fprintf(outfp, "\t%.4f", unknown_hap_freq_table_EW[iter->first]/sum_unknown); } else fprintf(outfp, "\t0.0000"); } fprintf(outfp, "\n"); iter++; } fprintf(outfp, "\n\n"); }