Merge pull request #50 from andrewjpage/only_acgt

allow pure bases and monomorphic, making it work with BEAST

Author: andrewjpage

VERSION

@@ -1 +1 @@
-2.2.2
+2.3.0

snp-sites.txt

@@ -5,7 +5,7 @@ SNP-SITES(1)
 
 NAME
 ----
-snp-sites - finds snp sites from a multi fasta alignment file
+snp-sites - finds SNP sites from a multi FASTA alignment file
 
 SYNOPSIS
 --------
@@ -33,6 +33,12 @@ OPTIONS
 *-p*::
 	Output a phylip file
 
+*-c*::
+	Only output columns containing exclusively ACGT
+
+*-b*::
+	Output monomorphic sites, useful for BEAST
+
 *-o*::
 	Specify an output filename
 
@@ -49,6 +55,8 @@ EXAMPLES
  snp-sites my-alignment.aln
 
  snp-sites my-gzipped-alignment.aln.gz
+ 
+ snp-sites -cb -o output_beast_format.aln inputfile.aln
 
 
 FORMAT OF THE INPUT FILE

src/alignment-file.c

@@ -105,7 +105,7 @@ void get_bases_for_each_snp(char filename[], char ** bases_for_snps)
   gzclose(fp);
 }
 
-void detect_snps(char filename[])
+void detect_snps(char filename[], int pure_mode, int output_monomorphic)
 {
   int i;
   int l;
@@ -134,16 +134,27 @@ void detect_snps(char filename[])
     }
     for(i = 0; i < length_of_genome; i++)
     {
+      if(first_sequence[i] == '#')
+      {
+        continue;
+      }
+        
       if(first_sequence[i] == 'N' && !is_unknown(seq->seq.s[i]))
       {
 	        first_sequence[i] = toupper(seq->seq.s[i]);
           pseudo_reference_sequence[i] = toupper(seq->seq.s[i]);
       }
-	  
-	  if(first_sequence[i] != '>' && !is_unknown(seq->seq.s[i]) && first_sequence[i] != 'N' && first_sequence[i] != toupper(seq->seq.s[i]))
+     
+      // in pure mode we only want /ACGT/i, if any other base is found the whole column is excluded
+      if(pure_mode && !is_pure(seq->seq.s[i]))
+      {
+        first_sequence[i] = '#';
+        continue;
+      }
+      
+	 	  if(first_sequence[i] != '>' && !is_unknown(seq->seq.s[i]) && first_sequence[i] != 'N' && first_sequence[i] != toupper(seq->seq.s[i]))
 	  {
 	      first_sequence[i] = '>';
-	      number_of_snps++;
 	  }
    }
 
@@ -157,18 +168,26 @@ void detect_snps(char filename[])
    number_of_samples++;
   }
 
+  for(i = 0; i < length_of_genome; i++)
+  {
+    if(first_sequence[i] == '>' || (output_monomorphic && first_sequence[i] != '#'))
+    {
+      number_of_snps++;
+    }
+  }
+  
   if(number_of_snps == 0)
   {
     fprintf(stderr, "Warning: No SNPs were detected so there is nothing to output.\n");
     fflush(stderr);
     exit(EXIT_FAILURE);
   }
-  
+
   int current_snp_index = 0;
   snp_locations = calloc(number_of_snps, sizeof(int));
   for(i = 0; i < length_of_genome; i++)
   {
-      if(first_sequence[i] == '>')
+      if(first_sequence[i] == '>' || (output_monomorphic && first_sequence[i] != '#'))
       {
           snp_locations[current_snp_index] = i;
           current_snp_index++;
@@ -192,3 +211,20 @@ int is_unknown(char base)
       return 0;
   }
 }
+
+int is_pure(char base)
+{
+  switch (base) {
+    case 'A':
+    case 'C':
+    case 'G':
+    case 'T':
+    case 'a':
+    case 'c':
+    case 'g':
+    case 't':
+      return 1;
+    default:
+      return 0;
+  }
+}

