merge snp2h5 changes from development branch

Author: gmcvicker

snp2h5/Makefile

@@ -1,5 +1,5 @@
 # Edit the following variables as needed
-HDF_INSTALL=$(HOME)/anaconda
+HDF_INSTALL=$(HOME)/anaconda2
 EXTLIB=-L$(HDF_INSTALL)/lib
 CC=gcc
 LIB=-lz -lm -lhdf5 -lhdf5_hl

snp2h5/err.c

@@ -20,7 +20,7 @@ void my_err(const char *format, ...) {
   }
 
   va_end(args);
-  abort();
+  exit(2);
 }
 
 

snp2h5/snp.h

@@ -2,14 +2,18 @@
 #define __SNP_H__
 
 
-#define SNP_MAX_ALLELE 32
-#define SNP_MAX_CHROM 32
+/* maximum sequence length of alleles (longer are truncated in snp_tab) */
+#define SNP_MAX_ALLELE 100
+
+/* maximum length of chromosome name */
+/* #define SNP_MAX_CHROM 32 */
+/* maximum length of SNP identifier */
 #define SNP_MAX_NAME 16
 
 
 typedef struct {
   char name[SNP_MAX_NAME];
-  char chrom[SNP_MAX_CHROM];
+  /*   char chrom[SNP_MAX_CHROM]; */
   long pos;
   char allele1[SNP_MAX_ALLELE];
   char allele2[SNP_MAX_ALLELE];

snp2h5/snp2h5.c

@@ -658,7 +658,8 @@ void parse_impute(Arguments *args, Chromosome *all_chroms, int n_chrom,
     /* open file, count number of lines */
     fprintf(stderr, "reading from file %s\n", imp_files[i]);
     gzf = util_must_gzopen(args->input_files[i], "rb");
-    set_file_info(gzf, args->input_files[i], args, &file_info, NULL, impute_info);
+    set_file_info(gzf, args->input_files[i], args, &file_info, NULL,
+		  impute_info);
 
     /* initialize output files and memory to hold genotypes,
      * haplotypes, etc
@@ -813,7 +814,7 @@ void parse_impute(Arguments *args, Chromosome *all_chroms, int n_chrom,
 
 	  line_num += 1;
 
-	  if((line_num % 1000) == 0) {
+	  if((line_num % 10000) == 0) {
 	    fprintf(stderr, ".");
 	  }
 	}
@@ -946,10 +947,9 @@ void parse_vcf(Arguments *args, Chromosome *all_chroms, int n_chrom,
 
     fprintf(stderr, "parsing file and writing to HDF5 files\n");
 
-    
     while(vcf_read_line(gzf, vcf, &snp,
 			geno_probs, haplotypes) != -1) {
-
+      
       if(geno_probs) {
 	write_h5matrix_row(gprob_info, row, geno_probs);
       }
@@ -975,7 +975,7 @@ void parse_vcf(Arguments *args, Chromosome *all_chroms, int n_chrom,
       if(snp_tab) {
 	snp_tab_append_row(snp_tab, &snp);
       }
-      
+
       row++;
       if((row % 1000) == 0) {
 	fprintf(stderr, ".");
@@ -1030,6 +1030,8 @@ int main(int argc, char **argv) {
 
   /* read chromosomes */
   all_chroms = chrom_read_file(args.chrom_file, &n_chrom);
+
+  fprintf(stderr, "long alleles will be truncated to %dbp\n", SNP_MAX_ALLELE);
 
   /* create new HDF5 file(s) */
   if(args.geno_prob_file) {

snp2h5/snptab.c

@@ -35,22 +35,30 @@ SNPTab *snp_tab_new(hid_t h5file, const char *chrom_name,
   /* set datatypes for each field */
   tab->name_type = H5Tcopy(H5T_C_S1);
   H5Tset_size(tab->name_type, SNP_MAX_NAME);
+  /* no longer store chromosome as each chromosome
+   * gets its own table 
+  /* tab->chrom_type = H5Tcopy(H5T_C_S1);
+   * H5Tset_size(tab->chrom_type, SNP_MAX_CHROM);
+   */
   tab->allele_type = H5Tcopy(H5T_C_S1);
-  H5Tset_size(tab->allele_type, SNP_MAX_ALLELE);    
+  H5Tset_size(tab->allele_type, SNP_MAX_ALLELE);
   tab->field_type[0] = tab->name_type; /* name */
+  /* tab->field_type[1] = tab->chrom_type; */ /* chromosome */
   tab->field_type[1] = H5T_NATIVE_LONG; /* pos */
   tab->field_type[2] = tab->allele_type; /* allele1 */
   tab->field_type[3] = tab->allele_type; /* allele2 */
 
   /* sizes of record and each field */
-  tab->record_size = sizeof(SNPTab);
+  tab->record_size = sizeof(SNP);
   tab->field_size[0] = sizeof(snp_desc.name);
+  /* tab->field_size[1] = sizeof(snp_desc.chrom); */
   tab->field_size[1] = sizeof(snp_desc.pos);
   tab->field_size[2] = sizeof(snp_desc.allele1);
   tab->field_size[3] = sizeof(snp_desc.allele2);
-
+    
   /* offsets of each field */
   tab->field_offset[0] = HOFFSET(SNP, name);
+  /* tab->field_offset[1] = HOFFSET(SNP, chrom); */
   tab->field_offset[1] = HOFFSET(SNP, pos);
   tab->field_offset[2] = HOFFSET(SNP, allele1);
   tab->field_offset[3] = HOFFSET(SNP, allele2);

snp2h5/snptab.h

@@ -6,16 +6,17 @@
 #include "snp.h"
 
 
-#define SNPTAB_NFIELDS 5
 #define SNPTAB_N_FIELDS 4
 
 /* chunk size affects performance a lot
  * small chunks = much faster writing of tables, but
  * worse compression
  */
+
 #define SNPTAB_CHUNK_SIZE 12
 
 
+
 /* SNPTab holds information about SNP table 
  * and datatypes of each field.
  */
@@ -30,6 +31,7 @@ typedef struct {
   size_t field_offset[SNPTAB_N_FIELDS];
 
   hid_t name_type;
+  hid_t chrom_type;
   hid_t allele_type;
 
   int compress;

snp2h5/util.c

@@ -71,18 +71,16 @@ char *util_read_entire_file(char *filename) {
  * before and after the count of newlines.
  */
 long util_fcount_lines(FILE *fh) {
-  char buf[UTIL_FGETS_BUF_SZ];
   long line_count;
-  int len;
+  char c;
 
   if(fseek(fh, 0L, SEEK_SET) != 0) {
     my_err("%s:%d: could not rewind filehandle", __FILE__, __LINE__);
   }
 
   line_count = 0;
-  while(fgets(buf, UTIL_FGETS_BUF_SZ, fh) != NULL) {
-    len = strlen(buf);
-    if(buf[len-1] == '\n') {
+  while((c = getc(fh)) != EOF) {
+    if(c == '\n') {
       line_count++;
     }
   }
@@ -95,6 +93,9 @@ long util_fcount_lines(FILE *fh) {
 }
 
 
+
+
+
 /**
  * write a subset of lines which are flagged as TRUE in the provided
  * array (that should have a length corresponding to the number lines
@@ -151,6 +152,7 @@ long util_count_lines(const char *filename) {
   return line_count;
 }
 
+
 /**
  * Counts the number of '\n' characters in the provided gzipped file.
  * The gzFile is rewound to the beginning of the file before and after the count of
@@ -572,40 +574,23 @@ char *util_str_ndup(const char *str, size_t n) {
  * This function still gives no indication that src string is
  * truncated, however.
  */
-size_t util_strncpy(char *dest, const char *src, size_t n) {
-  size_t i, len;
-  int src_end;
 
-  len = 0;
-  src_end = FALSE;
-  
-  for(i = 0; i < n-1; i++) {
-    if(src_end) {
-      /* Continue NULL padding to end of dest buf with '\0'.
-       * This is to have consistent behaviour with strncpy and also
-       * HDF5 API seems to expect this. WHen HDF5 string written
-       * with non-null chars past terminating '\0', entire length of
-       * string buffer is printed by pytables (including unprintable 
-       * binary chars) when string retrieved.
-       */
-      dest[i] = '\0';
-    } else {		 
+size_t util_strncpy(char *dest, const char *src, size_t n) {
+  size_t i = 0;
+  if(n > 0) {
+    while((i < n-1) && (src[i] != '\0')) {
       dest[i] = src[i];
-      if(src[i] == '\0') {
-	/* hit end of src string */
-	src_end = TRUE;
-      } else {
-	len += 1;
-      }
+      i++;
     }
+    /* padd to end with \0 */
+    memset(&dest[i], '\0', n-i);
+    /* dest[i] = '\0'; */
   }
-
-  dest[n-1] = '\0';
-  return len;
+	 
+  return i;
 }
 
 
-
 /**
  * Does an in-place uppercase of all characters in a string
  */

snp2h5/vcf.c

@@ -350,7 +350,11 @@ int vcf_read_line(gzFile vcf_fh, VCFInfo *vcf_info, SNP *snp,
   if(token == NULL) {
     my_err("expected at least %d tokens per line\n", n_fix_header);
   }
-  util_strncpy(snp->chrom, token, sizeof(snp->chrom));
+
+  /* we don't bother to store chromosome since we store 
+   * SNPs from each chromosome in their own table
+   */
+  /* util_strncpy(snp->chrom, token, sizeof(snp->chrom)); */
 
 
   /* pos */
@@ -372,26 +376,31 @@ int vcf_read_line(gzFile vcf_fh, VCFInfo *vcf_info, SNP *snp,
   if(token == NULL) {
     my_err("expected at least %d tokens per line\n", n_fix_header);
   }
-  vcf_info->ref_len = strlen(token);
   ref_len = util_strncpy(snp->allele1, token, sizeof(snp->allele1));
 
-  if(ref_len != vcf_info->ref_len) {
-    my_warn("truncating long allele (%ld bp) to %ld bp\n",
-	    vcf_info->ref_len, ref_len);
-  }
+  /* used to warn about truncations, but makes program too
+   * chatty if there are a lot of them
+   */
+  vcf_info->ref_len = 0;
+  /* vcf_info->ref_len = strlen(token); */
+  /* if(ref_len != vcf_info->ref_len) { */
+  /*   my_warn("truncating long allele (%ld bp) to %ld bp\n", */
+  /* 	    vcf_info->ref_len, ref_len); */
+  /* } */
 
   /* alt */
   token = strsep(&cur, delim);
   if(token == NULL) {
     my_err("expected at least %d tokens per line\n", n_fix_header);
   }
-  vcf_info->alt_len = strlen(token);
   alt_len = util_strncpy(snp->allele2, token, sizeof(snp->allele2));
-
-  if(alt_len != vcf_info->alt_len) {
-    my_warn("truncating long allele (%ld bp) to %ld bp\n",
-	    vcf_info->alt_len, alt_len);
-  }
+  
+  vcf_info->alt_len = 0;
+  /* vcf_info->alt_len = strlen(token); */
+  /* if(alt_len != vcf_info->alt_len) { */
+  /*   my_warn("truncating long allele (%ld bp) to %ld bp\n", */
+  /* 	    vcf_info->alt_len, alt_len); */
+  /* } */
 
   /* qual */
   token = strsep(&cur, delim);