Merge pull request #85 from pvanheus/rewrite_count_constant_sites

andrewjpage · web-flow · commit 657e74da6b5b · 2019-10-22T11:21:18.000+01:00
Rewrite how count_constant_sites is implemented, add test
diff --git a/VERSION b/VERSION
@@ -1 +1 @@
-2.5.0
+2.5.1
diff --git a/src/alignment-file.c b/src/alignment-file.c
@@ -105,15 +105,18 @@ void get_bases_for_each_snp(char filename[], char ** bases_for_snps)
   gzclose(fp);
 }
 
-void detect_snps(char filename[], int pure_mode, int output_monomorphic, int output_constant_site_counts)
+void detect_snps(char filename[], int pure_mode, int output_monomorphic) {
+    detect_snps_count_constant_sites(filename, pure_mode, output_monomorphic, NULL);
+}
+
+void detect_snps_count_constant_sites(char filename[], int pure_mode, int output_monomorphic, int* constant_site_counts)
 {
   int i;
   int l;
   number_of_snps = 0;
   number_of_samples = 0; 
   length_of_genome = 0;
   char * first_sequence;
-  int base_counts[] = {0, 0, 0, 0};
   /* array below allows quick mapping of A, C, T and G characters to indices in base_counts array */
   const int char_to_base_count_index[] = {-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 0, -1, 1, -1, -1, -1, 2, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 3};
 
@@ -195,15 +198,12 @@ void detect_snps(char filename[], int pure_mode, int output_monomorphic, int out
       {
           snp_locations[current_snp_index] = i;
           current_snp_index++;
-      } else if (is_pure(first_sequence[i])) {
-          base_counts[char_to_base_count_index[(int) toupper(first_sequence[i])]]++;
+      } else if (constant_site_counts != NULL && is_pure(first_sequence[i])) {
+          constant_site_counts[char_to_base_count_index[(int) toupper(first_sequence[i])]]++;
       }
 
   }
 
-  if (output_constant_site_counts)
-    printf("%d,%d,%d,%d\n", base_counts[0], base_counts[1], base_counts[2], base_counts[3]);
-
   free(first_sequence);
   kseq_destroy(seq);
   gzclose(fp);
diff --git a/src/alignment-file.h b/src/alignment-file.h
@@ -22,7 +22,8 @@
 
 #include "kseq.h"
 
-void detect_snps( char filename[], int pure_mode, int output_monomorphic, int output_constant_site_counts);
+void detect_snps( char filename[], int pure_mode, int output_monomorphic);
+void detect_snps_count_constant_sites(char filename[], int pure_mode, int output_monomorphic, int *constant_site_counts);
 void get_bases_for_each_snp(char filename[], char ** bases_for_snps);
 int is_unknown(char base);
 int get_length_of_genome();
diff --git a/src/main.c b/src/main.c
@@ -24,7 +24,6 @@
 #include <ctype.h>
 #include <unistd.h>
 #include <getopt.h>
-#include "alignment-file.h"
 #include "snp-sites.h"
 #include "config.h"
 
@@ -126,7 +125,7 @@ int main (int argc, char **argv) {
     strncpy(multi_fasta_filename, argv[optind], FILENAME_MAX); 
 
     if (output_constant_site_counts) {
-        detect_snps(multi_fasta_filename, pure_mode, output_monomorphic, output_constant_site_counts);
+        count_constant_sites(multi_fasta_filename, output_filename);
     } else if( pure_mode || output_monomorphic)
     {
       generate_snp_sites_with_ref_pure_mono(multi_fasta_filename,
diff --git a/src/snp-sites.c b/src/snp-sites.c
@@ -38,7 +38,7 @@ static int generate_snp_sites_generic(char filename[],
                                       int output_monomorphic)
 {
 	int i;
-	detect_snps(filename, pure_mode, output_monomorphic, 0);
+	detect_snps(filename, pure_mode, output_monomorphic);
 
   bases_for_snps =  calloc(get_number_of_snps()+1, sizeof(char*));
   
@@ -141,7 +141,32 @@ int generate_snp_sites_with_ref_pure_mono(char filename[],
                                     output_filename, output_reference, pure_mode, output_monomorphic);
 }
 
+void count_constant_sites(char multi_fasta_filename[], char output_filename[]) {
+    char cwd[100];
+    FILE *input_file;
+    FILE *output_file;
+    int *constant_site_counts = NULL;
 
+    output_file = (FILE *) fopen(output_filename, "w");
+    if (!output_file) {
+        fprintf(stderr, "ERROR: cannot open %s for writing: %s\n", output_filename, strerror(errno));
+        exit(EXIT_FAILURE);
+    }
+
+    constant_site_counts = (int *) calloc(4, sizeof(int));
+    if (constant_site_counts == NULL) {
+        fprintf(stderr, "ERROR: cannot allocated memory for constant_site_counts");
+        exit(EXIT_FAILURE);
+    }
+
+    detect_snps_count_constant_sites(multi_fasta_filename, 0, 0, constant_site_counts);
+
+
+    fprintf(output_file, "%d,%d,%d,%d\n", constant_site_counts[0], constant_site_counts[1],
+            constant_site_counts[2], constant_site_counts[3]);
+    fclose(output_file);
+    free(constant_site_counts);
+}
 // Inefficient
 void strip_directory_from_filename(char * input_filename, char * output_filename)
 {
diff --git a/src/snp-sites.h b/src/snp-sites.h
@@ -22,6 +22,9 @@
 #define _SNP_SITES_H_
 
 #include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <errno.h>
 
 int generate_snp_sites(char filename[],
                        int output_multi_fasta_file,
@@ -42,7 +45,9 @@ int generate_snp_sites_with_ref_pure_mono(char filename[],
                                 int output_reference,
                                 int pure_mode,
                                 int output_monomorphic);
-                                
+
+void count_constant_sites(char multi_fasta_filename[], char filename[]);
+
 void strip_directory_from_filename(char *input_filename,
                                    char *output_filename);
 
diff --git a/tests/check-snp-sites.c b/tests/check-snp-sites.c
@@ -183,64 +183,64 @@ END_TEST
 
 START_TEST (valid_genome_length)
 {
-  detect_snps("../tests/data/alignment_file_one_line_per_sequence.aln",0,0,0);
+  detect_snps("../tests/data/alignment_file_one_line_per_sequence.aln",0,0);
   fail_unless( get_length_of_genome() == 2000 );
 }
 END_TEST
 
 START_TEST (valid_genome_length_with_multiple_lines_per_sequence)
 {
-  detect_snps("../tests/data/alignment_file_multiple_lines_per_sequence.aln",0,0,0);
+  detect_snps("../tests/data/alignment_file_multiple_lines_per_sequence.aln",0,0);
   fail_unless( get_length_of_genome() == 2000 );
 }
 END_TEST
 
 START_TEST (valid_number_of_sequences_in_file)
 {
-  detect_snps("../tests/data/alignment_file_one_line_per_sequence.aln",0,0,0);
+  detect_snps("../tests/data/alignment_file_one_line_per_sequence.aln",0,0);
   fail_unless( get_number_of_samples() == 109 );
 }
 END_TEST
 
 START_TEST (valid_number_of_sequences_in_file_with_multiple_lines_per_sequence)
 {
-  detect_snps("../tests/data/alignment_file_multiple_lines_per_sequence.aln",0,0,0);
+  detect_snps("../tests/data/alignment_file_multiple_lines_per_sequence.aln",0,0);
   fail_unless( get_number_of_samples() == 109 );
 }
 END_TEST 
 
 START_TEST (number_of_snps_detected)
 {
-  detect_snps("../tests/data/alignment_file_multiple_lines_per_sequence.aln",0,0,0);
+  detect_snps("../tests/data/alignment_file_multiple_lines_per_sequence.aln",0,0);
   fail_unless( get_number_of_snps()  == 5);
 }
 END_TEST
 
 START_TEST (number_of_snps_detected_small)
 {
-  detect_snps("../tests/data/small_alignment.aln",0,0,0);
+  detect_snps("../tests/data/small_alignment.aln",0,0);
   fail_unless(  get_number_of_snps()  == 1);
 }
 END_TEST
   
 START_TEST (detect_snps_pure_mode)
 {
-  detect_snps("../tests/data/pure_mode_alignment.aln",1,0,0);
+  detect_snps("../tests/data/pure_mode_alignment.aln",1,0);
   fail_unless(  get_number_of_snps()  == 2);
 }
 END_TEST
   
   
 START_TEST (detect_snps_pure_mode_monomorphic)
 {
-  detect_snps("../tests/data/pure_mode_monomorphic_alignment.aln",1,1,0);
+  detect_snps("../tests/data/pure_mode_monomorphic_alignment.aln",1,1);
   fail_unless(  get_number_of_snps()  == 3);
 }
 END_TEST
 
 START_TEST (sample_names_from_alignment_file)
 {
-  detect_snps("../tests/data/small_alignment.aln",0,0,0);
+  detect_snps("../tests/data/small_alignment.aln",0,0);
   char ** current_sequence_names = get_sequence_names();
 
   fail_unless(strcmp(current_sequence_names[0],"reference_sequence") == 0);
@@ -269,6 +269,14 @@ START_TEST (check_strip_directory_from_filename_with_directory)
 }
 END_TEST
 
+START_TEST (check_count_constant_sites)
+{
+    count_constant_sites("../tests/data/small_alignment.aln", "small_alignment.constant_site_counts.txt");
+    fail_unless(compare_files("../tests/data/small_alignment.constant_site_counts.txt", "small_alignment.constant_site_counts.txt"));
+    remove("small_alignment.constant_site_counts.txt");
+}
+END_TEST
+
 Suite * snp_sites_suite (void)
 {
   Suite *s = suite_create ("Creating_SNP_Sites");
@@ -302,7 +310,8 @@ Suite * snp_sites_suite (void)
   tcase_add_test (tc_snp_sites, valid_phylip_plus_reference);
   tcase_add_test (tc_snp_sites, valid_alignment_with_pure_mode);
   tcase_add_test (tc_snp_sites, valid_alignment_with_monomorphic_sites);
-  
+  tcase_add_test (tc_snp_sites, check_count_constant_sites);
+
   tcase_add_exit_test(tc_snp_sites, invalid_with_uneven_file_lengths,EXIT_FAILURE);
   remove("uneven_alignment.aln.snp_sites.aln");
   suite_add_tcase (s, tc_snp_sites);
diff --git a/tests/data/small_alignment.constant_site_counts.txt b/tests/data/small_alignment.constant_site_counts.txt
@@ -0,0 +1 @@
+2,2,1,2

Original file line number	Diff line number	Diff line change
`@@ -183,64 +183,64 @@ END_TEST`
`183`	`183`
`184`	`184`	`START_TEST (valid_genome_length)`
`185`	`185`	`{`
`186`		`- detect_snps("../tests/data/alignment_file_one_line_per_sequence.aln",0,0,0);`
	`186`	`+ detect_snps("../tests/data/alignment_file_one_line_per_sequence.aln",0,0);`
`187`	`187`	`fail_unless( get_length_of_genome() == 2000 );`
`188`	`188`	`}`
`189`	`189`	`END_TEST`
`190`	`190`
`191`	`191`	`START_TEST (valid_genome_length_with_multiple_lines_per_sequence)`
`192`	`192`	`{`
`193`		`- detect_snps("../tests/data/alignment_file_multiple_lines_per_sequence.aln",0,0,0);`
	`193`	`+ detect_snps("../tests/data/alignment_file_multiple_lines_per_sequence.aln",0,0);`
`194`	`194`	`fail_unless( get_length_of_genome() == 2000 );`
`195`	`195`	`}`
`196`	`196`	`END_TEST`
`197`	`197`
`198`	`198`	`START_TEST (valid_number_of_sequences_in_file)`
`199`	`199`	`{`
`200`		`- detect_snps("../tests/data/alignment_file_one_line_per_sequence.aln",0,0,0);`
	`200`	`+ detect_snps("../tests/data/alignment_file_one_line_per_sequence.aln",0,0);`
`201`	`201`	`fail_unless( get_number_of_samples() == 109 );`
`202`	`202`	`}`
`203`	`203`	`END_TEST`
`204`	`204`
`205`	`205`	`START_TEST (valid_number_of_sequences_in_file_with_multiple_lines_per_sequence)`
`206`	`206`	`{`
`207`		`- detect_snps("../tests/data/alignment_file_multiple_lines_per_sequence.aln",0,0,0);`
	`207`	`+ detect_snps("../tests/data/alignment_file_multiple_lines_per_sequence.aln",0,0);`
`208`	`208`	`fail_unless( get_number_of_samples() == 109 );`
`209`	`209`	`}`
`210`	`210`	`END_TEST`
`211`	`211`
`212`	`212`	`START_TEST (number_of_snps_detected)`
`213`	`213`	`{`
`214`		`- detect_snps("../tests/data/alignment_file_multiple_lines_per_sequence.aln",0,0,0);`
	`214`	`+ detect_snps("../tests/data/alignment_file_multiple_lines_per_sequence.aln",0,0);`
`215`	`215`	`fail_unless( get_number_of_snps() == 5);`
`216`	`216`	`}`
`217`	`217`	`END_TEST`
`218`	`218`
`219`	`219`	`START_TEST (number_of_snps_detected_small)`
`220`	`220`	`{`
`221`		`- detect_snps("../tests/data/small_alignment.aln",0,0,0);`
	`221`	`+ detect_snps("../tests/data/small_alignment.aln",0,0);`
`222`	`222`	`fail_unless( get_number_of_snps() == 1);`
`223`	`223`	`}`
`224`	`224`	`END_TEST`
`225`	`225`
`226`	`226`	`START_TEST (detect_snps_pure_mode)`
`227`	`227`	`{`
`228`		`- detect_snps("../tests/data/pure_mode_alignment.aln",1,0,0);`
	`228`	`+ detect_snps("../tests/data/pure_mode_alignment.aln",1,0);`
`229`	`229`	`fail_unless( get_number_of_snps() == 2);`
`230`	`230`	`}`
`231`	`231`	`END_TEST`
`232`	`232`
`233`	`233`
`234`	`234`	`START_TEST (detect_snps_pure_mode_monomorphic)`
`235`	`235`	`{`
`236`		`- detect_snps("../tests/data/pure_mode_monomorphic_alignment.aln",1,1,0);`
	`236`	`+ detect_snps("../tests/data/pure_mode_monomorphic_alignment.aln",1,1);`
`237`	`237`	`fail_unless( get_number_of_snps() == 3);`
`238`	`238`	`}`
`239`	`239`	`END_TEST`
`240`	`240`
`241`	`241`	`START_TEST (sample_names_from_alignment_file)`
`242`	`242`	`{`
`243`		`- detect_snps("../tests/data/small_alignment.aln",0,0,0);`
	`243`	`+ detect_snps("../tests/data/small_alignment.aln",0,0);`
`244`	`244`	`char ** current_sequence_names = get_sequence_names();`
`245`	`245`
`246`	`246`	`fail_unless(strcmp(current_sequence_names[0],"reference_sequence") == 0);`
`@@ -269,6 +269,14 @@ START_TEST (check_strip_directory_from_filename_with_directory)`
`269`	`269`	`}`
`270`	`270`	`END_TEST`
`271`	`271`
	`272`	`+START_TEST (check_count_constant_sites)`
	`273`	`+{`
	`274`	`+ count_constant_sites("../tests/data/small_alignment.aln", "small_alignment.constant_site_counts.txt");`
	`275`	`+ fail_unless(compare_files("../tests/data/small_alignment.constant_site_counts.txt", "small_alignment.constant_site_counts.txt"));`
	`276`	`+ remove("small_alignment.constant_site_counts.txt");`
	`277`	`+}`
	`278`	`+END_TEST`
	`279`	`+`
`272`	`280`	`Suite * snp_sites_suite (void)`
`273`	`281`	`{`
`274`	`282`	`Suite *s = suite_create ("Creating_SNP_Sites");`
`@@ -302,7 +310,8 @@ Suite * snp_sites_suite (void)`
`302`	`310`	`tcase_add_test (tc_snp_sites, valid_phylip_plus_reference);`
`303`	`311`	`tcase_add_test (tc_snp_sites, valid_alignment_with_pure_mode);`
`304`	`312`	`tcase_add_test (tc_snp_sites, valid_alignment_with_monomorphic_sites);`
`305`		`-`
	`313`	`+ tcase_add_test (tc_snp_sites, check_count_constant_sites);`
	`314`	`+`
`306`	`315`	`tcase_add_exit_test(tc_snp_sites, invalid_with_uneven_file_lengths,EXIT_FAILURE);`
`307`	`316`	`remove("uneven_alignment.aln.snp_sites.aln");`
`308`	`317`	`suite_add_tcase (s, tc_snp_sites);`