--- /dev/null
+/*
+ * fasta.c
+ *
+ * $Id$
+ *
+ *****************************************************************************
+ *
+ * Copyright (c) 2004, Luke Sheneman
+ * All rights reserved.
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ *
+ * + Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * + Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions and the following disclaimer in
+ * the documentation and/or other materials provided with the
+ * distribution.
+ * + The names of its contributors may not be used to endorse or promote
+ * products derived from this software without specific prior
+ * written permission.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+ * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+ * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+ * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+ * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+ * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+ * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+ * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+ * POSSIBILITY OF SUCH DAMAGE.
+ *
+ *****************************************************************************
+ *
+ * Functions for parsing FASTA formatted alignment files
+ *
+ *****************************************************************************
+ *
+ * AUTHOR:
+ *
+ * Luke Sheneman
+ * sheneman@cs.uidaho.edu
+ *
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <ctype.h>
+
+#include "clearcut.h"
+#include "common.h"
+#include "fasta.h"
+
+
+#define NJ_NUM_DNA_AMBIGUITY_SYMS 14
+static const char NJ_dna_ambiguity_syms[NJ_NUM_DNA_AMBIGUITY_SYMS] =
+{
+ 'M', 'R', 'W', 'S', 'Y', 'K',
+ 'V', 'H', 'D', 'B', 'X', 'N',
+ '-', '.'
+};
+
+
+#define NJ_NUM_PROTEIN_AMBIGUITY_SYMS 6
+static const char NJ_protein_ambiguity_syms[NJ_NUM_PROTEIN_AMBIGUITY_SYMS] =
+{
+ 'X', 'B', 'Z', '*', '-', '.'
+};
+
+#define NJ_NUM_DNA_SYMS 5
+static const char NJ_dna_syms[NJ_NUM_DNA_SYMS] =
+{
+ 'A', 'G', 'C', 'T', 'U'
+};
+
+
+#define NJ_NUM_PROTEIN_SYMS 20
+static const char NJ_protein_syms[NJ_NUM_PROTEIN_SYMS] =
+{
+ 'A', 'R', 'N', 'D', 'C', 'Q', 'E', 'G', 'H', 'I',
+ 'L', 'K', 'M', 'F', 'P', 'S', 'T', 'W', 'Y', 'V'
+};
+
+
+
+
+
+/*
+ * NJ_is_whitespace() - Check to see if character is whitespace
+ *
+ * INPUTS:
+ * -------
+ * c -- character to check
+ *
+ * RETURNS:
+ * --------
+ * int -- 0 if not whitespace
+ * 1 if whitespace
+ */
+static inline
+int
+NJ_is_whitespace(char c) {
+ if( c == ' ' || /* space */
+ c == '\n' || /* newline */
+ c == '\r' || /* carriage-return */
+ c == '\v' || /* vertical tab */
+ c == '\f' || /* form feed */
+ c == '\t' ) { /* horizontal tab */
+ return(1);
+ } else {
+ return(0);
+ }
+}
+
+
+
+
+
+/*
+ * NJ_is_dna() -
+ *
+ * Determines if the given symbol is DNA
+ *
+ * RETURNS: 1 if DNA
+ * 0 if not DNA
+ *
+ */
+static inline
+int
+NJ_is_dna(char c) {
+
+ int i;
+ char up_c;
+
+ up_c = toupper(c);
+
+ for(i=0;i<NJ_NUM_DNA_SYMS;i++) {
+ if(up_c == NJ_dna_syms[i]) {
+ return(1);
+ }
+ }
+
+ return(0);
+}
+
+
+
+
+
+/*
+ * NJ_is_dna_ambiguity() -
+ *
+ * Determines if the given symbol is a
+ * DNA ambiguity code
+ *
+ * RETURNS: 1 if DNA Ambiguity Code
+ * 0 if not DNA Ambiguity Code
+ *
+ */
+static inline
+int
+NJ_is_dna_ambiguity(char c) {
+
+ int i;
+ char up_c;
+
+ up_c = toupper(c);
+
+ for(i=0;i<NJ_NUM_DNA_AMBIGUITY_SYMS;i++) {
+ if(up_c == NJ_dna_ambiguity_syms[i]) {
+ return(1);
+ }
+ }
+
+ return(0);
+}
+
+
+
+/*
+ * NJ_is_protein() -
+ *
+ * Determines if supplied symbol is amino acid symbol
+ *
+ * RETURNS: 1 if amino acid
+ * 0 if not amino acid
+ *
+ */
+static inline
+int
+NJ_is_protein(char c) {
+
+ int i;
+ char up_c;
+
+ up_c = toupper(c);
+
+ for(i=0;i<NJ_NUM_PROTEIN_SYMS;i++) {
+ if(up_c == NJ_protein_syms[i]) {
+ return(1);
+ }
+ }
+
+ return(0);
+}
+
+
+
+
+/*
+ * NJ_is_protein_ambiguity() -
+ *
+ * Determines if supplied symbol is amino acid ambiguity symbol
+ *
+ * RETURNS: 1 if amino acid ambiguity symbol
+ * 0 if not amino acid ambiguity symbol
+ *
+ */
+static inline
+int
+NJ_is_protein_ambiguity(char c) {
+
+ int i;
+ char up_c;
+
+ up_c = toupper(c);
+
+ for(i=0;i<NJ_NUM_PROTEIN_AMBIGUITY_SYMS;i++) {
+ if(up_c == NJ_protein_ambiguity_syms[i]) {
+ return(1);
+ }
+ }
+
+ return(0);
+}
+
+
+
+
+
+
+/*
+ * NJ_read_fasta() - A function for inputing FASTA sequences into an alignment
+ * data structure
+ *
+ *
+ * INPUTS:
+ * -------
+ * nj_args -- A pointer to a data structure containing command-line args
+ *
+ * RETURNS:
+ * --------
+ * NJ_alignment * -- A pointer to a multiple sequence alignment
+ *
+ * DESCRIPTION:
+ * ------------
+ *
+ *
+ * This function implements a state machine parser for parsing FASTA-formatted
+ * multiple sequence alignment files.
+ *
+ * Example Input:
+ *
+ * > sequence1
+ * ATAGATATAGATTAGAATAT----TATAGATAT----ATATAT-TTT-
+ * > sequence2
+ * --ATAGATA---ATATATATATTTT--GTCTCATAGT---ATATGCTT
+ * > sequence3
+ * TTATAGATA---ATATATATATTTTAAGTCTCATAGT-A-ATATGC--
+ *
+ * This function will parse alignments for DNA or protein, and will do
+ * so mindful of ambiguity codes for these kinds of sequences. All
+ * ambiguity codes are ignored by this program for the purposes of
+ * computing a distance matrix from a multiple alignment. By design,
+ * this program does not auto-detect DNA vs. Protein, and requires that
+ * the user explictly specify that on the command-line.
+ *
+ * Gaps can be represented either with the '-' or '.' characters.
+ *
+ * Newlines and other whitespace are allowed to be interspersed
+ * throughout the sequences.
+ *
+ * Taxon labels are required to be unique, and they must start with
+ * an alphabetic character (not a number, etc.). The parser will read
+ * the first token after the > character in the description line up until the
+ * first whitespace and use that for the taxon label.
+ *
+ * For example, in the line "> taxon1 is homo sapien", the taxon label will be
+ * "taxon1"
+ *
+ */
+NJ_alignment *
+NJ_read_fasta(NJ_ARGS *nj_args) {
+
+ FILE *fp = NULL;
+ char *buf = NULL;
+ char *ptr = NULL;
+ NJ_alignment *alignment = NULL;
+
+ char c;
+ int state;
+ long int index, x, seq;
+ long int i;
+ long int bufsize, nseqs = NJ_INITIAL_NSEQS;
+ int first_sequence_flag;
+
+
+
+
+ /*
+ * In this function, we implement a FASTA alignment parser which
+ * reads in an alignment character-by-character, maintaining state
+ * information which guides the parser.
+ *
+ * The program reads either DNA or Protein alignments. All title lines
+ * and sequences can be arbitrarily long. Gaps can be represented by
+ * "-" or "." characters.
+ *
+ * Ambiguity codes are also handled.
+ *
+ */
+
+ /*
+ * We can't handle reading fasta input unless the user explicity
+ * specifies the input type...just to be sure.
+ */
+ if( (!nj_args->dna_flag && !nj_args->protein_flag) ||
+ (nj_args->dna_flag && nj_args->protein_flag) ) {
+ fprintf(stderr, "Clearcut: Explicitly specify protein or DNA\n");
+ goto XIT_BAD;
+ }
+
+ /* open specified fasta file here */
+ if(nj_args->stdin_flag) {
+ fp = stdin;
+ } else {
+ fp = fopen(nj_args->infilename, "r");
+ if(!fp) {
+ fprintf(stderr, "Clearcut: Failed to open input FASTA file: %s\n", nj_args->infilename);
+ perror("Clearcut");
+ goto XIT_BAD;
+ }
+ }
+
+ /* allocate the initial buffer */
+ bufsize = NJ_INITIAL_BUFFER_SIZE;
+ buf = (char *)calloc(bufsize, sizeof(char));
+
+ /* allocate the alignment container here */
+ alignment = (NJ_alignment *)calloc(1, sizeof(NJ_alignment));
+
+ /* allocate initial title array */
+ // printf("allocating initial title array\n");
+ alignment->titles = (char **)calloc(NJ_INITIAL_NSEQS, sizeof(char *));
+
+ /* make sure that we successfully allocated memory */
+ if(!buf || !alignment || !alignment->titles) {
+ fprintf(stderr, "Clearcut: Memory allocation error in NJ_read_fasta()\n");
+ goto XIT_BAD;
+ }
+
+ /* a flag */
+ first_sequence_flag = 1;
+
+ index = 0; /* tracks the position in buffer */
+ x = 0; /* tracks the position on sequence */
+ seq = 0; /* tracks the active sequence */
+
+ /* intitial state of state machine */
+ state = NJ_FASTA_MODE_UNKNOWN;
+
+ while(1) {
+
+ /* get the next character */
+ c = fgetc(fp);
+ if(feof(fp)) {
+
+ if(state == NJ_FASTA_MODE_SEQUENCE) {
+ buf[index+1] = '\0';
+
+ /* copy buf to alignment */
+ for(i=1;i<=alignment->length;i++) {
+ alignment->data[seq*alignment->length+i-1] = buf[i];
+ }
+ }
+
+ break;
+ }
+
+ /* make sure our dynamic buffer is big enough */
+ if(index >= bufsize) {
+ bufsize *= 2;
+ buf = (char *)realloc(buf, bufsize);
+ if(!buf) {
+ fprintf(stderr, "Clearcut: Memory allocation error in NJ_read_fasta()\n");
+ goto XIT_BAD;
+ }
+ }
+
+ switch(state) {
+
+ case NJ_FASTA_MODE_UNKNOWN:
+
+ if(!NJ_is_whitespace(c)) {
+ if(c == '>') {
+ state = NJ_FASTA_MODE_TITLE;
+ buf[0] = '>';
+ } else {
+ goto XIT_BAD;
+ }
+ }
+
+ break;
+
+ case NJ_FASTA_MODE_TITLE:
+
+ if( c == '\n' ||
+ c == '\r' ) {
+
+ buf[index] = '\0';
+ state = NJ_FASTA_MODE_SEQUENCE;
+ index = 0;
+ x = -1;
+
+ /* make sure we've allocated enough space for titles and sequences */
+ if(seq == nseqs) {
+
+ // printf("realloc(). seq = %d, nseqs = %d\n", seq, nseqs);
+
+ nseqs *= 2;
+
+ alignment->titles = (char **)realloc(alignment->titles, nseqs*sizeof(char *));
+ if(!alignment->titles) {
+ fprintf(stderr, "Clearcut: Memory allocation error in NJ_read_fasta()\n");
+ goto XIT_BAD;
+ }
+
+ alignment->data = (char *)realloc(alignment->data, alignment->length*nseqs*sizeof(char));
+ if(!alignment->data) {
+ fprintf(stderr, "Clearcut: Allocation error in NJ_read_fasta()\n");
+ goto XIT_BAD;
+ }
+ }
+
+ // printf("Allocating %d bytes for title %d: %s\n", (int)strlen(buf), (int)seq, buf);
+
+ alignment->titles[seq] = (char *)calloc(strlen(buf), sizeof(char));
+ if(!alignment->titles[seq]) {
+ fprintf(stderr, "Clearcut: Memory allocation error in NJ_read_fasta()\n");
+ goto XIT_BAD;
+ }
+
+ /* lets forward to the first non-space (space/tab) character after the '>' */
+
+ if(first_sequence_flag) {
+ ptr = buf;
+ } else {
+ ptr = &buf[1];
+ }
+ while(*ptr == '\t' || *ptr == ' ') {
+ ptr++;
+ }
+ sscanf(ptr, "%s", alignment->titles[seq]); /* get the first word and use as the title */
+
+ alignment->nseq++;
+ }
+
+ buf[index++] = c;
+
+ break;
+
+
+ case NJ_FASTA_MODE_SEQUENCE:
+
+ if(c == '>') {
+
+ if(first_sequence_flag) {
+ first_sequence_flag = 0;
+
+ /* allocate our alignment data section here */
+ alignment->length = index-1;
+
+ nseqs = NJ_INITIAL_NSEQS;
+ alignment->data = (char *)calloc(alignment->length*nseqs, sizeof(char));
+ if(!alignment->data) {
+ fprintf(stderr, "Clearcut: Allocation error in NJ_read_fasta()\n");
+ goto XIT_BAD;
+ }
+ }
+
+ if(!first_sequence_flag) {
+ if(index-1 < alignment->length) {
+ fprintf(stderr, "Clearcut: Sequences must be of uniform length in alignment at sequence %ld\n", seq);
+ goto XIT_BAD;
+ }
+ }
+
+ /* re-allocate if necessary */
+ /*
+ if(seq >= nseqs) {
+ nseqs *= 2;
+ alignment->data = (char *)realloc(alignment->data, alignment->length*nseqs*sizeof(char));
+ if(!alignment->data) {
+ fprintf(stderr, "Clearcut: Allocation error in NJ_read_fasta()\n");
+ goto XIT_BAD;
+ }
+ }
+ */
+
+ /* copy buf to alignment */
+ for(i=1;i<=alignment->length;i++) {
+ alignment->data[seq*alignment->length+i-1] = buf[i];
+ }
+
+ state = NJ_FASTA_MODE_TITLE;
+ index = 1;
+ x = 1;
+
+ buf[0] = c;
+
+ seq++;
+
+ } else {
+
+ if(NJ_is_whitespace(c)) {
+ break;
+ }
+
+ if(!first_sequence_flag) {
+ if(index-1 >= alignment->length) {
+ fprintf(stderr, "Clearcut: Sequences must be of uniform length in alignment at sequence %ld\n", seq);
+ goto XIT_BAD;
+ }
+ }
+
+
+ /*
+ * Here we check to make sure that the symbol read is appropriate
+ * for the type of data the user specified. (dna or protein).
+ * We also handle ambiguity codes by converting them to a specific
+ * assigned ambiguity code character. Ambiguity codes are ignored
+ * when computing distances
+ */
+
+ if(nj_args->dna_flag) {
+ if(NJ_is_dna(c)) {
+ buf[index++] = toupper(c);
+ } else {
+ if(NJ_is_dna_ambiguity(c)) {
+ buf[index++] = NJ_AMBIGUITY_CHAR;
+ } else {
+ fprintf(stderr, "Clearcut: Unknown symbol '%c' in nucleotide sequence %ld.\n", c, seq);
+ goto XIT_BAD;
+ }
+ }
+ } else if(nj_args->protein_flag) {
+ if(NJ_is_protein(c)) {
+ buf[index++] = toupper(c);
+ } else {
+ if(NJ_is_protein_ambiguity(c)) {
+ buf[index++] = NJ_AMBIGUITY_CHAR;
+ } else {
+ fprintf(stderr, "Clearcut: Unknown symbol '%c' in protein sequence %ld.\n", c, seq);
+ goto XIT_BAD;
+ }
+ }
+ }
+
+ }
+
+ break;
+
+ default:
+ goto XIT_BAD;
+
+ break;
+
+ }
+ }
+
+ if(index-1 != alignment->length) {
+ fprintf(stderr, "Clearcut: Sequences must be of uniform length in alignment at sequence %ld\n", seq);
+ goto XIT_BAD;
+ }
+
+ /* check for duplicate taxon labels */
+ if(!NJ_taxaname_unique(alignment)) {
+ goto XIT_BAD;
+ }
+
+ return(alignment);
+
+
+ XIT_BAD:
+
+ if(fp) {
+ fprintf(stderr, "Clearcut: Fatal error parsing FASTA file at file offset %ld.\n", ftell(fp));
+ }
+
+ if(buf) {
+ free(buf);
+ }
+
+ NJ_free_alignment(alignment);
+
+ return(NULL);
+}
+
+
+
+
+/*
+ * NJ_print_alignment() - Print multiple sequence alignment (for debugging)
+ *
+ * INPUTS:
+ * -------
+ * alignment -- A pointer to the alignment
+ *
+ * RETURNS:
+ * --------
+ * NONE
+ *
+ */
+void
+NJ_print_alignment(NJ_alignment *alignment) {
+
+ long int i, j;
+
+ printf("nseq = %ld, length = %ld\n", alignment->nseq, alignment->length);
+
+ for(i=0;i<alignment->nseq;i++) {
+
+ printf("> %s\n", alignment->titles[i]);
+
+ for(j=0;j<alignment->length;j++) {
+ printf("%c", alignment->data[i*alignment->length+j]);
+ }
+
+ printf("\n");
+ }
+
+ return;
+}
+
+
+
+
+
+
+
+/*
+ *
+ * NJ_free_alignment() - Free all of the memory allocated for the
+ * multiple sequence alignment
+ *
+ * INPUTS:
+ * -------
+ * alignment -- A pointer to the multiple sequence alignment
+ *
+ * RETURNS:
+ * --------
+ * NONE
+ *
+ */
+void
+NJ_free_alignment(NJ_alignment *alignment) {
+
+ long int i;
+
+ if(alignment) {
+
+ /* free the allocated titles */
+ if(alignment->titles) {
+ for(i=0;i<alignment->nseq;i++) {
+
+ if(alignment->titles[i]) {
+ free(alignment->titles[i]);
+ }
+ }
+
+ free(alignment->titles);
+ }
+
+ /* free the alignment data */
+ if(alignment->data) {
+ free(alignment->data);
+ }
+
+ /* free the alignment itself */
+ free(alignment);
+ }
+
+ return;
+}
+
+
+
+
+/*
+ * NJ_taxaname_unique() - Check to see if taxanames are unique in alignment
+ *
+ * INPUTS:
+ * -------
+ * alignment -- a pointer to a multiple sequence alignment
+ *
+ * OUTPUTS:
+ * --------
+ * int -- 0 if all taxanames in alignment are unique
+ * 1 if all taxanames in alignment are NOT unique
+ *
+ *
+ * DESCRIPTION:
+ * ------------
+ *
+ * Check to see if the taxanames in the alignment are unique. It
+ * will be impossible to make sense of the final tree if the taxon
+ * labels are not unqiue.
+ *
+ */
+int
+NJ_taxaname_unique(NJ_alignment *alignment) {
+
+ long int i, j;
+
+ for(i=0;i<alignment->nseq;i++) {
+ for(j=i+1;j<alignment->nseq;j++) {
+ if(!strcmp(alignment->titles[i],
+ alignment->titles[j])) {
+ fprintf(stderr, "Clearcut: Taxa %ld and %ld (%s) do not have unique taxon labels.\n",
+ i, j, alignment->titles[i]);
+ return(0);
+ }
+ }
+ }
+
+ return(1);
+}
+
+
+void
+NJ_print_titles(NJ_alignment *alignment) {
+
+ int i;
+
+ for(i=0;i<alignment->nseq;i++) {
+ printf("%d: %s\n", i, alignment->titles[i]);
+ }
+
+ return;
+}
projects / mothur.git / blobdiff