Skip to content

Commit ccaecd1

Browse files
mengyaomengyao
committed
cpoly cpulldown cascertain memory error fixed
1 parent e42174e commit ccaecd1

6 files changed

Lines changed: 76 additions & 106 deletions

File tree

‎src/cascertain.c‎

Lines changed: 15 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,7 @@
22
* cascertain.c: Pull down the SNPs that match the ascertain criterion.
33
* Author: Nick Patterson
44
* Revised by: Mengyao Zhao
5-
* Last revise date: 2015-05-07
5+
* Last revise date: 2015-08-10
66
* Contact: mengyao_zhao@hms.harvard.edu
77
*/
88

@@ -30,13 +30,10 @@ typedef struct {
3030
char *table_path = NULL;
3131
char *regname = NULL ;
3232
char *snpname = NULL ;
33-
//char *iubfile = NULL;
34-
//char *iubmaskfile = NULL;
3533
char *iubfile = "/home/mz128/cteam/dblist/hetfa_postmigration.dblist" ;
3634
char *iubmaskfile = "/home/mz128/cteam/dblist/mask_postmigration.dblist" ;
3735

3836
char *parname = NULL ;
39-
//int pagesize = 20*1000*1000 ; // page size for getiub
4037
int pagesize = 1000*1000 ; // page size for getiub
4138
int minfilterval = 1 ;
4239

@@ -480,15 +477,17 @@ int loadfa(char **poplist, int npops, FATYPE ***pfainfo, char *reg, int lopos, i
480477
sprintf(region, "%s:%d-%d", reg, lo, hipos);
481478

482479

483-
if (db == 0) refname = strcat(table_path, "Href.fa");
484-
else getdbname(iubfile, "Href", &refname);
480+
if (db == 0) {
481+
strcpy(refname, table_path);
482+
strcat(refname, "Href.fa");
483+
} else getdbname(iubfile, "Href", &refname);
485484

486485
if (ncall==1) {
487486
ZALLOC(falist, npops, char *) ;
488487
ZALLOC(famasklist, npops, char *) ;
489488
if (db == 0) {
490-
numfalist = setfalist(poplist, npops, ".fa", falist) ;
491-
t = setfalist(poplist, npops, ".filter.fa", famasklist) ;
489+
numfalist = setfalist(poplist, npops, ".ccomp.fa.rz", falist) ;
490+
t = setfalist(poplist, npops, ".ccompmask.fa.rz", famasklist) ;
492491
} else {
493492
numfalist = getfalist(poplist, npops, iubfile, falist) ; // set falist with the absolute path of hetfa files in .dblist file; falist contains the iubfile names
494493
t = getfalist(poplist, npops, iubmaskfile, famasklist) ;
@@ -508,6 +507,7 @@ int loadfa(char **poplist, int npops, FATYPE ***pfainfo, char *reg, int lopos, i
508507
hasmask[k] = NO ;
509508
continue ;
510509
}
510+
// fprintf(stderr, "famasklist[%d]: %s\n", k, famasklist[k]);
511511
t = strcmp(famasklist[k], "NULL") ;
512512
if (t==0) {
513513
hasmask[k] = NO ;
@@ -534,7 +534,6 @@ int loadfa(char **poplist, int npops, FATYPE ***pfainfo, char *reg, int lopos, i
534534
printf("loading: %s\n", fapt -> faname) ;
535535
printnl() ;
536536
}
537-
// fprintf(stderr, "faname1: %s\n", fapt->faname);
538537

539538
fapt -> fai = fai_load(fapt -> faname) ;
540539
fapt -> popnum = k ;
@@ -631,7 +630,6 @@ int getiub(char *cc, char *ccmask, FATYPE **fainfo, char *reg, int pos)
631630
char regbuff[128] ;
632631
static long ncnt = 0 ;
633632
static long ncall = 0 ;
634-
//fprintf(stderr, "!in getiub!\n");
635633
++ncall ;
636634
fapt = fainfo[0] ;
637635
lastreg = fapt -> regname ;
@@ -661,7 +659,6 @@ int getiub(char *cc, char *ccmask, FATYPE **fainfo, char *reg, int pos)
661659
if (pos < lastlo) newpage = YES ;
662660
if (pos > lasthi) newpage = YES ;
663661
if (ncall == 1) newreg = YES ;
664-
//fprintf(stderr, "half getiub\n");
665662
if (newreg == YES) {
666663
fflush(stdout) ;
667664
freestring(&regname) ;
@@ -685,7 +682,6 @@ int getiub(char *cc, char *ccmask, FATYPE **fainfo, char *reg, int pos)
685682
if (pos < fapt -> lopos) return -2 ;
686683
if (pos > fapt -> hipos) return -2 ;
687684
cc[k] = getfacc(fapt, pos, 1) ; // genotype at pos; cc[k] is an iub code
688-
//fprintf(stderr, "pos: %d\tcc[%d]: %c\n", pos, k, cc[k]);
689685
if (hasmask[k]) ccmask[k] = getfacc(fapt, pos, 2) ; // mask at pos
690686
else ccmask[k] = '9' ;
691687
}
@@ -697,7 +693,6 @@ int getiub(char *cc, char *ccmask, FATYPE **fainfo, char *reg, int pos)
697693
}
698694

699695
++ncnt ;
700-
// fprintf(stderr, "ncnt: %d\n", ncnt);
701696
if (ncnt == 1) {
702697
printf("zz pos: %s %s\n", cc, ccmask) ;
703698
for (k=0; k<npops; ++k) {
@@ -800,13 +795,15 @@ int setfalist(char **poplist, int npops, char *dbfile, char **iublist) {
800795
int t;
801796
for (t = 0; t < npops; ++t) {
802797
iublist[t] = strdup(table_path);
803-
iublist[t] = (char*) realloc(iublist[t], 64);
798+
iublist[t] = (char*) realloc(iublist[t], strlen(iublist[t]) + strlen(poplist[t]) + strlen(dbfile) + 1);
804799
iublist[t] = strcat(iublist[t], poplist[t]);
805-
if ((!strcmp (poplist[t], "Chimp") || !strcmp (poplist[t], "Href")) && strcmp (dbfile, ".fa")) {
806-
free (iublist[t]);
800+
// fprintf(stderr, "poplist[%d]: %s\ndbfile: %s\n", t, poplist[t], dbfile);
801+
802+
if ((!strcmp (poplist[t], "Chimp") || !strcmp (poplist[t], "Href")) && !strcmp (dbfile, ".ccomp.fa.rz"))
803+
iublist[t] = strcat(iublist[t], ".fa");
804+
else if ((!strcmp (poplist[t], "Chimp") || !strcmp (poplist[t], "Href")) && !strcmp (dbfile, ".ccompmask.fa.rz"))
807805
iublist[t] = "NULL";
808-
} else
809-
iublist[t] = strcat(iublist[t], dbfile);
806+
else iublist[t] = strcat(iublist[t], dbfile);
810807
}
811808
return npops;
812809
}

‎src/cmakefilter.c‎

Lines changed: 20 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
// Revised by Mengyao Zhao on 2015-03-12
1+
// Revised by Mengyao Zhao on 2015-08-06
22

33
#include <libgen.h>
44
#include <nicksam.h>
@@ -699,6 +699,7 @@ int readfa(char **falist, char **fasta, int *flen, int n)
699699
flen[k] = len ;
700700
}
701701
}
702+
702703
void *setvcff(char *vcffile, char **vcfn)
703704
// vcfffile should not exist
704705
{
@@ -711,23 +712,27 @@ void *setvcff(char *vcffile, char **vcfn)
711712
sprintf(iname, "%s.gz", vcffile) ;
712713
if (ftest(iname)) {
713714
printf("got %s\n", iname) ;
715+
bname = basename(vcffile) ;
716+
sprintf(outname, "%s/%s:%s.vcf", wkdir,sampname,regname ) ;
717+
sprintf(ss, "cp %s %s.gz", iname, outname) ;
718+
printf("zzcopy: %s\n",ss) ;
719+
system(ss) ;
720+
sprintf(ss, "gunzip %s.gz", outname) ;
721+
system(ss) ;
722+
sprintf(ss, "chmod 664 %s", outname) ;
723+
system(ss) ;
724+
if (ftest(outname) == NO) fatalx("unzip fails\n") ;
714725
}
715-
else {
716-
fatalx("gzfetch fails\n") ;
726+
else {
727+
sprintf(outname, "%s.bgz", vcffile) ;
728+
// sprintf(iname, "%s.bgz", vcffile) ;
729+
// if (ftest(iname)) printf("got %s\n", iname) ;
730+
// else fatalx("gzfetch and bgzfetch fail\n") ;
731+
// else {
732+
// fatalx("gzfetch fails\n") ;
717733
}
718-
bname = basename(vcffile) ;
719-
sprintf(outname, "%s/%s:%s.vcf", wkdir,sampname,regname ) ;
720-
sprintf(ss, "cp %s %s.gz", iname, outname) ;
721-
printf("zzcopy: %s\n",ss) ;
722-
system(ss) ;
723-
sprintf(ss, "gunzip %s.gz", outname) ;
724-
system(ss) ;
725-
sprintf(ss, "chmod 664 %s", outname) ;
726-
system(ss) ;
727-
if (ftest(outname) == NO) fatalx("unzip fails\n") ;
728734

729735
*vcfn = strdup(outname) ;
730-
731736
}
732737

733738
int readvcf(char *vcffile, int **fasta, int reflen, int *flen)
@@ -748,6 +753,7 @@ int readvcf(char *vcffile, int **fasta, int reflen, int *flen)
748753

749754
if (ftest(vcffile) == NO) {
750755
setvcff(vcffile, &vcftmp) ;
756+
fprintf(stderr, "vcf: %s\n", vcftmp);
751757
openit(vcftmp, &vcffp, "r") ; // aborts on error
752758
}
753759
else {

‎src/cpoly.c‎

Lines changed: 21 additions & 48 deletions
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,7 @@
22
* cpoly.c: This program is used to extract heterozygote SNPs from multiple samples
33
* Author: Nick Patterson
44
* Revised by: Mengyao Zhao
5-
* Last revise date: 2015-05-07
5+
* Last revise date: 2015-08-10
66
* Contact: mengyao_zhao@hms.harvard.edu
77
*/
88

@@ -26,13 +26,7 @@ typedef struct {
2626
} ASC ;
2727

2828
char *table_path = NULL;
29-
//char *iubfile = NULL ;
30-
//char *iubmaskfile = NULL ;
31-
3229
char *regname = NULL ;
33-
//char *parflist = "/home/np29/biology/neander/nickdir/xwdir/may12src/parfxlm" ;
34-
//char *iubfile = "/home/np29/cteam/release/hetfaplus.dblist" ; // default database
35-
//char *iubmaskfile = "/home/np29/cteam/release/maskplus.dblist" ;
3630
char *iubfile = "/home/mz128/cteam/dblist/hetfa_postmigration.dblist" ;
3731
char *iubmaskfile = "/home/mz128/cteam/dblist/mask_postmigration.dblist" ;
3832
char *parname = NULL ;
@@ -166,7 +160,7 @@ int main(int argc, char **argv)
166160
else xchrom = atoi(regname) ;
167161
}
168162

169-
fprintf(stderr, "call numlines in [main]\n");
163+
// fprintf(stderr, "call numlines in [main]\n");
170164
npops = numlines(indivname) ;
171165
ZALLOC(poplist, npops, char *) ;
172166
npops = getss(poplist, indivname) ;
@@ -209,7 +203,6 @@ int main(int argc, char **argv)
209203
lo = lopos ;
210204
hi = MIN(hipos, lo+pagesize) ;
211205
loadfa(poplist, npops, &fainfo, reg, lo, hi) ;
212-
//fprintf(stderr, "main: fapt->fai: %p\n", fainfo[0]->fai);
213206
printf("npops: %d\n", npops) ;
214207
for (k=0; k< npops; ++k) {
215208
fapt = fainfo[k] ;
@@ -234,7 +227,6 @@ int main(int argc, char **argv)
234227
reg = regname ;
235228

236229
for (pos = lopos ; pos <= hipos; ++pos) {
237-
//t = getiub(cc, ccmask, fainfo, reg, pos) ;
238230
t = getiub(cc, ccmask, fainfo, ss, pos) ;
239231
if (t==-5) break ;
240232
if (t<0) continue ;
@@ -307,9 +299,7 @@ void prints(FILE *fff, int pos, char c1, char c2)
307299
fprintf(fff, "%15s ", sss) ;
308300
fprintf(fff, "%3s 0 %12d ", regname, pos) ;
309301
fprintf(fff, "%c %c\n", c1, c2) ;
310-
311302
return ;
312-
313303
}
314304

315305
int checktriallelic(char *pc1, char *pc2, char x1, char x2)
@@ -391,7 +381,7 @@ int getdbname(char *dbase, char *name, char **pfqname)
391381
char ***names ;
392382
int n, k, t, i ;
393383

394-
fprintf(stderr, "call numlines in [getdbname]\n");
384+
// fprintf(stderr, "call numlines in [getdbname]\n");
395385
n = numlines(dbase) ;
396386

397387
ZALLOC(names, 3, char **) ;
@@ -435,15 +425,19 @@ int loadfa(char **poplist, int npops, FATYPE ***pfainfo, char *reg, int lopos, i
435425
region = (char*)malloc((23+strlen(reg))*sizeof(char));
436426
sprintf(region, "%s:%d-%d", reg, lo, hipos);
437427

438-
if (db == 0) refname = strcat(table_path, "Href.fa");
439-
else getdbname(iubfile, "Href", &refname);
428+
if (db == 0) {
429+
strcpy(refname, table_path);
430+
strcat(refname, "Href.fa");
431+
} else getdbname(iubfile, "Href", &refname);
440432

441433
if (ncall==1) {
442434
ZALLOC(falist, npops, char *) ;
443435
ZALLOC(famasklist, npops, char *) ;
444436
if (db == 0) {
445-
numfalist = setfalist(poplist, npops, ".fa", falist) ;
446-
t = setfalist(poplist, npops, ".filter.fa", famasklist) ;
437+
numfalist = setfalist(poplist, npops, ".ccomp.fa.rz", falist) ;
438+
t = setfalist(poplist, npops, ".ccompmask.fa.rz", famasklist) ;
439+
//numfalist = setfalist(poplist, npops, ".fa", falist) ;
440+
//t = setfalist(poplist, npops, ".filter.fa", famasklist) ;
447441
} else {
448442
numfalist = getfalist(poplist, npops, iubfile, falist) ; // set falist with the absolute path of hetfa files in .dblist file
449443
t = getfalist(poplist, npops, iubmaskfile, famasklist) ;
@@ -524,7 +518,6 @@ int loadfa(char **poplist, int npops, FATYPE ***pfainfo, char *reg, int lopos, i
524518
if (len_s==0) fatalx("bad fetch %s %s\n", fapt->faname, region); // fetch fai
525519

526520
len = len_r < len_s ? len_r : len_s;
527-
// len = len_s;
528521
if (rz == 1) // raziped
529522
for (i = 0; i < len; ++i)
530523
if (fapt->rstring[i] == 'Q') fapt->rstring[i] = ref[i];
@@ -620,7 +613,7 @@ int getiub(char *cc, char *ccmask, FATYPE **fainfo, char *reg, int pos)
620613

621614
if (newreg == YES) {
622615

623-
printf("zznewrrr %s :: %s %d %d %d\n",lastreg, regbuff, pos, lastlo, lasthi) ; fflush(stdout) ;
616+
// printf("zznewrrr %s :: %s %d %d %d\n",lastreg, regbuff, pos, lastlo, lasthi) ; fflush(stdout) ;
624617
fflush(stdout) ;
625618
freestring(&regname) ;
626619

@@ -632,18 +625,17 @@ int getiub(char *cc, char *ccmask, FATYPE **fainfo, char *reg, int pos)
632625
}
633626

634627
if (newpage == YES) {
635-
fprintf(stderr, "pos-getiub: %d\n", pos);
628+
//fprintf(stderr, "pos-getiub: %d\n", pos);
636629
fflush(stdout) ;
637630
lastlo = pos ;
638631
lasthi = pos + pagesize ;
639632
lastlo = MAX(lastlo, lopos) ;
640633
lasthi = MAX(lasthi, hipos) ;
641634
lasthi = MIN(lasthi, lastlo+pagesize) ;
642-
printf("calling loodfa %s %d %d \n", regname, lastlo, lasthi) ;
635+
// printf("calling loodfa %s %d %d \n", regname, lastlo, lasthi) ;
643636
fflush(stdout) ;
644-
//fprintf(stderr, "getiub: fapt->fai: %p\n", fainfo[0]->fai);
645637
loadfa(poplist, npops, &fainfo, regname, lastlo, lasthi) ;
646-
printf("newpage: %d %p %d %d\n", pos, topheap(), lastlo, lasthi) ;
638+
// printf("newpage: %d %p %d %d\n", pos, topheap(), lastlo, lasthi) ;
647639

648640
fflush(stdout) ;
649641
}
@@ -769,13 +761,14 @@ int setfalist(char **poplist, int npops, char *dbfile, char **iublist) {
769761
int t;
770762
for (t = 0; t < npops; ++t) {
771763
iublist[t] = strdup(table_path);
772-
iublist[t] = (char*) realloc(iublist[t], 64);
764+
iublist[t] = (char*) realloc(iublist[t], strlen(iublist[t]) + strlen(poplist[t]) + strlen(dbfile) + 1);
773765
iublist[t] = strcat(iublist[t], poplist[t]);
774-
if ((!strcmp (poplist[t], "Chimp") || !strcmp (poplist[t], "Href")) && strcmp (dbfile, ".fa")) {
775-
free (iublist[t]);
766+
767+
if ((!strcmp (poplist[t], "Chimp") || !strcmp (poplist[t], "Href")) && !strcmp (dbfile, ".ccomp.fa.rz"))
768+
iublist[t] = strcat(iublist[t], ".fa");
769+
else if ((!strcmp (poplist[t], "Chimp") || !strcmp (poplist[t], "Href")) && !strcmp (dbfile, ".ccompmask.fa.rz"))
776770
iublist[t] = "NULL";
777-
} else
778-
iublist[t] = strcat(iublist[t], dbfile);
771+
else iublist[t] = strcat(iublist[t], dbfile);
779772
}
780773
return npops;
781774
}
@@ -817,25 +810,8 @@ int getfalist(char **poplist, int npops, char *dbfile, char **iublist)
817810

818811
fclose(fff) ;
819812
return nx ;
820-
821813
}
822-
/*
823-
char *myfai_fetch(faidx_t *fai, char *reg, int *plen)
824-
{
825-
char *treg, *s ;
826-
treg = strdup(reg) ;
827814

828-
if (fai==NULL) fatalx("(my_fai_fetch): fai NULL\n") ;
829-
830-
s = fai_fetch(fai, treg, plen) ;
831-
if (*plen > 0) {
832-
free(treg) ;
833-
return s ;
834-
}
835-
free(treg) ;
836-
return NULL ;
837-
}
838-
*/
839815
void clearfainfo(FATYPE *fapt, int mode)
840816
{
841817

@@ -924,7 +900,6 @@ int abx(int a, int b)
924900
}
925901

926902
int abxok(int abx, int abxmode) {
927-
928903
int t ;
929904

930905
if (abxmode >= 10) {
@@ -933,8 +908,6 @@ int abxok(int abx, int abxmode) {
933908
return NO ;
934909
}
935910

936-
937-
938911
switch (abxmode) {
939912
case 0:
940913
return YES ;

0 commit comments

Comments
 (0)