main.c
/*
* Wellcome Trust Sanger Institute
* Copyright (C) 2011 Wellcome Trust Sanger Institute
*
* This program is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License
* as published by the Free Software Foundation; either version 2
* of the License, or (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <getopt.h>
#include <unistd.h>
#include "snp_sites.h"
#include "gubbins.h"
#include "../config.h"
#include "string_cat.h"
#define MAX_FILENAME_SIZE 1024
const char* program_name;
// Assumptions:
// The sequences in the multi fasta alignment file are the same length
// Your only interested in SNPs, INDELS are ignored
// The first sequence is chosen as the reference sequence
// If there is an indel in the reference sequence, the first normal base found in another strain is used.
void print_usage(FILE* stream, int exit_code)
{
fprintf (stream, "This program is not supposed to be directly run. Use run_gubbins.py instead\n");
fprintf (stream, "Usage: %s [options] alignment_file\n", program_name);
fprintf (stream, "Version: %s\n", PACKAGE_VERSION);
fprintf (stream,
" -r detect recombinations mode\n"
" -t Newick tree file\n"
" -v VCF file\n"
" -f Original Multifasta file\n"
" -m Min SNPs for identifying a recombination block\n"
" -a Min window size\n"
" -b Max window size\n"
" -p p value for detecting recombinations\n"
" -i p value ratio for trimming recombinations\n"
" -h Display this usage information.\n\n"
);
exit (exit_code);
}
int check_file_exists_or_exit(char * filename)
{
if (access( filename, F_OK ) != -1 )
{
return 1;
}
else
{
printf("Error: File '%s' does not exist\n",filename);
print_usage(stderr, EXIT_FAILURE);
return 0;
}
}
int main (argc, argv) int argc; char **argv;
{
int c;
char multi_fasta_filename[MAX_FILENAME_SIZE] = {""};
char vcf_filename[MAX_FILENAME_SIZE] = {""};
char tree_filename[MAX_FILENAME_SIZE] = {""};
char original_multi_fasta_filename[MAX_FILENAME_SIZE] = {""};
int recombination_flag = 0 ;
int min_snps = 3;
int window_min = 100;
int window_max = 10000;
float uncorrected_p_value = 0.05;
float trimming_ratio = 1.0;
int extensive_search_flag = 0;
program_name = argv[0];
while (1)
{
static struct option long_options[] =
{
{"help", no_argument, 0, 'h'},
{"recombination", no_argument, 0, 'r'},
{"vcf", required_argument, 0, 'v'},
{"tree", required_argument, 0, 't'},
{"original_multifasta", required_argument, 0, 'f'},
{"min_snps", required_argument, 0, 'm'},
{"window_min", required_argument, 0, 'a'},
{"window_max", required_argument, 0, 'b'},
{"p_value", required_argument, 0, 'p'},
{"trimming_ratio", required_argument, 0, 'i'},
{"extended_search", required_argument, 0, 'x'},
{0, 0, 0, 0}
};
/* getopt_long stores the option index here. */
int option_index = 0;
c = getopt_long (argc, argv, "hrxv:f:t:m:a:b:p:i:",
long_options, &option_index);
/* Detect the end of the options. */
if (c == -1)
break;
switch (c)
{
case 0:
/* If this option set a flag, do nothing else now. */
if (long_options[option_index].flag != 0)
break;
printf ("option %s", long_options[option_index].name);
if (optarg)
printf (" with arg %s", optarg);
printf ("\n");
break;
case 'h':
print_usage(stdout, EXIT_SUCCESS);
case 'r':
recombination_flag = 1;
break;
case 'x':
extensive_search_flag = 1;
break;
case 'f':
memcpy(original_multi_fasta_filename, optarg, size_of_string(optarg) +1);
break;
case 'v':
memcpy(vcf_filename, optarg, size_of_string(optarg) +1);
break;
case 'm':
min_snps = atoi(optarg);
break;
case 'a':
window_min = atoi(optarg);
break;
case 'b':
window_max = atoi(optarg);
break;
case 'p':
uncorrected_p_value = atof(optarg);
break;
case 'i':
trimming_ratio = atof(optarg);
break;
case 't':
memcpy(tree_filename, optarg, size_of_string(optarg) +1);
break;
case '?':
/* getopt_long already printed an error message. */
break;
default:
abort ();
}
}
/* Print any remaining command line arguments (not options). */
if (optind < argc)
{
memcpy(multi_fasta_filename, argv[optind], size_of_string(argv[optind]) +1);
optind++;
}
check_file_exists_or_exit(multi_fasta_filename);
if(recombination_flag == 1)
{
check_file_exists_or_exit(vcf_filename);
check_file_exists_or_exit(tree_filename);
check_file_exists_or_exit(original_multi_fasta_filename);
run_gubbins(vcf_filename,
tree_filename,
multi_fasta_filename,
min_snps,
original_multi_fasta_filename,
window_min,
window_max,
uncorrected_p_value,
trimming_ratio,
extensive_search_flag);
}
else
{
generate_snp_sites(multi_fasta_filename, 0, ".gaps");
generate_snp_sites(multi_fasta_filename, 1, "");
}
exit(EXIT_SUCCESS);
}