src/alignment-file.h

@@ -22,7 +22,7 @@
 
 #include "kseq.h"
 
-void detect_snps( char filename[]);
+void detect_snps( char filename[], int pure_mode, int output_monomorphic);
 void get_bases_for_each_snp(char filename[], char ** bases_for_snps);
 int is_unknown(char base);
 int get_length_of_genome();
@@ -31,6 +31,7 @@ int get_number_of_snps();
 char ** get_sequence_names();
 int * get_snp_locations();
 char * get_pseudo_reference_sequence();
+int is_pure(char base);
 
 #define MAX_SAMPLE_NAME_SIZE 2048
 #define DEFAULT_NUM_SAMPLES 65536

src/main.c

@@ -32,16 +32,21 @@
 
 static void print_usage()
 {
-	printf("Usage: snp-sites [-mvph] [-o output_filename] <file>\n");
-	printf("This program finds snp sites from a multi FASTA alignment file.\n");
-	printf(" -r		output internal pseudo reference sequence\n");
-  printf(" -m		output a multi fasta alignment file (default)\n");
-	printf(" -v		output a VCF file\n");
-	printf(" -p		output a phylip file\n");
-	printf(" -o		specify an output filename\n");
-	printf(" -h		this help message\n");
-	printf(" -V		print version and exit\n");
-	printf(" <file>		input alignment file which can optionally be gzipped\n\n");
+  printf("Usage: snp-sites [-rmvpcbhV] [-o output_filename] <file>\n");
+  printf("This program finds snp sites from a multi FASTA alignment file.\n");
+  printf(" -r     output internal pseudo reference sequence\n");
+  printf(" -m     output a multi fasta alignment file (default)\n");
+  printf(" -v     output a VCF file\n");
+  printf(" -p     output a phylip file\n");
+  printf(" -o STR specify an output filename\n");
+  printf(" -c     only output columns containing exclusively ACGT\n");
+  printf(" -b     output monomorphic sites, used for BEAST\n");
+  printf(" -h     this help message\n");
+  printf(" -V     print version and exit\n");
+  printf(" <file> input alignment file which can optionally be gzipped\n\n");
+  
+  printf("Example: creating files for BEAST\n");
+  printf(" snp-sites -cb -o outputfile.aln inputfile.aln\n\n");
 
   printf("If you use this program, please cite:\n");
   printf("\"SNP-sites: rapid efficient extraction of SNPs from multi-FASTA alignments\",\n");
@@ -65,8 +70,10 @@ int main (int argc, char **argv) {
   int output_vcf_file = 0;
   int output_phylip_file = 0;
   int output_reference = 0;
-	
-	 while ((c = getopt (argc, argv, "mvrpo:V")) != -1)
+	int pure_mode = 0;
+  int output_monomorphic =0;
+  
+	 while ((c = getopt (argc, argv, "mvrbpco:V")) != -1)
       switch (c)
         {
         case 'm':
@@ -83,6 +90,12 @@ int main (int argc, char **argv) {
           break;
         case 'r':
           output_reference = 1;
+          break;
+        case 'c':
+          pure_mode = 1;
+          break;
+        case 'b':
+          output_monomorphic = 1;
           break;
 	      case 'o':
           strncpy(output_filename, optarg, FILENAME_MAX);
@@ -105,7 +118,15 @@ int main (int argc, char **argv) {
     }
 
     strncpy(multi_fasta_filename, argv[optind], FILENAME_MAX); 
-    if (output_reference) {
+    
+    if( pure_mode || output_monomorphic)
+    {
+      generate_snp_sites_with_ref_pure_mono(multi_fasta_filename,
+                                  output_multi_fasta_file,
+                                  output_vcf_file, output_phylip_file,
+                                  output_filename,output_reference,pure_mode,output_monomorphic);
+    }
+    else if(output_reference ) {
       generate_snp_sites_with_ref(multi_fasta_filename,
                                   output_multi_fasta_file,
                                   output_vcf_file, output_phylip_file,

src/snp-sites.c

@@ -32,10 +32,10 @@ static int generate_snp_sites_generic(char filename[],
                                       int output_vcf_file,
                                       int output_phylip_file,
                                       char output_filename[],
-                                      int output_reference)
+                                      int output_reference, int pure_mode, int output_monomorphic)
 {
 	int i;
-	detect_snps(filename);
+	detect_snps(filename, pure_mode, output_monomorphic);
 
 	char* bases_for_snps[get_number_of_snps()];
 
@@ -112,7 +112,7 @@ int generate_snp_sites(char filename[], int output_multi_fasta_file,
 {
   return generate_snp_sites_generic(filename, output_multi_fasta_file,
                                     output_vcf_file, output_phylip_file,
-                                    output_filename, 0);
+                                    output_filename, 0,0,0);
 }
 
 int generate_snp_sites_with_ref(char filename[], int output_multi_fasta_file,
@@ -121,9 +121,24 @@ int generate_snp_sites_with_ref(char filename[], int output_multi_fasta_file,
 {
     return generate_snp_sites_generic(filename, output_multi_fasta_file,
                                       output_vcf_file, output_phylip_file,
-                                      output_filename, 1);
+                                      output_filename, 1,0,0);
 }
 
+int generate_snp_sites_with_ref_pure_mono(char filename[],
+                                int output_multi_fasta_file,
+                                int output_vcf_file,
+                                int output_phylip_file,
+                                char output_filename[],
+                                int output_reference,
+                                int pure_mode,
+                                int output_monomorphic)
+{
+  return generate_snp_sites_generic(filename, output_multi_fasta_file,
+                                    output_vcf_file, output_phylip_file,
+                                    output_filename, output_reference, pure_mode, output_monomorphic);
+}
+
+
 // Inefficient
 void strip_directory_from_filename(char * input_filename, char * output_filename)
 {

src/snp-sites.h

@@ -33,6 +33,16 @@ int generate_snp_sites_with_ref(char filename[],
                                 int output_vcf_file,
                                 int output_phylip_file,
                                 char output_filename[]);
+                                
+int generate_snp_sites_with_ref_pure_mono(char filename[],
+                                int output_multi_fasta_file,
+                                int output_vcf_file,
+                                int output_phylip_file,
+                                char output_filename[],
+                                int output_reference,
+                                int pure_mode,
+                                int output_monomorphic);
+                                
 void strip_directory_from_filename(char *input_filename,
                                    char *output_filename);
 

src/vcf.c

@@ -113,6 +113,7 @@ void output_vcf_row(FILE * vcf_file_pointer, char * bases_for_snp, int snp_locat
 	free(alt_bases);
 
 	fprintf( vcf_file_pointer, "\n");	
+  
 }
 
 
@@ -145,14 +146,15 @@ char * alternative_bases(char reference_base, char * bases_for_snp, int number_o
 			}
 		}
 	}
+  alt_bases[num_alt_bases] = '\0';
 	return alt_bases;
 }
 
 char * format_allele_index(char base, char reference_base, char * alt_bases)
 {
 	int length_of_alt_bases = strlen(alt_bases);
 	assert(length_of_alt_bases < 100);
-	char * result = calloc(3, sizeof(char));
+	char * result = calloc(5, sizeof(char));
 	int index;
 
   if(is_unknown(base))
@@ -183,6 +185,14 @@ char * format_alternative_bases(char * alt_bases)
 {
 	int number_of_alt_bases = strlen(alt_bases);
 	assert( number_of_alt_bases < MAXIMUM_NUMBER_OF_ALT_BASES );
+  
+  if(number_of_alt_bases ==0)
+  {
+    char * formatted_alt_bases = calloc(2 + 1, sizeof(char));
+    formatted_alt_bases[0] = '.';
+    return formatted_alt_bases;
+  }
+
 	char * formatted_alt_bases = calloc(number_of_alt_bases*2 + 1, sizeof(char));
 	int i;
 	formatted_alt_bases[0] = alt_bases[0];