diff --git a/KEEPfileinputnew b/KEEPfileinputnew
new file mode 100644
index 0000000000000000000000000000000000000000..80a445a869fd1fce38f525486b73234632894d1a
Binary files /dev/null and b/KEEPfileinputnew differ
diff --git a/allass.h b/allass.h
index ba153d23ed4a3033af3f16f69c617429182b966e..9f6278ac844a0e2f9d4f116d899f7e49f4e101dc 100755
--- a/allass.h
+++ b/allass.h
@@ -1,14 +1,11 @@
-
-
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
-#define MAX_INDS 2000
-#define MAX_LOCI 150000
-#define MAX_LINE 500000
+#define MAX_LOCI 100000
+#define MAX_LINE 50000
#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
@@ -18,6 +15,8 @@
#define NDIV (1+(IM-1)/NTAB)
#define EPS2 1.2e-7
#define RNMX (1.0-EPS2)
+#define SQR(a) (a)*(a)
+#define FMAX(a,b) ((a) > (b)? (a): (b))
/*>>>>>>>>DEFINE GLOBAL VARIABLES, STRUCTURES, & TYPES<<<<<<<<<*/
/*LOCI AND LDU DATA FOR MERGING ------------------------------------------*/
@@ -37,81 +36,6 @@ typedef lds *ldsPtr;
/*FIELDS---------------------------------------------------*/
-struct field{
-char locus[20];
-int col1;
-int col2;
-struct field *nextPtr;
-};
-
-typedef struct field fields;
-typedef fields *fieldsPtr;
-
-/*Haplotypes : up to loci ------------------------------------------*/
-
-struct hap{
-char allele[MAX_LOCI][6];
-double count;
-int flag;
-struct hap *nextPtr;
-};
-
-typedef struct hap haps;
-typedef haps *hapsPtr;
-
-/*For pooling haplotype counts - whole file ------------------------------------------*/
-
-struct all{
-char locus1[15];
-char locus2[15];
-double kb1;
-double kb2;
-double a;
-double b;
-double c;
-double d;
-int flag;
-struct all *nextPtr;
-};
-
-typedef struct all alls;
-typedef alls *allsPtr;
-
-
-
-/*For pooling haplotype counts ------------------------------------------*/
-/*
-struct poo{
-char locus1[15];
-char locus2[15];
-double kb1;
-double kb2;
-double a;
-double b;
-double c;
-double d;
-struct poo *nextPtr;
-};
-
-typedef struct poo poos;
-typedef poos *poosPtr;
-*/
-/*Allele frequencies -------------------------------------------------*/
-
-struct al{
-char locus[20];
-int index;
-int idis;
-char allele[6];
-double count;
-double count2;
-double it;
-double freq;
-struct al *nextPtr;
-};
-typedef struct al als;
-typedef als *alsPtr;
-
/*abcd table ----------------------------------------------------------*/
struct ab{
@@ -126,9 +50,6 @@ double ldu;
double ldu1;
double ldu2;
double n;
-double chi;
-double freq1;
-double freq2;
int i1;
int i2;
int flag;
@@ -141,15 +62,11 @@ typedef abcds *abcdPtr;
/*atoi table (of 9 counts for diallelic diplotypes) */
struct ai{
-char locus[20];
char locus1[20];
char locus2[20];
-char ckb[25];
-char ckb1[25];
-char ckb2[25];
double kb1;
double kb2;
-double aitab[50];
+int aitab[15];
double n;
double freq1;
double freq2;
@@ -190,32 +107,19 @@ typedef hugs *hugsPtr;
/*Mapc data---------------------------------------------------*/
struct map2{
-char locus[15];
+char locus[20];
double kb;
int order;
-int col1;
-int col2;
double freq;
+double chi;
struct map2 *nextPtr;
};
typedef struct map2 mapkb2;
typedef mapkb2 *mapPtr2;
-/*Input fields---------------------------------------------------*/
-struct fields{
-char field[15];
-int col1;
-int col2;
-struct fields *nextPtr;
-};
-typedef struct fields fields1;
-typedef fields1 *fieldPtr;
/*ped data---------------------------------------------------*/
struct ped{
-char FAM[15];
-char ID[15];
-int N;
int GEN[MAX_LOCI][2];
int numg;
struct ped *nextPtr;
@@ -226,137 +130,93 @@ typedef peds *pedPtr;
/*********************************************************************************/
/*FILES.....*/
FILE *output_f;
-FILE *output_fi;
-FILE *output_f5;
-FILE *output_f3;
-FILE *output_f2;
-FILE *fped;
FILE *output_ft;
-FILE *output_pf;
FILE *job_fp;
+FILE *output_f2;
+FILE *fped;
/*POINTERS.....*/
ldsPtr ldPtr,ldsstartPtr;
-fieldsPtr g_p1,fields_startPtr;
-hapsPtr hap_startPtr,g_hp1;
-alsPtr alstartPtr;
abcdPtr abcdstartPtr,g_abcd;
aisPtr aistartPtr,g_aiPtr,g_ai2Ptr;
-allsPtr all_startPtr;
mapPtr2 g_p12,map_startPtr2;
-pedPtr ped_startPtr;
-fieldPtr g_f1,field_startPtr;
+pedPtr g_ped1,ped_startPtr;
+hugsPtr bigstart;
/********************/
-double g_lnl,g_location[MAX_LOCI][3];
+double g_maf, g_hwp;
+double g_flag, g_chi1,g_prob,global_p, g_lnl,g_location[MAX_LOCI][3];
double g_alllocation[MAX_LOCI][4];
double global_m,global_l,global_e;
-double g_kee,g_kll,g_ktt,g_kss,g_kel,g_ket,g_kem,g_kes,g_klt,g_klm,g_kls,g_kts,g_kms;
+double g_kel,g_kem,g_klm;
double g_uei,g_uli,g_umi,g_pred;
-double g_meank,g_meanchi,g_startlnl;
+double g_startlnl;
double interv[MAX_LOCI][7],g_rho,g_rhoi;
double cinterv[MAX_LOCI][7];
-double g_best,g_sumki2,g_sumki,g_cov,g_covi;
-double g_seE,g_seL,g_seM,g_bestlik,g_bestE,g_bestL,g_bestM;
-double g_bestseE,g_bestseL,g_bestseM;
+double g_best,g_sumki2,g_sumki;
+double g_max,g_sed;
+double g_df,g_V,g_ekeep,g_mkeep;
/********************/
-char ped_file[20];
-char g_datfile[15];
-char g_posneg[6];
char g_temp[MAX_LINE];
char gg_temp[500];
char terfile[20],intefile[20],datfile[20];
char g_loci[MAX_LOCI][20];
char g_allloci[MAX_LOCI][20];
-char g_disease[15],g_jobfile[20];
/********************/
-int g_bad,g_n,g_nloci,g_nallloci,g_ind,g_niter,g_finish;
-int ldflag2,ldflag;
+int g_bad,g_n,g_nloci,g_nallloci,g_niter,g_finish;
+int ldflag;
int itlast, iteast,itmast,g_nki;
int g_iter[5], global_fin,g_calls,g_ind;
-int ite,itl,itm,itr,itd,g_sig;
+int ite,itl,itm,g_sig;
int g_estimatel;
+int g_opt,path;
+int g_npar,g_int;
+int g_nhole;
+int g_initial;
+float ran1(long * );
+long idum;
/********************/
-
-void chi(double,double);
-void final();
+double **dmatrix(long nrl, long nrh, long ncl, long nch);
void jobin();
-void file_input(char *);
-void pool(char *);
-void file_input2(char *);
-void make_map();
-void getfields();
-void make_hap();
-void get_tab();
void getpa(char[500]);
-double gammln(double);
-double *dvector(long nl, long nh);
-double **dmatrix(long nrl, long nrh, long ncl, long nch);
-int **imatrix(long nrl, long nrh, long ncl, long nch);
-void runewt();
-void cal_SE2(double **,int);
-void runewt2();
-void intfile2();
-void intfile3();
-#define SQR(a) (a)*(a)
-#define FMAX(a,b) ((a) > (b)? (a): (b))
-void nrerror();
-double *dvector();
-double **dmatrix();
int *ivector();
+double bound(int,double[]);
void free_dmatrix();
+void cal_SE2(double **,int);
+void get_tab();
+void getfields();
double func(double[]);
-double bound(int,double[]);
void dfunc(double[],double[]);
-void writeint();
-void metricho(double,double,double,double);
+double *dvector(long nl, long nh);
+void nrerror();
+void runewt2();
void diplo(char *);
-void matrix();
-void fill_str2(char[20],char[20]);
-void update();
void ludcmp(double **,int,int *,double *);
-void meankchi(char[20]);
-void get_predl();
+void fill_gai(int, int, int, int, int, int, int, int, int, double, double, char[20], char[20],double, double);
void make_map4();
void readped();
-void writ();
-void keep_it();
-void keep_it2();
+void metricho(double,double,double,double);
+void meankchi(char[20]);
+void update();
void one_interval(int , intsPtr,double *, double *, double *);
-void locus_list();
-int g_nhole;
-double g_max,g_pp,g_ppi,g_sed;
+void keep_it2();
void fast();
-hugsPtr bigstart;
-int g_npar,g_int;
-double g_df,g_V,g_ekeep,g_mkeep;
-void directmap();
-int g_code,g_opt,path;
-void hapallele();
-void haplo();
-char g_SNP[20];
void quicklike(double,double,double);
-int g_initial;
-void writeter();
-void writeterfin();
-void printinputmap();
-void reorder();
-void reorderall();
+void keep_it();
+void get_predl();
+void locus_list();
void directmap2();
void updat1();
+void reorderall();
void runewt3();
-void file_inputg();
-void multallele1();
-void diplo1(char *);
-void getfreqsh();
-void getfreqsh2();
-void hapallele2();
+void reorder();
+void newentry();
+void newentry2();
+void file_inputnew();
void getfreqs();
void multallele();
void segments();
-void ldmap2();
-void mergeloc();
void lubksb(double **, int, int *, double []);
void intfile4();
void intfileseg();
@@ -364,18 +224,10 @@ void ldmapseg();
void writeterfinseg();
void intfile6();
void free_dvector(double *, long, long);
-void printmap();
-float ran1(long * );
-char g_buffer[MAX_INDS][MAX_LINE];
-int g_used[MAX_INDS];
-long idum;
-int g_nrec;
-int g_number;
-
void dfpmin(double p[], int n, double gtol, int *iter, double *fret,
double(*func)(double []), void (*dfunc)(double [], double []));
-
-void intfile7();
-void file_inputg2();
-
-void simboot();
+double gammln(double);
+void gcf(double *,double,double,double *,double *);
+double gammp(double,double);
+void gser(double *,double,double,double *,double *);
+void chi(double,double);
diff --git a/cal_SE2.c b/cal_SE2.c
index b3aaed4a540fb986683113d470c03bb880d8cf38..c93a774275a0f035927bda63406bb1e3ad6d3a48 100755
--- a/cal_SE2.c
+++ b/cal_SE2.c
@@ -9,10 +9,6 @@ double big=9999999.99999,**temp, **y, d, *col;
char blank[14];
if(ldflag==2)return;
-
-g_seE=999.999;
-g_seM=999.999;
-g_seL=999.999;
strcpy(blank," ");
indx=ivector(1,n);
col=dvector(1,n);
@@ -36,29 +32,24 @@ for(j=1; j<=n; j++)
j=1;
fprintf(output_f,"\n\nStdrd Errors:");
-if(g_iter[1]==1&&(y[j][j]>=0.)&&(y[j][j]<big)){g_seE=sqrt(y[j][j]);fprintf(output_f,"%13.5f",sqrt(y[j][j]));j++;
+if(g_iter[1]==1&&(y[j][j]>=0.)&&(y[j][j]<big)){fprintf(output_f,"%13.5f",sqrt(y[j][j]));j++;
}
else
{
fprintf(output_f," %s",blank);
}
-if(g_iter[2]==1&&(y[j][j]>=0.)&&(y[j][j]<big)){g_seL=sqrt(y[j][j]);fprintf(output_f," %13.5f",sqrt(y[j][j]));j++;}
+if(g_iter[2]==1&&(y[j][j]>=0.)&&(y[j][j]<big)){fprintf(output_f," %13.5f",sqrt(y[j][j]));j++;}
else
{
fprintf(output_f," %s",blank);
}
-if(g_iter[3]==1&&(y[j][j]>=0.)&&(y[j][j]<big)){g_seM=sqrt(y[j][j]);fprintf(output_f," %13.5f",sqrt(y[j][j]));j++;}
+if(g_iter[3]==1&&(y[j][j]>=0.)&&(y[j][j]<big)){fprintf(output_f," %13.5f",sqrt(y[j][j]));j++;}
else
{
fprintf(output_f," %s",blank);
}
-/*fprintf(output_f,"\n\n Kel Kem Klm ");
-fprintf(output_f,"\n%13.5f %13.5f %13.5f ",g_kel,g_kem,g_klm);
-*/
-
fprintf(output_f,"\n\nN (Number of pairs)= %7d Predicted L (Lp)= %13.5f",
g_nki,sqrt(2./3.14159265)*(g_sumki2/g_nki)/(g_sumki/g_nki));
-
}
diff --git a/chi.c b/chi.c
new file mode 100644
index 0000000000000000000000000000000000000000..f4b43a5e9f10645f7496f52ec6f9acd845a81b36
--- /dev/null
+++ b/chi.c
@@ -0,0 +1,135 @@
+#include "allass.h"
+
+void chi(double df, double chi)
+{
+double thep,prob,chi1;
+double xp,c0,c1,c2,d1,d2,d3,t,t2,t3;
+
+t=gammp(df,chi);
+prob=1.-global_p;
+c0 = 2.515517;
+c1 = 0.802853;
+c2 = 0.010328;
+d1 = 1.432788;
+d2 = 0.189269;
+d3 = 0.001308;
+thep = prob/2;
+if(g_flag==1)t=sqrt(log(1./(thep*thep)));
+t2 = t*t;
+t3 = t2*t;
+
+xp = t - (c0+c1*t+c2*t2)/ (1.+d1*t+d2*t2+d3*t3) ;
+chi1 = xp*xp;
+
+g_chi1=chi1;
+g_prob=2.*thep;
+}
+
+/*-----------------------------------------------------------------------------------
+Function name: gammp
+Date last modified: 7/2/95
+Description: Chi square function
+Variables passed: a and x
+Global variables used: -
+Variables returned: result
+------------------------------------------------------------------------------------*/
+
+double gammp(double a , double x)
+
+{
+
+void gcf(double *gammcf,double a, double x,double *gln,double *bigp);
+void gser(double *gamser,double a, double x, double *gln, double *bigp);
+double gamser,gammcf,gln,bigp;
+a=a/2.;
+x=x/2.;
+
+if(x<0.0 || a<=0.0){printf("\nInvalid arguments in gammp");}
+if(x<(a+1.0)){
+gser(&gamser,a,x,&gln,&bigp);
+global_p = bigp;
+g_flag=1;
+return gamser;
+
+} else {
+ gcf(&gammcf,a,x,&gln,&bigp);
+ global_p = 1.-bigp;
+ g_flag=0;
+ return gammcf;
+
+}
+}
+
+#define ITMAX 50
+#define EPS 3.0e-7
+void gser(double *gamser,double a,double x,double *gln,double *bigp)
+{
+double gammln(double xx);
+int n;
+double sum,del,ap;
+*gln=gammln(a);
+if(x<=0.0){
+if(x<0.0){printf("x less than 0 in routine gser");return;}
+
+*gamser=0.0;
+return;
+}else{
+ap=a;
+del=sum=1.0/a;
+for(n=1;n<=ITMAX;n++){
+++ap;
+del *= x/ap;
+sum += del;
+if(fabs(del)<fabs(sum)*EPS){
+*bigp=sum*exp(-x+a*log(x)-*gln);
+return;
+}}
+printf("a too large, ITMAX too small in routine gser");
+return;
+}
+}
+
+#define ITMAX 50
+#define EPS 3.0e-7
+#define FPMIN 1.0e-300
+void gcf(double *gammcf,double a,double x,double *gln,double *bigp)
+{
+ double gammln(double xx);
+ int i;
+ double an,b,c,d,del,h;
+ *gln=gammln(a);
+ b=x+1.0-a;
+ c=1.0/FPMIN;
+ d=1.0/b;
+h=d;
+ for(i=1;i<=ITMAX;i++){
+ an=-i*(i-a);
+ b+=2.0;
+ d=an*d+b;
+ if(fabs(d)<FPMIN)d=FPMIN;
+ c=b+an/c;
+ if(fabs(c)<FPMIN)c=FPMIN;
+ d=1.0/d;
+ del=d*c;
+ h*=del;
+ if(fabs(del-1.0)<EPS) break;
+}
+if(i>ITMAX){printf("a too large, ITMAX too small in gcf");}
+*bigp=exp(-x+a*log(x)-(*gln))*h;
+*gammcf = sqrt((-2.* (log(h)-x + a*log(x)- *gln ))+log(4.));
+}
+
+double gammln(double xx)
+{
+ double x,y,tmp,ser;
+static double cof[6]={76.18009172947146,-86.50532032941677,
+24.01409824083091,-1.231739572450155,
+0.1208650973866179e-2,-0.5395239384953e-5};
+int j;
+y=x=xx;
+tmp=x+5.5;
+tmp -=(x+0.5)*log(tmp);
+ser=1.000000000190015;
+for(j=0;j<=5;j++)ser+=cof[j]/++y;
+return -tmp+log(2.5066282746310005*ser/x);
+}
diff --git a/diplo.c b/diplo.c
index c2dc382fd9bbb8f6cecc1fa56f27fe65dbc411a5..980e1eb009a804c128fd6b93e532a9944ae150ef 100755
--- a/diplo.c
+++ b/diplo.c
@@ -5,21 +5,18 @@ void diplo(char *outputfile)
/*Needs 9 counts from a 3*3 genotype table */
char xflag=' ';
int numy=0,nrec=0,times,flag,kk;
-double nd,kb,kb2,kb1,f11o,D,x11,x12,x21,x22,f11,f12,f21,f22,temp,sum,n11,n12,n13,n21,n22,n23,n31,n32,n33,Q,R;
+double kb,kb2,kb1,f11o,D,x11,x12,x21,x22,f11,f12,f21,f22,temp,sum,n11,n12,n13,n21,n22,n23,n31,n32,n33,Q,R;
double a,b,c,d;
output_f2=NULL;
-output_f3=NULL;
if((output_f2=fopen(outputfile,"w"))==NULL)
{printf("\nCannot open intermediate file\n");exit(1);}
-/*hill.out file gives the haplotype counts A,B,C,D - can be used for pooling across
-populations*/
g_aiPtr=aistartPtr;
while(g_aiPtr!=NULL)
{
-kb1=atof(g_aiPtr->ckb1);
-kb2=atof(g_aiPtr->ckb2);
+kb1=g_aiPtr->kb1;
+kb2=g_aiPtr->kb2;
kb=fabs(kb1-kb2);
sum=0;
for(kk=0;kk<9;kk++){ if(kk<9)sum=sum+g_aiPtr->aitab[kk]; }
@@ -72,7 +69,6 @@ if(times<15000&&fabs(f11o-f11)>0.000000001) { f11o=f11; goto top; }
/**********************************************/
/*CHECK, GET f11 under E.M. algorithm*/
-nd=sum;
numy++;
/**********************************************/
a=f11;
@@ -114,16 +110,13 @@ metricho(Q,R,sum,D);
if(g_rho>1.0)g_rho=1.0;
if(g_rho<0.0)g_rho=0.0;
/***********************************/
-
if(nrec>0)
{
-fprintf(output_f2,"\n%-15s %-15s %10.3f %10.3f %12.10f %10.3f %8.2f %5.0f %12.10f %12.10f %12.10f max_int %6d max_kb %8.2f %12.10f %12.10f",
-g_aiPtr->locus1,g_aiPtr->locus2,kb1,kb2,g_rho,g_rhoi,(g_rho*g_rho)*g_rhoi,sum,Q,R,D,g_int,g_max,g_aiPtr->freq1,g_aiPtr->freq2);
+fprintf(output_f2,"\n%-15s %-15s %10.3f %10.3f %12.10f %10.3f " , g_aiPtr->locus1,g_aiPtr->locus2,kb1,kb2,g_rho,g_rhoi);
}
if(nrec==0)
{
-fprintf(output_f2,"%-15s %-15s %10.3f %10.3f %12.10f %10.3f %8.2f %5.0f %12.10f %12.10f %12.10f max_int %6d max_Kb %8.2f %12.10f %12.10f",
-g_aiPtr->locus1,g_aiPtr->locus2,kb1,kb2,g_rho,g_rhoi,(g_rho*g_rho)*g_rhoi,sum,Q,R,D,g_int,g_max,g_aiPtr->freq1,g_aiPtr->freq2);
+fprintf(output_f2,"%-15s %-15s %10.3f %10.3f %12.10f %10.3f ", g_aiPtr->locus1,g_aiPtr->locus2,kb1,kb2,g_rho,g_rhoi);
}
nrec++;
diff --git a/file_inputnew.c b/file_inputnew.c
new file mode 100644
index 0000000000000000000000000000000000000000..87ca58c7e37a553751a6d9d67098cf2a40e4f951
--- /dev/null
+++ b/file_inputnew.c
@@ -0,0 +1,114 @@
+#include "allass.h"
+
+void file_inputnew()
+{
+/*READ IN GENOTYPE FILE*/
+/*This is expecting a genotype file in PLINK TPED format which has chromosome number, SNP name, genetic distance in cM, and base pair position followed by geontypes for all individuals in one line. Therefore very long files of many SNPs can be handled. The chrom number and cM locations are not used by this program */
+char buffer[MAX_LINE];
+int j,iset,tok, alkeep1, alkeep2, ii, ikeep=0,iii,nloci=0;
+char temp1[50];
+int i;
+
+g_ped1=NULL;
+g_p12=NULL;map_startPtr2=NULL;
+nloci=0;
+if((fped=fopen(datfile,"r"))==NULL){printf("\n Cannot open genotype file ");exit(1);}
+/******************************************************************************************/
+/******************************************************************************************/
+/******************************************************************************************/
+while(fgets(buffer, MAX_LINE, fped)!=NULL)
+ {
+ for(iii=0;iii<MAX_LINE;iii++){if(buffer[iii]=='\0'){ikeep=iii;goto lab1;} }
+ /*************************************************************************/
+ /*************************************************************************/
+lab1: buffer[ikeep+1]=' ';
+ if(nloci==0)
+ {
+ for(ii=0;ii<50;ii++){temp1[ii]='\0';}
+ j=0;
+ iset=0;
+ alkeep1=0;
+ alkeep2=0;
+ tok=0;
+ /*****/
+ for(i=0;i<=ikeep;i++)
+ {
+ if((buffer[i]!=' ')&&(buffer[i+1]!=' ')) { temp1[j]=buffer[i];j++; }
+ if((buffer[i]!=' ')&&(buffer[i+1]==' '))
+ {
+ temp1[j]=buffer[i];j++;
+ tok++;
+ if(tok==2)
+ {
+ newentry2();
+ strcpy(g_p12->locus, temp1);
+ }
+ if(tok==4)
+ {
+ g_p12->kb=atof(temp1)/1000.;
+ }
+ if(tok>4)
+ { if(iset==0){alkeep1=atoi(temp1);iset=1;goto skip;}
+ if(iset==1){alkeep2=atoi(temp1);newentry();
+ g_ped1->GEN[nloci][0]=alkeep1; g_ped1->GEN[nloci][1]=alkeep2;iset=0; }
+ }
+skip:; for(ii=0;ii<50;ii++){temp1[ii]='\0';}
+ j=0;
+ }
+
+ }/*for loop*/
+/***IF FIRST LINE USE TO CREATE WHOLE PED STRUCTURE */
+ } /*if nloci==0*/
+
+ /*************************************************************************/
+ /*************************************************************************/
+ if(nloci>0)
+ {
+ g_ped1=ped_startPtr;
+ for(ii=0;ii<50;ii++){temp1[ii]='\0';}
+ j=0;
+ iset=0;
+ alkeep1=0;
+ alkeep2=0;
+ tok=0;
+ /*****/
+ for(i=0;i<=ikeep;i++)
+ {
+ if((buffer[i]!=' ')&&(buffer[i+1]!=' ')) { temp1[j]=buffer[i];j++; }
+ if((buffer[i]!=' ')&&(buffer[i+1]==' '))
+ {
+ temp1[j]=buffer[i];j++;
+ tok++;
+ if(tok==2)
+ {
+ newentry2();
+ strcpy(g_p12->locus, temp1);
+ }
+ if(tok==4)
+ {
+ g_p12->kb=atof(temp1)/1000.;
+ }
+ if(tok>4)
+ {
+ if(iset==0){alkeep1=atoi(temp1);iset=1;goto skips;}
+ if(iset==1){alkeep2=atoi(temp1);
+ g_ped1->GEN[nloci][0]=alkeep1;
+ g_ped1->GEN[nloci][1]=alkeep2;g_ped1=g_ped1->nextPtr;iset=0; }
+ }
+skips:; for(ii=0;ii<50;ii++){temp1[ii]='\0';}
+ j=0;
+ }
+
+ }/*for loop*/
+ } /* if nloci >0*/
+ /*************************************************************************/
+ /*************************************************************************/
+ nloci++;
+ } /*fgets*/
+
+
+
+g_ped1=ped_startPtr;
+g_ped1->numg=nloci-1;
+}
+
diff --git a/fill_gai.c b/fill_gai.c
new file mode 100644
index 0000000000000000000000000000000000000000..d70edacdd283667248b8733e5789562827e79460
--- /dev/null
+++ b/fill_gai.c
@@ -0,0 +1,54 @@
+#include "allass.h"
+void fill_gai(int n11, int n12, int n13, int n21, int n22, int n23, int n31, int n32, int n33, double kb1, double kb2, char locus1[20], char locus2[20], double freq1, double freq2)
+{
+double n;
+n=n11+n12+n13+n21+n22+n23+n31+n32+n33;
+
+if(g_aiPtr==NULL)
+ {
+ aistartPtr=(aisPtr)malloc(sizeof(ais));
+ aistartPtr->nextPtr=NULL;
+ g_aiPtr=aistartPtr;
+ g_aiPtr->kb1=kb1;
+ g_aiPtr->kb2=kb2;
+ g_aiPtr->n=n;
+ g_aiPtr->aitab[0]=n11;
+ g_aiPtr->aitab[1]=n12;
+ g_aiPtr->aitab[2]=n13;
+ g_aiPtr->aitab[3]=n21;
+ g_aiPtr->aitab[4]=n22;
+ g_aiPtr->aitab[5]=n23;
+ g_aiPtr->aitab[6]=n31;
+ g_aiPtr->aitab[7]=n32;
+ g_aiPtr->aitab[8]=n33;
+ strcpy(g_aiPtr->locus1, locus1);
+ strcpy(g_aiPtr->locus2, locus2);
+ g_aiPtr->freq1=freq1;
+ g_aiPtr->freq2=freq2;
+ }
+
+
+else /* If not the first record in the ais structure */
+ {
+ g_aiPtr->nextPtr=(aisPtr)malloc(sizeof(ais));
+ g_aiPtr=g_aiPtr->nextPtr;
+ g_aiPtr->kb1=kb1;
+ g_aiPtr->kb2=kb2;
+ g_aiPtr->n=n;
+ g_aiPtr->aitab[0]=n11;
+ g_aiPtr->aitab[1]=n12;
+ g_aiPtr->aitab[2]=n13;
+ g_aiPtr->aitab[3]=n21;
+ g_aiPtr->aitab[4]=n22;
+ g_aiPtr->aitab[5]=n23;
+ g_aiPtr->aitab[6]=n31;
+ g_aiPtr->aitab[7]=n32;
+ g_aiPtr->aitab[8]=n33;
+ strcpy(g_aiPtr->locus1, locus1);
+ strcpy(g_aiPtr->locus2, locus2);
+ g_aiPtr->freq1=freq1;
+ g_aiPtr->freq2=freq2;
+
+ }
+
+ }
diff --git a/getfreqs.c b/getfreqs.c
index 52bc5abc0ca4cacdc35e433cd06ebc4920f9dc2d..49d07358b4ad119931442ab2b36162c81e4042a3 100755
--- a/getfreqs.c
+++ b/getfreqs.c
@@ -2,19 +2,21 @@
void getfreqs()
{
-int r,nr,N,
+int j,nr1,nr2,r,nr,N,
numg,nals[MAX_LOCI],i,altab[2][2][MAX_LOCI];
pedPtr ped1;
mapPtr2 p1;
-double fmin,f;
+double chis[MAX_LOCI],df=1.,a,b,c,xchi,etot,tot,f,n1=0,n2=0,f1,f2,e11,e12,e22,n11,n12,n22;
/*numg is the number of loci in the file */
/*nals is the rth allele at locus i */
/*altab - r alleles, 0=allele label, 1 = allele count ith locus */
-
+/*N is the count of all alleles */
+/*nr is the number of rare alleles */
for(i=0;i<MAX_LOCI;i++)
{
nals[i]=0;
+ chis[i]=0.;
altab[0][0][i]=0;
altab[0][1][i]=0;
altab[1][0][i]=0;
@@ -26,46 +28,90 @@ numg=ped1->numg;
/******************************************************************************/
for(i=0;i<=numg;i++)
{
+xchi=tot=etot=0.;
+n11=n12=n22=0;
ped1=ped_startPtr;
while(ped1!=NULL)
{
- if(ped1->N==999)ped1->N=1;
+ /***Count genotypes******/
+ if((ped1->GEN[i][0]!=0)&&(ped1->GEN[i][1]!=0))
+ {
+ if((ped1->GEN[i][0]==1)&&(ped1->GEN[i][1]==1))n11=n11+1.;
+ if((ped1->GEN[i][0]==1)&&(ped1->GEN[i][1]==2))n12=n12+1.;
+ if((ped1->GEN[i][0]==2)&&(ped1->GEN[i][1]==1))n12=n12+1.;
+ if((ped1->GEN[i][0]==2)&&(ped1->GEN[i][1]==2))n22=n22+1.;
+ }
+
+ /***Count genotypes******/
for(r=0;r<2;r++)
{
if(ped1->GEN[i][0]!=0)
- {if(altab[r][0][i]==ped1->GEN[i][0]) {altab[r][1][i]=altab[r][1][i]+ped1->N;goto lab1;}}
+ {if(altab[r][0][i]==ped1->GEN[i][0]) {altab[r][1][i]++;goto lab1;}}
}
if(ped1->GEN[i][0]!=0)
{altab[nals[i]][0][i]=ped1->GEN[i][0];
- altab[nals[i]][1][i]=ped1->N;
+ altab[nals[i]][1][i]=1;
nals[i]++;}
+
lab1: for(r=0;r<2;r++)
{
if(ped1->GEN[i][1]!=0)
- {if(altab[r][0][i]==ped1->GEN[i][1]) {altab[r][1][i]=altab[r][1][i]+ped1->N;goto lab2;}}
+ {if(altab[r][0][i]==ped1->GEN[i][1]) {altab[r][1][i]++;goto lab2;}}
}
if(ped1->GEN[i][1]!=0)
- {altab[nals[i]][0][i]=ped1->GEN[i][1]; altab[nals[i]][1][i]=ped1->N; nals[i]++; }
+ {altab[nals[i]][0][i]=ped1->GEN[i][1];
+ altab[nals[i]][1][i]=1;
+ nals[i]++; }
+
lab2: ped1=ped1->nextPtr;
}
-}
+tot=n11+n12+n22;
+if((altab[0][0][i])==1) {n1=altab[0][1][i]; n2=altab[1][1][i]; }
+if((altab[1][0][i])==1) {n2=altab[0][1][i]; n1=altab[1][1][i]; }
+f1=n1/(n1+n2);
+f2=n2/(n1+n2);
+e11=(f1*f1)*tot;
+e12=2*(f1*f2)*tot;
+e22=(f2*f2)*tot;
+etot=e11+e12+e22;
+if(e11==0)e11=0.5;
+if(e12==0)e12=0.5;
+if(e22==0)e22=0.5;
+a=((n11-e11)*(n11-e11))/e11;
+b=((n12-e12)*(n12-e12))/e12;
+c=((n22-e22)*(n22-e22))/e22;
+xchi=a+b+c;
+g_prob=1.;
+if(xchi>0.)chi(df,xchi);
+/*fprintf(output_f,"\nExp: %10.3f %10.3f %10.3f %10.3f Obs: %10.3f %10.3f %10.3f %10.3f Chi: %10.3f P=%10.6f ", e11,e12,e22,etot,n11,n12,n22,tot,xchi,g_prob);
+*/
+chis[i]=g_prob;
+
+} /*outer loop over all markers*/
/*****************************************************************************/
-/*LOOP OVER ALL LOCI */
+fprintf(output_f,"\n\nTotal number of loci in file=%d\n",numg+1);
+fprintf(output_f,"\nCut-off used for MAF = %f",g_maf);
+fprintf(output_f,"\nCut-off used for Hardy-Weinberg test (Pvalue) = %f ",g_hwp);
+fprintf(output_f,"\nLoci excluded with low MAF and/or significant Hardy-Weinberg test ");
+fprintf(output_f,"\n Locus Kb N_min_als N_als MAF HW_Pvalue ");
p1=map_startPtr2;
+j=0;
for(i=0;i<=numg;i++)
{
- fmin=1.;
N=0;
for(r=0;r<2;r++)N=N+altab[r][1][i];
-for(r=0;r<2;r++)
- {
- nr=altab[r][1][i];
- f=(double)nr/(double)N;
- if (f<fmin)fmin=f;
-
- }
-p1->freq=fmin;
+ nr1=altab[0][1][i];
+ nr2=altab[1][1][i];
+ nr=nr1;
+ if(nr2<nr1)nr=nr2;
+ f=(double)nr/(double)N;
+ p1->freq=f;
+ p1->chi=chis[i];
+if((p1->freq<=g_maf)&&(p1->chi>g_hwp)){fprintf(output_f,"\n%20s %10.3f %5d %5d %10.4f ",p1->locus, p1->kb,nr,N,p1->freq);}
+if((p1->freq>g_maf)&&(p1->chi<=g_hwp)){fprintf(output_f,"\n%20s %10.3f %5d %5d %10.4f ",p1->locus, p1->kb,nr,N,p1->chi); }
+if((p1->freq<=g_maf)&&(p1->chi<=g_hwp)){fprintf(output_f,"\n%20s %10.3f %5d %5d %10.4f %10.4f ",p1->locus, p1->kb,nr,N,p1->freq,p1->chi); }
p1=p1->nextPtr;
}
+fflush(output_f);
}
diff --git a/intfile4.c b/intfile4.c
index d7a3e287e4262bd1545e1f8dd0262f1ad6746167..10b97a937ce1cf66d07f84e57988159e381de211 100755
--- a/intfile4.c
+++ b/intfile4.c
@@ -12,20 +12,19 @@ int ic=0,i,j;
char temp[20];
char buff[200];
char cbuff[200];
-char mint[30];
FILE *fp;
abcdPtr abcd_p1;
char locs[20],loc1[20],loc2[20];
-double f1,f2,kb1,kb2;
+double kb1,kb2;
g_nallloci=0;
-if((output_ft=fopen(intefile,"w"))==NULL){ printf("\nCannot open map (output) file");exit(1);}
-if((fp=fopen(datfile,"r"))==NULL){ printf("\nCannot open intermediate file for reading");exit(1);}
-fprintf(output_f,"\nReading the intermediate data file: %s ",datfile);
-fprintf(output_f,"\nWriting the map file: %s ",intefile);
-fprintf(output_ft,"\n#Reading the intermediate data file: %s ",datfile);
-fprintf(output_ft,"\n#Writing the map file: %s ",intefile);
+if((output_ft=fopen(terfile,"w"))==NULL){ printf("\nCannot open map (output) file");exit(1);}
+if((fp=fopen(intefile,"r"))==NULL){ printf("\nCannot open intermediate file for reading");exit(1);}
+fprintf(output_f,"\n\nReading the intermediate data file: %s ",intefile);
+fprintf(output_f,"\nWriting the map file: %s ",terfile);
+fprintf(output_ft,"\n#Reading the intermediate data file: %s ",intefile);
+fprintf(output_ft,"\n#Writing the map file: %s ",terfile);
abcd_p1=NULL;
@@ -34,33 +33,6 @@ while(fgets(buff,200,fp)!=NULL)
if(buff[0]=='\n')goto endit2;
ic++;
strcpy(cbuff,buff);
-/*******************/
- if(ic==1){
- fprintf(output_f,"\n\nFirst line of the intermediate file:");
- fprintf(output_f,"\n%s",cbuff);
- for(i=0;i<30;i++)mint[i]='\0';
-
- for(i=0;i<200;i++)
- {
- if(cbuff[i]=='n'&&cbuff[i+1]=='t')
- {
- mint[0]=cbuff[i+2]; mint[1]=cbuff[i+3]; mint[2]=cbuff[i+4]; mint[3]=cbuff[i+5]; mint[4]=cbuff[i+6];
- mint[5]=cbuff[i+7]; mint[6]=cbuff[i+8]; mint[7]=cbuff[i+9]; mint[8]='\0';
- }
- }
- fprintf(output_f,"\n\nThe maximum number of intervals for computing each epsilon is: %s ",mint);
-
- for(i=0;i<30;i++)mint[i]='\0';
- for(i=0;i<200;i++)
- {
- if(cbuff[i]=='K'&&cbuff[i+1]=='b')
- {
- mint[0]=cbuff[i+2]; mint[1]=cbuff[i+3]; mint[2]=cbuff[i+4]; mint[3]=cbuff[i+5]; mint[4]=cbuff[i+6];
- mint[5]=cbuff[i+7]; mint[6]=cbuff[i+8]; mint[7]=cbuff[i+9]; mint[8]=cbuff[i+10]; mint[9]=cbuff[i+11];mint[10]='\0';
- }
- }
- fprintf(output_f,"\n\nThe maximum distance (kb) between any pair is set to %s ",mint);
- }
/*******************/
for(i=0;i<20;i++) { loc1[i]='\0'; loc2[i]='\0'; }
kb1=kb2=0.;
@@ -87,24 +59,12 @@ j=0;
for(i=43;i<53;i++){temp[j]=buff[i];j++;}
kb2=atof(temp);
-for(i=0;i<20;i++)temp[i]='\0';
-j=0;
-for(i=163;i<175;i++){temp[j]=buff[i];j++;}
-f1=atof(temp);
-
-for(i=0;i<20;i++)temp[i]='\0';
-j=0;
-for(i=177;i<188;i++){temp[j]=buff[i];j++;}
-f2=atof(temp);
-
/*MAKE A LOCUS LIST*/
for(i=0;i<20;i++)locs[i]='\0';
for(i=0;i<20;i++) { if(loc1[i]=='.')goto gogo; locs[i]=loc1[i]; }
gogo:for(i=0;i<g_nallloci;i++) { if(strcmp(locs,g_allloci[i])==0)goto endit; }
strcpy(g_allloci[g_nallloci],locs);
g_alllocation[g_nallloci][1]=kb1;
-g_alllocation[g_nallloci][2]=f1;
-
g_nallloci++;
endit: for(i=0;i<20;i++)locs[i]='\0';
@@ -112,7 +72,6 @@ for(i=0;i<20;i++) { if(loc2[i]=='.')goto gogo2; locs[i]=loc2[i]; }
gogo2:for(i=0;i<g_nallloci;i++) { if(strcmp(locs,g_allloci[i])==0)goto endit2; }
strcpy(g_allloci[g_nallloci],locs);
g_alllocation[g_nallloci][1]=kb2;
-g_alllocation[g_nallloci][2]=f2;
g_nallloci++;
diff --git a/intfile6.c b/intfile6.c
index ae502540c9cb336a9778ac62340e6e04b073d3a4..c90329f95b4266790f36d1fe75e39511e94ce6bd 100755
--- a/intfile6.c
+++ b/intfile6.c
@@ -12,17 +12,15 @@ int ic=0,i,j;
char temp[20];
char buff[200];
char cbuff[200];
-char mint[30];
FILE *fp;
abcdPtr abcd_p1;
char locs[20],loc1[20],loc2[20];
-double f1,f2,chi,n,kb1,kb2,p,k;
+double chi,n,kb1,kb2,p,k;
g_nloci=0;
-
-if((fp=fopen(datfile,"r"))==NULL){ printf("\nCannot open intermediate file for reading");exit(1);}
-fprintf(output_f,"\nReading the intermediate data file: %s ",datfile);
-fprintf(output_f,"\nWriting the map file: %s ",intefile);
+if((fp=fopen(intefile,"r"))==NULL){ printf("\nCannot open intermediate file for reading");exit(1);}
+fprintf(output_f,"\nReading the intermediate data file: %s ",intefile);
+fprintf(output_f,"\nWriting the map file: %s ",terfile);
abcd_p1=NULL;
@@ -31,33 +29,6 @@ while(fgets(buff,200,fp)!=NULL)
if(buff[0]=='\n')goto skip;
ic++;
strcpy(cbuff,buff);
-/*******************/
- if(ic==1){
- fprintf(output_f,"\n\nFirst line of the intermediate file:");
- fprintf(output_f,"\n%s",cbuff);
- for(i=0;i<30;i++)mint[i]='\0';
-
- for(i=0;i<200;i++)
- {
- if(cbuff[i]=='n'&&cbuff[i+1]=='t')
- {
- mint[0]=cbuff[i+2]; mint[1]=cbuff[i+3]; mint[2]=cbuff[i+4]; mint[3]=cbuff[i+5]; mint[4]=cbuff[i+6];
- mint[5]=cbuff[i+7]; mint[6]=cbuff[i+8]; mint[7]=cbuff[i+9]; mint[8]='\0';
- }
- }
- fprintf(output_f,"\n\nThe maximum number of intervals for computing each epsilon is: %s ",mint);
-
- for(i=0;i<30;i++)mint[i]='\0';
- for(i=0;i<200;i++)
- {
- if(cbuff[i]=='K'&&cbuff[i+1]=='b')
- {
- mint[0]=cbuff[i+2]; mint[1]=cbuff[i+3]; mint[2]=cbuff[i+4]; mint[3]=cbuff[i+5]; mint[4]=cbuff[i+6];
- mint[5]=cbuff[i+7]; mint[6]=cbuff[i+8]; mint[7]=cbuff[i+9]; mint[8]=cbuff[i+10]; mint[9]=cbuff[i+11];mint[10]='\0';
- }
- }
- fprintf(output_f,"\n\nThe maximum distance (kb) between any pair is set to %s ",mint);
- }
/*******************/
for(i=0;i<20;i++) { loc1[i]='\0'; loc2[i]='\0'; }
kb1=kb2=p=k=chi=n=0.;
@@ -93,28 +64,6 @@ for(i=0;i<20;i++)temp[i]='\0';
j=0;
for(i=68;i<78;i++){temp[j]=buff[i];j++;}
k=atof(temp);
-
-for(i=0;i<20;i++)temp[i]='\0';
-j=0;
-for(i=79;i<87;i++){temp[j]=buff[i];j++;}
-chi=atof(temp);
-
-for(i=0;i<20;i++)temp[i]='\0';
-j=0;
-for(i=88;i<93;i++){temp[j]=buff[i];j++;}
-n=atof(temp);
-
-
-for(i=0;i<20;i++)temp[i]='\0';
-j=0;
-for(i=163;i<175;i++){temp[j]=buff[i];j++;}
-f1=atof(temp);
-
-for(i=0;i<20;i++)temp[i]='\0';
-j=0;
-for(i=177;i<188;i++){temp[j]=buff[i];j++;}
-f2=atof(temp);
-
/*MAKE A LOCUS LIST*/
for(i=0;i<20;i++)locs[i]='\0';
for(i=0;i<20;i++) { if(loc1[i]=='.')goto gogo; locs[i]=loc1[i]; }
@@ -146,9 +95,6 @@ abcd_p1->ldu2=0.;
abcd_p1->kb1=kb1;
abcd_p1->kb2=kb2;
abcd_p1->kb=fabs(kb2-kb1);
-abcd_p1->chi=chi;
-abcd_p1->freq1=f1;
-abcd_p1->freq2=f2;
abcd_p1->flag=1;
}
else
@@ -165,17 +111,12 @@ abcd_p1->ldu2=0.;
abcd_p1->kb1=kb1;
abcd_p1->kb2=kb2;
abcd_p1->kb=fabs(kb2-kb1);
-abcd_p1->chi=chi;
-abcd_p1->freq1=f1;
-abcd_p1->freq2=f2;
abcd_p1->flag=1;
}
skip:i=1;
} /*end while*/
-/*locus_list();
-*/
fclose(fp);
fflush(output_f);
}
diff --git a/intfileseg.c b/intfileseg.c
index fefe93c1e60a80f9d7ab4a2f313cf12f7d181be3..c5de4d503d315438ca075644b04694b4cb31f40c 100755
--- a/intfileseg.c
+++ b/intfileseg.c
@@ -12,12 +12,12 @@ char buff[200];
FILE *fp;
abcdPtr abcd_p1;
char locs[20],loc1[20],loc2[20];
-double f1,f2,chi,n,kb1,kb2,p,k;
+double chi,n,kb1,kb2,p,k;
g_nloci=0;
abcd_p1=NULL;
abcdstartPtr=NULL;
-if((fp=fopen(datfile,"r"))==NULL){ printf("\nCannot open intermediate file for reading");exit(1);}
+if((fp=fopen(intefile,"r"))==NULL){ printf("\nCannot open intermediate file for reading");exit(1);}
while(fgets(buff,200,fp)!=NULL)
@@ -61,29 +61,6 @@ j=0;
for(i=68;i<78;i++){temp[j]=buff[i];j++;}
k=atof(temp);
-for(i=0;i<20;i++)temp[i]='\0';
-j=0;
-for(i=79;i<87;i++){temp[j]=buff[i];j++;}
-chi=atof(temp);
-
-for(i=0;i<20;i++)temp[i]='\0';
-j=0;
-for(i=88;i<93;i++){temp[j]=buff[i];j++;}
-n=atof(temp);
-
-
-for(i=0;i<20;i++)temp[i]='\0';
-j=0;
-for(i=163;i<175;i++){temp[j]=buff[i];j++;}
-f1=atof(temp);
-
-for(i=0;i<20;i++)temp[i]='\0';
-j=0;
-for(i=177;i<188;i++){temp[j]=buff[i];j++;}
-f2=atof(temp);
-
-
-
/*************************************************************************/
/*************************************************************************/
/*TEST TO SEE IF WE WANT THIS PAIR*/
@@ -128,9 +105,6 @@ abcd_p1->ldu2=0.;
abcd_p1->kb1=kb1;
abcd_p1->kb2=kb2;
abcd_p1->kb=fabs(kb2-kb1);
-abcd_p1->chi=chi;
-abcd_p1->freq1=f1;
-abcd_p1->freq2=f2;
abcd_p1->flag=1;
}
else
@@ -147,9 +121,6 @@ abcd_p1->ldu2=0.;
abcd_p1->kb1=kb1;
abcd_p1->kb2=kb2;
abcd_p1->kb=fabs(kb2-kb1);
-abcd_p1->chi=chi;
-abcd_p1->freq1=f1;
-abcd_p1->freq2=f2;
abcd_p1->flag=1;
}
diff --git a/ldmap b/ldmap
index 2f11d21ae33c62185b558a593fe50b8932b6a6ed..b747ae2df8d107f1ced854c312c46412614eaa00 100755
--- a/ldmap
+++ b/ldmap
@@ -4,166 +4,58 @@ clear
echo " "
echo " "
echo " L D M A P "
-echo "Construction of linkage disequilibrium maps from haplotype and diplotype data..."
-echo " VERSION 1.0 12/11/04 "
+echo "Construction of linkage disequilibrium maps from diplotype data..."
+echo " VERSION 2.0 November 2014 "
MODE=0
-while [ $MODE -ge -1 -a $MODE -lt 9 ]
+while [ $MODE -ge -1 -a $MODE -lt 2 ]
do
echo " "
echo " "
echo ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"
echo "Enter number to select function "
-echo " 1 Haplotypes with SNPs - construct intermediate file "
+echo " 1 Construct LD map in segments from intermediate file "
-echo " 2 Diplotypes with SNPs - construct intermediate file "
-
-echo " 3 Construct LD map in segments from intermediate file "
-
-echo " 4 Construct LD map from intermediate file "
-
-echo " 5 Add LDU locations to a genotype file "
-
-echo " 6 Construct intermediate file from pooled haplotype counts "
-
-echo " 7 Test HWE in genotype file and obtain allele frequencies "
-
-echo " 8 Combine two intermediate files "
-
-echo " 9 quit "
+echo " 2 quit "
echo " "
read MODE
-while [ $MODE -ge 1 -a $MODE -lt 9 ]
+while [ $MODE -ge 1 -a $MODE -lt 2 ]
do
- if [ $MODE -eq 1 -o $MODE -eq 2 ]
+ if [ $MODE -eq 1 ]
then
- echo "*Define input and output files*"
+ echo "*Construct an LDU map*"
echo " "
echo "Input data file is: "
read DATA
- echo "Intermediate output file is: "
- read TERSE
- echo " CP file (to run the program) is: "
- read CP
- rm -f $CP
- echo "Writing an executable 'cp' file........"
- echo "ldmapper1" $DATA $TERSE $MODE >> $CP
- chmod +x $CP
- echo "Run this job (y/n) ?"
- read run
- case $run in y) /bin/sh $CP ;;
- esac
- MODE=0
- fi
-
- if [ $MODE -eq 3 ]
- then
- TMODE=MODE
- echo "Construct an LDU map in overlapping segments "
+ echo "Intermediate file is: "
+ read INT
echo "Job file is: "
read JOB
- echo "Intermediate (input) file is: "
- read DATA
- echo "Terse output file is: "
- read TERSE
- echo "Output file is: "
- read PROLIX
-
-# echo "Accept defaults throughout ?"
-# read yorn
-
-# case $yorn in
-# y)
-# TMODE=100
-#echo " CP file (to run the program) is: "
-#read CP
-#rm -f $CP
-# echo "Writing an executable 'cp' file........"
-# echo "ldmapper1" $JOB $DATA $TERSE $PROLIX $TMODE >> $CP
-# chmod +x $CP
-# ;;
-
-# n)
- ldmapper1 $JOB $DATA $TERSE $PROLIX $MODE
-# ;;
-
-# esac
- MODE=0
- fi
-
- if [ $MODE -eq 4 ]
- then
- echo "Construct an LDU map "
- echo "Define input and output files"
- echo " "
- echo "Job file is: "
- read JOB
- echo "Intermediate (input) file is: "
- read DATA
echo "Map (output) file is: "
- read TERSE
- echo "Analysis (output) file is: "
- read PROLIX
- ldmapper1 $JOB $DATA $TERSE $PROLIX $MODE
- MODE=0
- fi
-
- if [ $MODE -eq 5 ]
- then
- echo "Add LDU locations to a genotype data file "
- ldmapper1 $DATA $MODE
- MODE=0
- fi
-
- if [ $MODE -eq 6 ]
- then
- echo "Construct intermediate file from pooled counts "
- echo "Pooled count file name is: "
- read DATA
- ldmapper1 $DATA $MODE
- MODE=0
- fi
-
- if [ $MODE -eq 7 ]
- then
- echo "Testing HWE for genotypes"
- echo " "
- echo "Genotype file is "
- read FAM
- echo "Prolix file is "
+ read TERSE
+ echo "Long output file is: "
read PROLIX
- echo "CP file name "
+ echo "MAF cut-off (eg. 0.05) is: "
+ read MAF
+ echo "HWE P-value cut-off (eg. 0.001) is: "
+ read HWE
+ echo "Script file (to run the program) is: "
read CP
rm -f $CP
- echo "Writing an executable 'cp' file........"
- echo " "
- echo "betasnp1" $FAM $PROLIX 11 >> $CP
- chmod +x $CP
+ echo "Writing an executable script file........"
+ echo "ldmapper1" $DATA $INT $JOB $TERSE $PROLIX $MAF $HWE >> $CP
+ chmod +x $CP
echo "Run this job (y/n) ?"
read run
- case $run in
- y) /bin/sh $CP ;;
+ case $run in y) /bin/sh $CP ;;
esac
MODE=0
fi
- if [ $MODE -eq 8 ]
- then
- echo "Combine two intermediate files "
- echo "Define input and output files"
- echo " "
- echo "Concatenated intermediate (input) file is: "
- read DATA
- echo "Combined intermediate (output) file is: "
- read PROLIX
- ldmapper1 $DATA $PROLIX $MODE
- MODE=0
- fi
-
done
done
diff --git a/ldmap_manual.pdf b/ldmap_manual.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..843b1cb07c6c818cc1e9961284a05c3a4be3a2a5
Binary files /dev/null and b/ldmap_manual.pdf differ
diff --git a/ldmap_to_csv.pl b/ldmap_to_csv.pl
new file mode 100644
index 0000000000000000000000000000000000000000..ab7914d426ec6ede2f92e98cc753f271a7ab6eb4
--- /dev/null
+++ b/ldmap_to_csv.pl
@@ -0,0 +1,16 @@
+use warnings; use strict;
+
+## Script to convert LDmap output to CSV
+## RJPengelly, 2014, r.j.pengelly@soton.ac.uk
+
+open FILE, $ARGV[0];
+open OUT, ">$ARGV[0]" . '.csv';
+
+print OUT 'N,Locus,KBmap,LDUmap';
+while (<FILE>) {
+ if ($_ =~ /#/) {next;}
+ $_ =~ s/^\s+//;
+ $_ =~ s/\s+$//;
+ $_ =~ s/\s+/\,/g;
+ print OUT "$_\n";
+}
diff --git a/ldmapper.c b/ldmapper.c
index 7bb7351a1bc4cb5e9b8c2b0015c3567c86583569..90726158654b0490ebee1608955f475398414310 100755
--- a/ldmapper.c
+++ b/ldmapper.c
@@ -1,131 +1,41 @@
+
#include "allass.h"
int main(int argc, char *argv[])
{
-char datfile2[50],temp[20],outputfile[20],jobfile[20],flagg[5];
-/* time_t start; */
+char outputfile[20],jobfile[20], temp1[20],temp2[20];
+/*This is the maximum number of intervals used:*/
g_int=100;
+/*This is the maximum distance between pairs of SNPs used:*/
g_max=500.;
+/*This is the Hardy-Wienberg cut off */
+g_hwp=0.001;
+/*This is the MAF cut off */
+g_maf=0.05;
+
g_finish=0;
g_calls=0;
-ldflag2=0;
-if(argc==2)strcpy(flagg,argv[1]);
-if(argc==3)strcpy(flagg,argv[2]);
-if(argc==4)strcpy(flagg,argv[3]);
-if(argc==5)strcpy(flagg,argv[4]);
-if(argc==6)strcpy(flagg,argv[5]);
-if(argc==7)strcpy(flagg,argv[6]);
-g_opt=atoi(flagg);
+strcpy(datfile,argv[1]);strcpy(intefile,argv[2]); strcpy(jobfile,argv[3]);strcpy(terfile,argv[4]);strcpy(outputfile,argv[5]);
+strcpy(temp1,argv[6]);strcpy(temp2,argv[7]);
-if(g_opt==1||g_opt==2) {strcpy(datfile,argv[1]);strcpy(intefile,argv[2]);}
+g_maf=atof(temp1);
+g_hwp=atof(temp2);
-if(g_opt!=1&&g_opt!=2&&g_opt!=3)
-{
-if(argc==3){strcpy(datfile,argv[1]);}
-if(argc==4){strcpy(intefile,argv[1]);strcpy(outputfile,argv[2]);}
-if(argc==5){strcpy(jobfile,argv[1]);strcpy(intefile,argv[2]);strcpy(outputfile,argv[3]);}
-if(argc==6){strcpy(jobfile,argv[1]);strcpy(datfile,argv[2]);strcpy(intefile,argv[3]);strcpy(outputfile,argv[4]);}
-}
+if ((output_f=fopen(outputfile, "w"))==NULL){ printf("\nCannot open outputfile.\n\n"); exit(1);}
/*////////////////////////////////////////////////////////////////////////////*/
-/***OPTION 1 AND 2 ******/
-
-/*Make int files ******/
-if(g_opt==1||g_opt==2)
- {
-
+/*Make int file ******/
printf("\n**NOTE**: The current settings of the program can accomodate up to %d loci and files up to %d characters wide ",MAX_LOCI,MAX_LINE);
-strcpy(ped_file,datfile);
-strcpy(datfile2,"temp.dat");
-if((output_f3=fopen(datfile2,"w"))==NULL){printf("\nCannot open temp.dat file");exit(1);}
-fprintf(output_f3,"\n\n\nASSOCIATION IN SNP DIPLOTYPES (LDMAP)\n\n\n ");
-fflush(output_f3);
-
-
-if(g_opt==1)
-{
-printf("\nAre alleles coded as 0,1 or 1,2 ? (if 0,1 enter 0, if 1,2 enter 1)");
-scanf("%s",temp);
-if(temp[0]=='1')g_code=0;
-if(temp[0]=='0')g_code=1;
-}
-file_inputg();
-
-if(g_opt==1) { if(g_code==0){getfreqsh2();hapallele();} if(g_code==1){getfreqsh();hapallele2();} }
-
-if(g_opt==2){getfreqs();multallele();}
-fflush(output_f3);
-fclose(output_f3);
-file_input2(datfile2);
-if(g_opt==1)haplo(intefile);
-if(g_opt==2)diplo(intefile);
-}
+file_inputnew();
+getfreqs();
+multallele();
+diplo(intefile);
/*////////////////////////////////////////////////////////////////////////////*/
-if(g_opt==3)
-{
/*CONSTRUCT MAP IN SECTIONS */
-strcpy(jobfile,argv[1]);
-strcpy(g_jobfile,jobfile);
if ((job_fp=fopen(jobfile, "r"))==NULL){ printf("\nCannot open jobfile.\n\n"); exit(1);}
-
-strcpy(datfile,argv[2]);
-strcpy(intefile,argv[3]);
-strcpy(outputfile,argv[4]);
- if ((output_f=fopen(outputfile, "w"))==NULL){ printf("\nCannot open outputfile.\n\n"); exit(1);}
segments();
-}
-/*////////////////////////////////////////////////////////////////////////////*/
-/***OPTION 4 ******/
-/*Construct LD map ******/
-if(g_opt==4)
- {
- strcpy(g_jobfile,jobfile);
- if ((job_fp=fopen(jobfile, "r"))==NULL){ printf("\nCannot open jobfile.\n\n"); exit(1);}
- if ((output_f=fopen(outputfile, "w"))==NULL){ printf("\nCannot open outputfile.\n\n"); exit(1);}
- printf("\n\n\nL D M A P - VERSION 1.0 CREATED 12/11/04 \n\n\n");
- fprintf(output_f,"\n\n\nL D M A P - VERSION 1.0 CREATED 12/11/04 \n\n\n");
- intfile3();
- locus_list();
- ldmap2();
- }
-
/*////////////////////////////////////////////////////////////////////////////*/
-/***OPTION 5 ******/
-
-if(g_opt==5)mergeloc();
-/*////////////////////////////////////////////////////////////////////////////*/
-/*////////////////////////////////////////////////////////////////////////////*/
-/***OPTION 6 ******/
-/*Process pooled count file (hill.out) ******/
-if(g_opt==6)
-{
-printf("\nEnter name of intermediate file to be constructed :");
-scanf("%s",intefile);
-if ((output_f=fopen(intefile, "w"))==NULL){ printf("\nCannot open intermed file.\n\n"); exit(1);}
- pool(datfile);
-}
-
-if(g_opt==8)
-{
-strcpy(datfile,outputfile);
-intfile7();
-
-}
-/*////////////////////////////////////////////////////////////////////////////*/
-/*////////////////////////////////////////////////////////////////////////////*/
-/*////////////////////////////////////////////////////////////////////////////*/
-if(g_opt==9)
-{
-/*GENERATE A SUB-SAMPLE DROPPING ONE INDIVIDUAL */
-strcpy(ped_file,argv[1]);strcpy(outputfile,argv[2]);
-if ((output_f=fopen(outputfile, "w"))==NULL){ printf("\nCannot open outputfile.\n\n"); exit(1);}
-file_inputg2();
-fclose(output_fi);
-}
-
-
-
return 0;
}
diff --git a/ldmapper1 b/ldmapper1
index 89f65b16801122e823e6853f3b92eac0380dcac0..beccfbe56541c26c1a0ed617f6acfc5e6ada2cf7 100755
Binary files a/ldmapper1 and b/ldmapper1 differ
diff --git a/ldmapseg.c b/ldmapseg.c
index de6405e4fb611e3bc3c551763f50c88c6ad9d982..2f7683a14e68634b1ace48c096ded0220d3f1d36 100755
--- a/ldmapseg.c
+++ b/ldmapseg.c
@@ -14,22 +14,17 @@ abcdPtr abcd_p1;
hugsPtr hug1;
intsPtr intop;
int ntimes=0,maxit,ix,i,jth,j,ii;
-/* set the vv variables..*/
-/* double vved_old,vved_new; */
double vvthee,vved,vvdd;
int vvjj,vvleft,vvright;
intsPtr vvintop;
hugsPtr vvhug1;
-/** aman added 0926 ***/
-
g_estimatel=0;
g_initial=0;
conv=0.1;
minkb=large;
g_best=large;
maxit=20000;
-ldflag2=1;
itlast=0; iteast=0; itmast=0;
useld=0;
intervals=g_int;
@@ -39,11 +34,10 @@ stope=1;
if(path==0)fprintf(output_f,"\n\n\n******** Fitting pairwise data to kb map ******** ");
if(path==1)fprintf(output_f,"\n\n\n******** Creating an LD map from a Kb map ******** ");
if(path==2)fprintf(output_f,"\n\n\n******** Fitting pairwise data to LD map ******** ");
-if(path==3)fprintf(output_f,"\n\n\n******** Creating an LD map from an existing LD map ******** ");
/******************************************************************************************************/
start1=time(NULL);
if(path==2)istop=1;
-if((path==2)||(path==3)) { useld=1; directmap2();}
+if((path==2)) { useld=1; directmap2();}
/*IF NEEDED INSERT EXISTING LDU MAP AT THIS POINT*/
if(useld==1)
@@ -101,9 +95,6 @@ orig_e=global_e;
swept=1./orig_e;
/****************************************************************************************************/
-/*
-reorder();
-*/
/****************************************************************************************************/
maxkb=maxkb-minkb;
@@ -119,7 +110,7 @@ fprintf(output_f,"\n\nN(number of pairs)=%15d m(number of SNPs)=%15d df=%14.1f
fflush(output_f);
if(path==0)return;
-if(istop==1) { /*printmap();*/ goto fin; }
+if(istop==1) { goto fin; }
/****************************************************************************************************/
jth=0;
@@ -154,8 +145,6 @@ for(ii=0;ii<g_nloci;ii++)
}
}
-if(path==3)printinputmap();
-
fprintf(output_f,"\n\n\n***** CONSTRUCTING LDU MAP.................\n ");
if(path==1)
{
@@ -166,17 +155,6 @@ fprintf(output_f,"\n %10.6f %10.6f %10.6f %16.8f %12.5f ",
global_e,global_l,global_m,maxkb*global_e,g_lnl);
}
-if(path==3)
-{
-if(g_estimatel==0)global_l=g_pred;
-quicklike(global_e,global_l,global_m);
-
-fprintf(output_f,"\n\n\nIter E L M n_holes LDU_length -2lnlk");
-fprintf(output_f,
-"\n %10.6f %10.6f %10.6f %12.5f ",
-global_e,global_l,global_m,g_lnl);
-}
-
g_startlnl=g_lnl;
fflush(output_f);
@@ -218,7 +196,6 @@ hug1=bigstart;
ii=0;
allldu=0;
-/********** aman added ...computing the first value of ed 0925 *************/
/*THIS IS NOT LAYERED ANY MORE */
vvhug1=bigstart;
diff --git a/locus_list.c b/locus_list.c
index b04e55eda845fb32172fe758c5923499fd8cb15e..dbf015f47ae2c1cef6eba20ff0af21a785a49357 100755
--- a/locus_list.c
+++ b/locus_list.c
@@ -5,7 +5,7 @@
void locus_list()
{
int i,ii;
-double kb1,kb2,ldu1,ldu2,freq1,freq2;
+double kb1,kb2,ldu1,ldu2;
abcdPtr abcd_p1;
char locs[20];
/*************GET THE LOCATIONS FOR ALL IN THE LOCUS LIST ***/
@@ -21,19 +21,17 @@ kb1=0;
kb2=0;
ldu1=0;
ldu2=0;
-freq1=0;
-freq2=0;
abcd_p1=abcdstartPtr;
while(abcd_p1!=NULL)
{
for(i=0;i<20;i++)locs[i]='\0';
for(i=0;i<20;i++) {
-if(abcd_p1->locus[i]=='.')goto gogo4; strcpy(locs,abcd_p1->locus);kb1=abcd_p1->kb1;ldu1=abcd_p1->ldu1;freq1=abcd_p1->freq1; }
+if(abcd_p1->locus[i]=='.')goto gogo4; strcpy(locs,abcd_p1->locus);kb1=abcd_p1->kb1;ldu1=abcd_p1->ldu1;}
gogo4: for(ii=0;ii<g_nloci;ii++)
{
-if(strcmp(g_loci[ii],locs)==0){g_location[ii][0]=ldu1;g_location[ii][1]=kb1;g_location[ii][2]=freq1; }
+if(strcmp(g_loci[ii],locs)==0){g_location[ii][0]=ldu1;g_location[ii][1]=kb1;}
}
abcd_p1=abcd_p1->nextPtr;
}
@@ -42,10 +40,10 @@ abcd_p1=abcdstartPtr;
while(abcd_p1!=NULL)
{
for(i=0;i<20;i++)locs[i]='\0';
-for(i=0;i<20;i++) { if(abcd_p1->locus2[i]=='.')goto gogo3; strcpy(locs,abcd_p1->locus2);kb2=abcd_p1->kb2;ldu2=abcd_p1->ldu2;freq2=abcd_p1->freq2; }
+for(i=0;i<20;i++) { if(abcd_p1->locus2[i]=='.')goto gogo3; strcpy(locs,abcd_p1->locus2);kb2=abcd_p1->kb2;ldu2=abcd_p1->ldu2;}
gogo3: for(ii=0;ii<g_nloci;ii++)
{
-if(strcmp(g_loci[ii],locs)==0){g_location[ii][0]=ldu2;g_location[ii][1]=kb2;g_location[ii][2]=freq2;}
+if(strcmp(g_loci[ii],locs)==0){g_location[ii][0]=ldu2;g_location[ii][1]=kb2;}
}
abcd_p1=abcd_p1->nextPtr;
}
diff --git a/makefile b/makefile
index 16d6c2042d25b647ae497552bb66080f1315ffb0..1830f2a6a38a6a386d2f35734f273e0f213c845f 100755
--- a/makefile
+++ b/makefile
@@ -1,12 +1,12 @@
# @(#)makefile 1.15
SOURCES= \
- simboot.o file_inputg2.o ran1.o intfile7.o getfreqsh2.o getfreqsh.o holes.o intfile6.o directmap2.o updat1.o writeterfinseg.o keep_it2.o ldmapseg.o reorderall.o intfile4.o reorder.o writeterfin.o writeter.o printinputmap.o getfreqs.o intfileseg.o segments.o quicklike.o hapallele2.o diplo1.o multallele1.o haplo.o hapallele.o printmap.o locus_list.o fast.o directmap.o one_interval.o pool.o ludcmp.o keep_it.o update.o final.o get_predl.o nrutil.o ldmap3.o fill_str2.o mergeloc.o meankchi.o multallele.o file_inputg.o runewt3.o ldmap2.o bound.o intfile3.o metricho.o diplo.o file_input2.o func.o dfunc.o dfpmin.o lnsrch.o cal_SE2.o getpa.o jobin.o runewt2.o ldmapper.o
+ chi.o ran1.o newentry2.o newentry.o file_inputnew.o fill_gai.o intfile6.o directmap2.o updat1.o writeterfinseg.o keep_it2.o ldmapseg.o reorderall.o intfile4.o reorder.o getfreqs.o intfileseg.o segments.o quicklike.o locus_list.o fast.o one_interval.o ludcmp.o update.o get_predl.o nrutil.o meankchi.o multallele.o runewt3.o bound.o metricho.o diplo.o func.o dfunc.o dfpmin.o lnsrch.o cal_SE2.o getpa.o jobin.o runewt2.o ldmapper.o
PROGRAM= ldmapper1
# gcc can take "-Wall -Werror" for checking
-CFLAGS= -g -O
+CFLAGS= -g -O -Wall -Werror
CPPFLAGS= $(INCLUDEDIRS)
LDFLAGS=
LINTFLAGS=
diff --git a/meankchi.c b/meankchi.c
index 930296e653963f48f74829bb4633e8339f95aa81..aa59f902476f8d39721a55b750c9df5758b93ca0 100755
--- a/meankchi.c
+++ b/meankchi.c
@@ -6,8 +6,6 @@
/*GET LOCUS SPECIFIC CHI SQUARE ETC */
double lnl=0., xx,d,k,n=0,po,pe,E,L,M;
- g_meanchi=0;
- g_meank=0;
g_abcd=abcdstartPtr;
while(g_abcd!=NULL)
{
@@ -27,7 +25,4 @@
}
g_abcd=g_abcd->nextPtr;
}
- g_meanchi=-2.*lnl;
- g_meank=n;
-
}
diff --git a/metricho.c b/metricho.c
index 67a6fdeaf9827bf7469e310ac21a1ac567807869..e829b8c0dcc8bee2c3530cf6240a68e11ab21542 100755
--- a/metricho.c
+++ b/metricho.c
@@ -12,12 +12,6 @@ D=DDD;
kd=n/(Q*(1.-Q)*R*(1.-R));
/*The metrics and their information */
-
-/*Covariance */
-C=1.;
-g_cov=D/C;
-g_covi=(C*C)*kd;
-
/*Association rho */
C=(Q*(1.-R));
g_rho=D/C;
diff --git a/multallele.c b/multallele.c
index ccc42a9bcc01ddac8539f2de2a47d69fe9366e32..0e59eb43edddf22ebadcdde70e359a01cf502d77 100755
--- a/multallele.c
+++ b/multallele.c
@@ -4,9 +4,8 @@ void multallele()
{
/*SET UP 9 CELL TABLES FOR DIPLOTYPE DATA (ASSOCIATION) -FOR SNPS ONLY */
/*SNPS MUST BE CODED 1,2 */
-double diff,kb,kb1,kb2,sum,maxkb;
-char temp[20];
-int theint,i1,i2,idiff,N,j,numg,i,n11,n12,n13;
+double diff,kb1,kb2,sum,maxkb;
+int theint,i1,i2,idiff,N,j,i,n11,n12,n13;
int n21,n22,n23;
int n31,n32,n33;
pedPtr ped1;
@@ -14,6 +13,7 @@ mapPtr2 map1,map2;
maxkb=g_max;
+/*
printf("\n\nThe maximum distance between a pair of SNPs to be used is: %10.4f kb, do you wish to change this? (y/n) ",maxkb);
scanf("%s",temp);
if(temp[0]=='y')
@@ -23,8 +23,9 @@ scanf("%s",temp);
maxkb=atof(temp);
g_max=maxkb;
}
-
+*/
theint=g_int;
+/*
printf("\nThe maximum number of intervals between pairs of SNPs to be used is: %d, do you wish to change this? (y/n) ",theint);
scanf("%s",temp);
if(temp[0]=='y')
@@ -34,6 +35,7 @@ scanf("%s",temp);
theint=atoi(temp);
g_int=theint;
}
+*/
/*******************/
/*******************/
i=0;
@@ -46,34 +48,6 @@ map1=map1->nextPtr;
}
/*******************/
/*******************/
-
-
-ped1=ped_startPtr;
-numg=ped1->numg;
-
-/************************************************************************/
-/************************************************************************/
-
-fprintf(output_f3,"\n\n\nNumber of loci =%d ",numg+1);
-fprintf(output_f3,"\n------------------------------------------------------------------------------------------------------------------------");
-
-fprintf(output_f3,"\n(1111,1-4) ");
-fprintf(output_f3,"\n(1112,6-9) ");
-fprintf(output_f3,"\n(1122,11-14) ");
-fprintf(output_f3,"\n(1211,16-19) ");
-fprintf(output_f3,"\n(1212,21-24) ");
-fprintf(output_f3,"\n(1222,26-29) ");
-fprintf(output_f3,"\n(2211,31-34) ");
-fprintf(output_f3,"\n(2221,36-39) ");
-fprintf(output_f3,"\n(2222,41-44) ");
-fprintf(output_f3,"\n(N,46-49) ");
-fprintf(output_f3,"\n(locus1,52-63) ");
-fprintf(output_f3,"\n(locus2,69-80) ");
-fprintf(output_f3,"\n(kb1,86-97) ");
-fprintf(output_f3,"\n(kb2,103-114) ");
-fprintf(output_f3,"\n(freq1,120-131) ");
-fprintf(output_f3,"\n(freq2,134-145) ");
-fprintf(output_f3,"\n------------------------------------------------------------------------------------------------------------------------");
/************************************/
i=0;
map1=map_startPtr2;
@@ -84,39 +58,22 @@ j=i+1;
map2=map1->nextPtr;
while(map2!=NULL)
{
-
-
-
-kb1=map1->kb;
-kb2=map2->kb;
-i2=map2->order;
-
-diff=fabs(kb1-kb2);
-if(diff>maxkb)goto skip;
-
-idiff=abs(i1-i2);
-
-
-if(idiff>theint+1)goto skip;
+ kb1=map1->kb;
+ kb2=map2->kb;
+ i2=map2->order;
+ diff=fabs(kb1-kb2);
+ if(diff>maxkb)goto skip;
+ idiff=abs(i1-i2);
+ if(idiff>theint+1)goto skip;
/************************************************************************/
/*THIS SECTION REFERS TO A PAIR OF LOCI **/
/*FOR THIS PAIR HAVE TO LOOK AT EACH ALLELE COMBINATION */
/************************************************************************/
-n11=0;
-n12=0;
-n13=0;
-n21=0;
-n22=0;
-n23=0;
-n31=0;
-n32=0;
-n33=0;
+n11=0; n12=0; n13=0; n21=0; n22=0; n23=0; n31=0; n32=0; n33=0;
ped1=ped_startPtr;
while(ped1!=NULL)
{
- if(ped1->N==999)ped1->N=1;
- if(ped1->N==0)ped1->N=1;
- N=ped1->N;
+ N=1;
if(ped1->GEN[i][0]!=0&&ped1->GEN[i][1]!=0&&ped1->GEN[j][0]!=0&&ped1->GEN[j][1]!=0)
{
if((ped1->GEN[i][0]==1)&&(ped1->GEN[i][1]==1))
@@ -159,17 +116,19 @@ sum=n11+n12+n13+n21+n22+n23+n31+n32+n33;
kb1=map1->kb;
kb2=map2->kb;
-kb=fabs(kb1-kb2);
-map1->locus[13]='\0';
-map2->locus[13]='\0';
if(sum>0.)
{
-if(map1->freq>0.05&&map2->freq>0.05)
+/*Select pairs of loci with sample MAF >0.05 for both*/
+
+if(map1->chi>g_hwp&&map2->chi>g_hwp)
+{
+if(map1->freq>g_maf&&map2->freq>g_maf)
{
-fprintf(output_f3,"\n%4d %4d %4d %4d %4d %4d %4d %4d %4d ",n11,n12,n13,n21,n22,n23,n31,n32,n33);
-fprintf(output_f3,"%4d %12s %12s %12.4f %12.4f %12.10f %12.10f",
-n11+n12+n13+n21+n22+n23+ n31+n32+n33,map1->locus,map2->locus,kb1,kb2,map1->freq,map2->freq);
+fill_gai(n11,n12,n13,n21,n22,n23,n31,n32,n33,kb1,kb2,map1->locus, map2->locus, map1->freq, map2->freq);
}
+}
+
+
}
skip:map2=map2->nextPtr;
j++;
@@ -178,6 +137,4 @@ i++;
map1=map1->nextPtr;
}
-
-
}
diff --git a/newentry.c b/newentry.c
new file mode 100644
index 0000000000000000000000000000000000000000..85e3c6b825fae6f3bf664e1a0c191ac0ab105890
--- /dev/null
+++ b/newentry.c
@@ -0,0 +1,20 @@
+#include "allass.h"
+
+void newentry()
+{
+int i;
+/******************************************************************************************/
+ if(g_ped1==NULL)
+ {
+ ped_startPtr=(pedPtr)malloc(sizeof(peds));
+ ped_startPtr->nextPtr=NULL;
+ g_ped1=ped_startPtr;
+ for(i=0;i<MAX_LOCI;i++){g_ped1->GEN[i][0]=0;g_ped1->GEN[i][1]=0;}
+ }
+ else
+ {
+ g_ped1->nextPtr=(pedPtr)malloc(sizeof(peds));
+ g_ped1=g_ped1->nextPtr;
+ for(i=0;i<MAX_LOCI;i++){g_ped1->GEN[i][0]=0;g_ped1->GEN[i][1]=0;}
+ }
+}
diff --git a/newentry2.c b/newentry2.c
new file mode 100644
index 0000000000000000000000000000000000000000..ce382a6d2c2a31758da2bec21fb620a103981441
--- /dev/null
+++ b/newentry2.c
@@ -0,0 +1,17 @@
+#include "allass.h"
+
+void newentry2()
+{
+/*Make new entry in map structure */
+ if(g_p12==NULL)
+ {
+ map_startPtr2=(mapPtr2)malloc(sizeof(mapkb2));
+ map_startPtr2->nextPtr=NULL;
+ g_p12=map_startPtr2;
+ }
+ else
+ {
+ g_p12->nextPtr=(mapPtr2)malloc(sizeof(mapkb2));
+ g_p12=g_p12->nextPtr;
+ }
+}
diff --git a/one_interval.c b/one_interval.c
index c141405e40066f0159730ec008fe3e8daef5df17..0e4d3c8fe6df8a5945d24ff8789ed21a208fc54a 100755
--- a/one_interval.c
+++ b/one_interval.c
@@ -5,9 +5,7 @@ void one_interval(int ii, intsPtr intop,double *u, double *bigk, double *lnl )
/*INPUT Malecot M and L, INPUT AND OUTPUT u,bigk,lnl*/
/*GLOBAL - NEEDS TO BE PASSED = interv ARRAY + global_l, global_m, g_nloci*/
-/* double thee,u1,bigk1,lnl1=0.,term0,term,l,m,ed,dd,pexp,k,p; */
double u1,bigk1,lnl1=0.,term0,term,l,m,ed,dd,pexp,k,p;
-/*int jj,left,right; */
u1=0.;
bigk1=0.;
diff --git a/runewt2.c b/runewt2.c
index a11a52992ef37205c493a33e7076be4717297107..5f370493f6fcdf88d15551bd3e1116ef0d893f3f 100755
--- a/runewt2.c
+++ b/runewt2.c
@@ -20,7 +20,7 @@ fprintf(output_f,
"\nIT -2lnlk E L M Kee Kll Kmm ");
fflush(output_f);
g_lnl = 999999.;
-g_kel=0.; g_kem=0.; g_ket=0.; g_kes=0.; g_klm=0.; g_klt=0.; g_kls=0.; g_kms=0.; g_kts=0.;
+g_kel=0.; g_kem=0.; g_klm=0.;
n=0;
global_fin=0;
g_iter[1]=0;
diff --git a/runewt3.c b/runewt3.c
index 6dc519a7a0b883c4bdcb23b05d9a2a5aad392207..00632abbfa1031ad4e5a4bb9a70c547645b1bc89 100755
--- a/runewt3.c
+++ b/runewt3.c
@@ -13,7 +13,7 @@ fprintf(output_f,"\n------------------------------------------------------------
fprintf(output_f,"\nIT -2lnlk E L M S Kee Kll Kmm Kss ");
fflush(output_f);
}
-g_kel=0.; g_kem=0.; g_ket=0.; g_kes=0.; g_klm=0.; g_klt=0.; g_kls=0.; g_kms=0.; g_kts=0.;
+g_kel=0.; g_kem=0.; g_klm=0.;
n=0;
global_fin=0;
diff --git a/segments.c b/segments.c
index 70f50a5053d53bfa23c41d104d100b98e41b7291..7f407e4838b384ee6a659006c0258c53a5d40494 100755
--- a/segments.c
+++ b/segments.c
@@ -5,7 +5,6 @@
void segments()
{
/*READ IN A PORTION OF AN INT FILE - FOR SEGMENTAL MAP CONSTRUCTION */
-char temp[10];
int j,i;
int nseg=-1,bign;
int overlap=-1,pos,segsize,a,b,ifin=0;
@@ -16,8 +15,8 @@ the segmental construction of an ld map*/
printf("\n....Compiling locus list....");
intfile4();
-printf("\nThe input map has %d loci in total ",g_nallloci);
-fprintf(output_f,"\nThe input map has %d loci in total ",g_nallloci);
+printf("\nThe map has %d loci in total ",g_nallloci);
+fprintf(output_f,"\nThe map has %d loci in total ",g_nallloci);
/*!!!!!!!!!!!!!!!!!!!!*/
/*MODIFY FOR IRIDIS - JULY 09 */
@@ -28,14 +27,6 @@ if(nseg<1)nseg=1;
fprintf(output_f,"\nCompiling the map in %d segments ",nseg);
printf("\nCompiling the map in %d segments ",nseg);
-/*
-printf("\nFor assembling the LD map how many segments are required ?");
-scanf("%s", temp);
-nseg=atoi(temp);
-*/
-
-
-
for(i=0;i<g_nallloci;i++) { g_alllocation[i][0]=large; }
bign=g_nallloci;
@@ -45,13 +36,6 @@ segsize=segsize+1; /*ADD 1 TO MANAGE ROUNDING */
printf("\nThis has approximately %d loci per segment",segsize-1);
fprintf(output_f,"\nThis has approximately %d loci per segment",segsize-1);
-
-/*
-printf("\nWhat overlap is required ? ");
-scanf("%s", temp);
-overlap=atoi(temp);
-*/
-
overlap=25;
fprintf(output_f,"\nThe overlap is %d ",overlap);
@@ -131,8 +115,6 @@ for(i=0;i<g_nallloci;i++)
/*NOW DEAL WITH THE NON OVERLAP SECTION WHICH FOLLOWS - NEED A SECOND SMALL OFFSET TO ACCOUNT
FOR SMALL DISCREPANCIES IN SCALES IN THE OVERLAP SECTION - MAY BE NEGATIVE*/
offset=big1-big2;
-/*fprintf(output_f,"\nSECOND OFFSET = %f big1 = %f big2 = %f ",offset,big1,big2);
-*/
for(i=0;i<g_nallloci;i++)
{
diff --git a/writeterfinseg.c b/writeterfinseg.c
index 1f600735344d8defd3899619a73d3179d3a22d86..f9eb1f6602569ecb51964d84b5bb9fdfc9615704 100755
--- a/writeterfinseg.c
+++ b/writeterfinseg.c
@@ -7,7 +7,7 @@ char flank1[20],copy1[20],flank2[20],copy2[20];
int jth,ii;
/* removed */
-if((output_ft=fopen(intefile,"w"))==NULL){ printf("\nCannot open map (output) file");exit(1);}
+if((output_ft=fopen(terfile,"w"))==NULL){ printf("\nCannot open map (output) file");exit(1);}
fprintf(output_ft,"# ");
fprintf(output_ft,"\n# LDU MAP PARAMETERS..............................");
fprintf(output_ft,"\n# E= %10.6f L= %10.6f M=%10.6f -2*lnlk= %12.5f",global_e,global_l,global_m,g_lnl);