From 5a34c36f495773a8693fa49875577ea21e1e8931 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Thu, 30 Apr 2026 09:56:46 +0000 Subject: [PATCH 1/2] feat: complete rewrite implementing all improvements from plan Agent-Logs-Url: https://github.com/ChrisCreevey/catsequences/sessions/d6cfecf1-be96-4edd-82eb-e4900e7f01b7 Co-authored-by: ChrisCreevey <4115315+ChrisCreevey@users.noreply.github.com> --- README.md | 159 ++- catsequences.c | 1149 ++++++++++++++++------ tests/expected/basic_aln.fas | 10 + tests/expected/basic_occupancy.txt | 8 + tests/expected/basic_partitions.txt | 3 + tests/expected/direct_aln.fas | 10 + tests/expected/direct_partitions.txt | 3 + tests/expected/exclude.txt | 2 + tests/expected/exclude_aln.fas | 6 + tests/expected/include.txt | 2 + tests/expected/include_aln.fas | 4 + tests/expected/iqtree_partitions.txt | 3 + tests/expected/missing_n.fas | 10 + tests/expected/nexus_aln.nex | 12 + tests/expected/nexus_part_partitions.txt | 6 + tests/expected/occ_matrix.csv | 6 + tests/expected/phylip_aln.phy | 6 + tests/expected/raxml_partitions.txt | 3 + tests/run_tests.sh | 162 +++ 19 files changed, 1227 insertions(+), 337 deletions(-) create mode 100644 tests/expected/basic_aln.fas create mode 100644 tests/expected/basic_occupancy.txt create mode 100644 tests/expected/basic_partitions.txt create mode 100644 tests/expected/direct_aln.fas create mode 100644 tests/expected/direct_partitions.txt create mode 100644 tests/expected/exclude.txt create mode 100644 tests/expected/exclude_aln.fas create mode 100644 tests/expected/include.txt create mode 100644 tests/expected/include_aln.fas create mode 100644 tests/expected/iqtree_partitions.txt create mode 100644 tests/expected/missing_n.fas create mode 100644 tests/expected/nexus_aln.nex create mode 100644 tests/expected/nexus_part_partitions.txt create mode 100644 tests/expected/occ_matrix.csv create mode 100644 tests/expected/phylip_aln.phy create mode 100644 tests/expected/raxml_partitions.txt create mode 100755 tests/run_tests.sh diff --git a/README.md b/README.md index 5da0252..64bd469 100644 --- a/README.md +++ b/README.md @@ -1,91 +1,164 @@ # catsequences -### A tool for concatenating multiple fasta alignments for supermatrix phylogenetic analyses +### A tool for concatenating multiple FASTA alignments for supermatrix phylogenetic analyses Please reference this DOI: [![DOI](https://zenodo.org/badge/79135872.svg)](https://zenodo.org/badge/latestdoi/79135872) +--- -Usage: `catsequences FILENAME` +## Installation -Where: - - `FILENAME` contains a list of files with sequences to be concatenated. - - The concatenated alignments are outputted to a file called `allseqs.fas`. - - Partition information will be written to the file `allseqs.partitions.txt` +Download `catsequences.c` or clone the repository: -The sequence names are read up to the first `.` `|` or ` ` (space) or the end of the name. +```bash +git clone https://github.com/ChrisCreevey/catsequences.git +``` -**NOTE**: You can change the delimiters that are used to split the name by editing the code - see below for details +Build: +```bash +cc catsequences.c -o catsequences -lm +``` ----------------------------------- +Copy the executable to somewhere on your `PATH` (e.g. `~/bin`). -# To install: +--- -Download the file `catsequences.c` and or if you have git installed use the command: -``` -git clone https://github.com/ChrisCreevey/catsequences.git -``` +## Basic usage -To build catsequences type: ``` -cc catsequences.c -o catsequences -lm +catsequences [OPTIONS] LISTFILE +catsequences [OPTIONS] FILE1.fasta FILE2.fasta ... ``` -Copy the executable `catsequences` to somewhere on your path like `~/bin` to have access to it from anywhere. +**List-file mode** (original behaviour): --------------------------------- +```bash +ls *.fasta > list +catsequences list +``` + +**Direct FASTA mode** (pass files on the command line): -# To use: +```bash +catsequences gene1.fasta gene2.fasta gene3.fasta +``` -Usage: +Glob wildcards are also accepted inside a list file: ``` -catsequences FILENAME +genes/*.fasta ``` - where `FILENAME` is a file containing a list of files with sequences to be concatenated. +--- + +## Options + +| Flag | Description | +|------|-------------| +| `-d DELIMS` | Name-delimiter characters (default: `" \|."`) — truncate the species name at the first occurrence of any character in this string | +| `-o PREFIX` | Output filename prefix (default: `allseqs`) | +| `-m CHAR` | Missing-data character (default: `?`; common alternatives: `-`, `N`) | +| `-q` | Quiet — suppress all progress messages | +| `-v` | Verbose — print per-file scan/read progress | +| `-f FORMAT` | Partition file format: `generic` (default), `raxml`, `iqtree`, `nexus` | +| `-F FORMAT` | Alignment output format: `fasta` (default), `nexus`, `phylip` | +| `-O` | Write an occupancy matrix CSV | +| `--min-taxa N` | Skip any gene file where fewer than `N` × 100 % of species are present (e.g. `0.5` = require at least 50 % occupancy) | +| `--on-duplicate MODE` | What to do when a species name appears more than once in a file: `warn` (default — keep first, print a message), `error` (abort), `keep-last` (overwrite with the last copy) | +| `--include FILE` | Only output species whose names appear in `FILE` (one name per line, without `>`) | +| `--exclude FILE` | Exclude species whose names appear in `FILE` (one name per line, without `>`) | + +--- + +## Outputs + +| File | Description | +|------|-------------| +| `.fas` / `.nex` / `.phy` | Concatenated alignment | +| `.partitions.txt` | Partition coordinates | +| `.occupancy.txt` | Per-species and overall occupancy summary | +| `.occupancy.csv` | Occupancy matrix — present (1) / absent (0) per gene per species (written only with `-O`) | - - The concatenated alignments are outputted to a file called `allseqs.fas` - - - Partition information will be written to the file `allseqs.partitions.txt` +Species that are absent from a gene file are padded with the missing-data character for that gene's full length. - - The sequence names are read up to the first `.` `|` or ` ` (space) (or the end of the name) +--- -The easiest way to create the list input file is to use a command like: +## Partition format examples + +### generic (default) ``` -ls *.fasta > list +gene1.fasta=1-500; +gene2.fasta=501-1200; +``` + +### raxml (`-f raxml`) +``` +DNA, gene1 = 1-500 +DNA, gene2 = 501-1200 +``` + +### iqtree (`-f iqtree`) +``` +gene1 = 1-500 +gene2 = 501-1200 ``` -assuming that all your input file have the extension ".fasta" -The software will put in '?' for the entire length of any gene for which the species is missing. +### nexus (`-f nexus`) +``` +#NEXUS +begin sets; +charset gene1 = 1-500; +charset gene2 = 501-1200; +end; +``` +--- -# Trouble Shooting +## Examples -## If the total concatenated seuquences length is larger than defined in the code you recieve an ERROR saying: +Run with a list file, writing RAxML partitions and a NEXUS alignment: +```bash +catsequences -f raxml -F nexus -o mymatrix list ``` -ERROR: The total concatenated sequence length is greater than the maximum length defined (1000000)! - To fix this increase the value of MAXLEN on line 15 of 'catsequences.c` and recompile the tool. + +Require at least 50 % of species per gene, write occupancy CSV: + +```bash +catsequences --min-taxa 0.5 -O -o filtered list ``` -This line you want to change is: +Include only a subset of species: +```bash +echo -e "Homo_sapiens\nPan_troglodytes\nGorilla_gorilla" > primates.txt +catsequences --include primates.txt -o primates list ``` -#define MAXLEN 1000000 + +Change the delimiter so only `|` splits names (not `.` or space): + +```bash +catsequences -d '|' list ``` -Although this means that the total length of all the seqeunces you are trying to concatentate is greater than 1 million bp. -Please make sure that this is as expected before changing the code. +--- -## I need a different delimiter character in the sequence names than the default "|", "." or " " (space) +## Testing -Add (or remove) characters between the `" "` as needed on line 16 of `catsequences.c`: +Build and run the regression test suite: -``` -#define NAMEDELIMITER " |." +```bash +cc catsequences.c -o catsequences -lm +bash tests/run_tests.sh ``` -Then recompile (build) the tool. +--- +## Notes +- Species name lengths and the number of species and genes are all determined + dynamically at run time — there are no hard-coded limits requiring recompilation. +- The tool reads each input file twice: once to discover species and measure + sequence lengths, and once to read sequence data. This keeps peak memory + proportional to `(number of species) × (total concatenated length)`. diff --git a/catsequences.c b/catsequences.c index 1e5da56..c242591 100755 --- a/catsequences.c +++ b/catsequences.c @@ -1,313 +1,874 @@ +/* + * catsequences - concatenate multiple FASTA alignments for supermatrix + * phylogenetic analyses. + * + * Usage: + * catsequences [OPTIONS] LISTFILE + * catsequences [OPTIONS] FILE1.fasta FILE2.fasta ... + * + * The list file contains one alignment filename per line (glob wildcards + * are supported). Alternatively pass the FASTA files directly on the + * command line. + */ + #include #include #include +#include +#include +#include #ifndef TRUE -#define TRUE 1 -#endif -#ifndef FALSE -#define FALSE 0 -#endif -#ifndef TEMP -#define TEMP 2 +# define TRUE 1 +# define FALSE 0 #endif -#define MAXLEN 1000000 -#define NAMEDELIMITER " |." /* species name is initial string up to first ' ', '|' or '.' */ +/* Growing-array initial capacities */ +#define INIT_SPECIES 256 +#define INIT_FILES 256 +#define INIT_BUFLEN 4096 +/* Partition output formats */ +#define PFMT_GENERIC 0 +#define PFMT_RAXML 1 +#define PFMT_IQTREE 2 +#define PFMT_NEXUS 3 -void clean_exit(int x); +/* Alignment output formats */ +#define AFMT_FASTA 0 +#define AFMT_NEXUS 1 +#define AFMT_PHYLIP 2 -char **species_names = NULL, ***sequences = NULL; -int *seqlens = NULL, numspecies = 0, numseqs = 0, num_files=0; +/* Duplicate-handling modes */ +#define DUP_WARN 0 +#define DUP_ERROR 1 +#define DUP_KEEPLAST 2 -int main(int argc, char *argv[]) - { - FILE *infile1 = NULL, *infile2 = NULL, *outfile2 = NULL; - char *filename = NULL, string[MAXLEN], string2[MAXLEN], c = '\0'; - int i = 0, j = 0, k=0, found = 0, aln_len =0, error = FALSE, filenum=0, previous=FALSE, totlen=0; - - if(argc < 2) - { - printf("\ncatsequences is designed to concatenate sequence alignments. \n\nUsage: catsequences filename\n\tWhere filename contains a list of files with sequences to be concatenated\n\tThe concatenated alignments are outputted to a file called \"allseqs.fas\"\n\tPartition information will be written to the file \"allseqs.partitions.txt\"\n\n\n\tThe sequence names are read up to the first \".\" or \"|\" (or the end of the name)\n\n"); - exit(1); +/* ── global arrays (freed at exit) ───────────────────────────── */ +static char **species_names = NULL; +static int numspecies = 0; +static int species_cap = 0; + +static char **filelist = NULL; +static int num_files = 0; +static int files_cap = 0; + +static int *seqlens = NULL; /* actual sequence length per file */ +static int *file_active = NULL; /* 1 if file passes filters */ +static int *file_species_count = NULL; /* species found per file */ + +static char ***sequences = NULL; /* sequences[file][species] */ + +/* ── runtime options ─────────────────────────────────────────── */ +static const char *name_delimiters = " |."; +static char missing_char = '?'; +static const char *out_prefix = "allseqs"; +static int quiet = FALSE; +static int verbose = FALSE; +static int part_format = PFMT_GENERIC; +static int aln_format = AFMT_FASTA; +static int dup_mode = DUP_WARN; +static double min_taxa_frac = 0.0; +static const char *include_file = NULL; +static const char *exclude_file = NULL; +static int write_occ_csv = FALSE; + +static char **include_list = NULL; +static int include_count = 0; +static char **exclude_list = NULL; +static int exclude_count = 0; + +/* ═══════════════════════ helpers ════════════════════════════════ */ + +static void die(const char *msg) +{ + fprintf(stderr, "Error: %s\n", msg); + exit(EXIT_FAILURE); +} + +/* + * Read one line from fp into *buf (grows via realloc as needed). + * Strips the trailing \r and/or \n. + * Returns the number of characters stored (>=0), or -1 on EOF with no data. + */ +static int read_line(FILE *fp, char **buf, size_t *cap) +{ + size_t len = 0; + int c; + + if (*buf == NULL) { + *cap = INIT_BUFLEN; + *buf = malloc(*cap); + if (!*buf) die("out of memory"); + } + + for (;;) { + c = getc(fp); + if (c == EOF) { + if (len == 0) return -1; + break; + } + if (c == '\n') break; + if (c == '\r') { + int nxt = getc(fp); + if (nxt != '\n' && nxt != EOF) ungetc(nxt, fp); + break; + } + if (len + 1 >= *cap) { + *cap *= 2; + char *tmp = realloc(*buf, *cap); + if (!tmp) die("out of memory"); + *buf = tmp; + } + (*buf)[len++] = (char)c; + } + (*buf)[len] = '\0'; + return (int)len; +} + +/* Return basename of path with the extension stripped. */ +static void gene_name(const char *path, char *out, size_t outsz) +{ + const char *base = strrchr(path, '/'); + base = base ? base + 1 : path; + const char *dot = strrchr(base, '.'); + size_t len = dot ? (size_t)(dot - base) : strlen(base); + if (len >= outsz) len = outsz - 1; + memcpy(out, base, len); + out[len] = '\0'; +} + +/* Load one name per line from a file into *list / *count. */ +static void load_name_list(const char *path, char ***list, int *count) +{ + FILE *fp = fopen(path, "r"); + if (!fp) { + fprintf(stderr, "Error: Cannot open %s: %s\n", path, strerror(errno)); + exit(EXIT_FAILURE); + } + int cap = 64; + *list = malloc((size_t)cap * sizeof(char *)); + if (!*list) die("out of memory"); + *count = 0; + + char *buf = NULL; + size_t bufcap = 0; + while (read_line(fp, &buf, &bufcap) >= 0) { + if (buf[0] == '\0') continue; + if (*count >= cap) { + cap *= 2; + char **tmp = realloc(*list, (size_t)cap * sizeof(char *)); + if (!tmp) die("out of memory"); + *list = tmp; + } + (*list)[*count] = strdup(buf); + if (!(*list)[*count]) die("out of memory"); + (*count)++; + } + free(buf); + fclose(fp); +} + +/* + * Return TRUE if name (with or without leading '>') matches a list entry. + * List entries are stored without '>'. + */ +static int in_list(const char *name, char **list, int count) +{ + const char *nm = (name[0] == '>') ? name + 1 : name; + for (int i = 0; i < count; i++) + if (strcmp(nm, list[i]) == 0) return TRUE; + return FALSE; +} + +/* Append one path to the global filelist, growing as needed. */ +static void add_file(const char *path) +{ + if (num_files >= files_cap) { + files_cap = files_cap ? files_cap * 2 : INIT_FILES; + char **tmp = realloc(filelist, (size_t)files_cap * sizeof(char *)); + if (!tmp) die("out of memory"); + filelist = tmp; + } + filelist[num_files] = strdup(path); + if (!filelist[num_files]) die("out of memory"); + num_files++; +} + +/* Expand a glob pattern and add all matches; treat as literal if no match. */ +static void add_pattern(const char *pattern) +{ + glob_t g; + int r = glob(pattern, GLOB_TILDE, NULL, &g); + if (r == 0) { + for (size_t i = 0; i < g.gl_pathc; i++) + add_file(g.gl_pathv[i]); + globfree(&g); + } else { + if (r != GLOB_NOMATCH) globfree(&g); + add_file(pattern); + } +} + +/* Read a list file and expand each line as a glob pattern. */ +static void load_file_list(const char *listpath) +{ + FILE *fp = fopen(listpath, "r"); + if (!fp) { + fprintf(stderr, "Error: Cannot open list file %s: %s\n", + listpath, strerror(errno)); + exit(EXIT_FAILURE); + } + char *buf = NULL; + size_t cap = 0; + while (read_line(fp, &buf, &cap) >= 0) { + if (buf[0] == '\0') continue; + add_pattern(buf); + } + free(buf); + fclose(fp); +} + +/* Pre-flight: verify every file is readable; report all missing at once. */ +static void preflight_check(void) +{ + int bad = FALSE; + for (int i = 0; i < num_files; i++) { + FILE *fp = fopen(filelist[i], "r"); + if (!fp) { + fprintf(stderr, "Error: Cannot open %s: %s\n", + filelist[i], strerror(errno)); + bad = TRUE; + } else { + fclose(fp); + } + } + if (bad) exit(EXIT_FAILURE); +} + +/* Find species by name, or add a new entry. Returns the index. */ +static int find_or_add_species(const char *name) +{ + for (int i = 0; i < numspecies; i++) + if (strcmp(name, species_names[i]) == 0) return i; + + if (numspecies >= species_cap) { + species_cap = species_cap ? species_cap * 2 : INIT_SPECIES; + char **tmp = realloc(species_names, + (size_t)species_cap * sizeof(char *)); + if (!tmp) die("out of memory"); + species_names = tmp; + } + species_names[numspecies] = strdup(name); + if (!species_names[numspecies]) die("out of memory"); + return numspecies++; +} + +/* ═════════════════════════ first pass ═══════════════════════════ + * Scan each FASTA file to build the species list and measure + * sequence lengths. No sequences are stored yet. + * ═══════════════════════════════════════════════════════════════ */ +static void first_pass(void) +{ + seqlens = calloc((size_t)num_files, sizeof(int)); + file_active = malloc((size_t)num_files * sizeof(int)); + file_species_count = calloc((size_t)num_files, sizeof(int)); + if (!seqlens || !file_active || !file_species_count) die("out of memory"); + for (int i = 0; i < num_files; i++) file_active[i] = TRUE; + + char *buf = NULL; + size_t cap = 0; + + for (int j = 0; j < num_files; j++) { + FILE *fp = fopen(filelist[j], "r"); + if (!fp) { + fprintf(stderr, "Error: Cannot open %s\n", filelist[j]); + exit(EXIT_FAILURE); + } + + int aln_len = -1; + int in_seq = FALSE; + int seq_len = 0; + + while (read_line(fp, &buf, &cap) >= 0) { + if (buf[0] == '>') { + /* Commit length of the sequence just finished */ + if (in_seq) { + if (aln_len == -1) { + aln_len = seq_len; + seqlens[j] = aln_len; + } else if (seq_len != aln_len) { + fprintf(stderr, + "ERROR: sequences differ in length in %s" + " (%d vs %d)\n", + filelist[j], seq_len, aln_len); + exit(EXIT_FAILURE); + } + seq_len = 0; + } + in_seq = TRUE; + + /* Truncate the header at the first delimiter (after '>') */ + buf[1 + strcspn(buf + 1, name_delimiters)] = '\0'; + find_or_add_species(buf); + file_species_count[j]++; + + } else if (in_seq) { + /* Count non-whitespace characters in sequence lines */ + for (const char *p = buf; *p; p++) + if (*p != ' ' && *p != '\t') seq_len++; + } + } + + /* Commit the last sequence in the file */ + if (in_seq) { + if (aln_len == -1) { + aln_len = seq_len; + seqlens[j] = aln_len; + } else if (seq_len != aln_len) { + fprintf(stderr, + "ERROR: sequences differ in length in %s (%d vs %d)\n", + filelist[j], seq_len, aln_len); + exit(EXIT_FAILURE); + } + } + + fclose(fp); + + if (verbose) + printf("Scanned %s: %d bp, %d species\n", + filelist[j], seqlens[j], file_species_count[j]); + } + free(buf); + + /* Apply --min-taxa filter now that we know the total species count */ + if (min_taxa_frac > 0.0 && numspecies > 0) { + for (int j = 0; j < num_files; j++) { + double frac = (double)file_species_count[j] / (double)numspecies; + if (frac < min_taxa_frac) { + file_active[j] = FALSE; + if (!quiet) + printf("Skipping %s: %.1f%% species occupancy" + " < required %.1f%%\n", + filelist[j], frac * 100.0, + min_taxa_frac * 100.0); + } } - - - /* open file with list of files to be concatenated */ - if((infile1 = fopen(argv[1], "r")) == NULL) /* check to see if the file is there */ - { /* Open the fundamental tree file */ - fprintf(stderr, "Error: Cannot open list file %s\n", argv[1]); - exit(1); - } - - outfile2 = fopen("allseqs.partitions.txt", "w"); - filename=malloc(1000*sizeof(char)); - filename[0]='\0'; - string[0] ='\0'; - string2[0] = '\0'; - - - /* assign the arry to hold the names of the species (max 10000 in this build) */ - species_names = malloc(10000*sizeof(char*)); - if(species_names == NULL) clean_exit(1); - - /* Count the number of file to be read in */ - i=0; - while(!feof(infile1)) - { - c= getc(infile1); - if(c == '\n' || c == '\r') i++; - } - rewind(infile1); - num_files = i; - - - /* assign the array to hold the length of the aligned sequences in each file */ - seqlens = malloc(num_files*sizeof(int)); - if(seqlens == NULL) clean_exit(6); - j=0; - - while(!feof(infile1) && !error) - { - fscanf(infile1, "%s\n", filename); /* Get the name of the sequene file tobe opened */ - if((infile2 = fopen(filename, "r")) == NULL) /* check to see if the file is there */ - { /* Open the fundamental tree file */ - fprintf(stderr, "Error: Cannot open alignment file %s\n", filename); - exit(1); - } - - /*infile2 = fopen(filename, "r"); */ /* open the sequence file */ - /*printf("%s\n", filename);*/ - - fprintf(outfile2, "%s\t=\t%d-", filename, totlen-j+1); - - aln_len = -1; - /* read in the names and see how many species we need for this gene */ - c = getc(infile2); - while(!feof(infile2) && !error) - { - if(c == '>') /* read in the names are */ - { - i=0; - string[i] = c; - i++; - while(!feof(infile2) && (c=getc(infile2)) != '\n' && c != '\r' && c != '.') - { - if(c != '\r' && c != '\n' && !feof(infile2)) - { - string[i] = c; - i++; - } - if(i == MAXLEN) - { - printf("ERROR: The total concatenated sequence length is greater than the maximum length defined (%d)!\n", MAXLEN); - printf("\tTo fix this increase the value of MAXLEN on line 15 of 'catsequences.c` and recompile the tool\n"); - exit(1); - } - } - if(c == '.') - { - while(!feof(infile2) && (c=getc(infile2)) != '\n' && c != '\r'); - } - string[i] = '\0'; - - - - if(string[0] == '>') - { - string[strcspn(string, NAMEDELIMITER)] = '\0'; /* Truncate the name at the first instnace of any of the NAME delimiters (to capture the species name and exclude gene names for instnace) */ - found = -1; - for(i=0; i') - { - string[strcspn(string, NAMEDELIMITER)] = '\0'; /* species name is initial string up to [ |.] */ - found = -1; - for(i=0; i 0) - { - previous=TRUE; - printf("Found an extra copy of %s in file %s.... ignoring this copy\n", species_names[i], filename); - } - i=numspecies; - } - } - if(found == -1) - { - printf("error: %s not found\n", string); - exit(EXIT_FAILURE); - } - } - else - { - if(!previous) - { - strcat(sequences[filenum][found], string); - if(strlen(sequences[filenum][found]) > seqlens[filenum]) printf("%lu/%d\n", strlen(sequences[filenum][found]), seqlens[filenum]); - } - } - - - i=0; - } - fclose(infile2); - filenum++; + } +} + +/* ════════════════════════ second pass ═══════════════════════════ + * Read the actual sequence data into memory. + * ═══════════════════════════════════════════════════════════════ */ +static void second_pass(void) +{ + sequences = calloc((size_t)num_files, sizeof(char **)); + if (!sequences) die("out of memory"); + + for (int j = 0; j < num_files; j++) { + sequences[j] = calloc((size_t)numspecies, sizeof(char *)); + if (!sequences[j]) die("out of memory"); + for (int k = 0; k < numspecies; k++) { + sequences[j][k] = malloc((size_t)(seqlens[j] + 1)); + if (!sequences[j][k]) die("out of memory"); + sequences[j][k][0] = '\0'; + } + } + + char *buf = NULL; + size_t cap = 0; + + for (int j = 0; j < num_files; j++) { + FILE *fp = fopen(filelist[j], "r"); + if (!fp) { + fprintf(stderr, "Error: Cannot reopen %s\n", filelist[j]); + exit(EXIT_FAILURE); } - fclose(infile1); - - - -/* for(i=0; i= 0) { + if (buf[0] == '>') { + buf[1 + strcspn(buf + 1, name_delimiters)] = '\0'; + + found = -1; + previous = FALSE; + for (int i = 0; i < numspecies; i++) { + if (strcmp(buf, species_names[i]) == 0) { + found = i; + break; + } + } + if (found < 0) { + fprintf(stderr, + "Error: species '%s' not in master list\n", buf); + exit(EXIT_FAILURE); } - else - { -/* printf("-%d\t", seqlens[j]-1); */ - for(k=0; k<(seqlens[j]-1); k++) - { - fprintf(infile1, "?"); + + /* Handle duplicates */ + if (sequences[j][found][0] != '\0') { + switch (dup_mode) { + case DUP_WARN: + printf("Duplicate %s in %s: ignoring extra copy\n", + species_names[found] + 1, filelist[j]); + previous = TRUE; + break; + case DUP_ERROR: + fprintf(stderr, + "Error: duplicate species %s in %s\n", + species_names[found] + 1, filelist[j]); + exit(EXIT_FAILURE); + case DUP_KEEPLAST: + sequences[j][found][0] = '\0'; + break; } } - - + + } else if (found >= 0 && !previous) { + /* Append sequence line, guarding against overflow */ + size_t have = strlen(sequences[j][found]); + size_t room = (size_t)seqlens[j] - have; + strncat(sequences[j][found], buf, room); } - fprintf(infile1, "\n"); - /* printf("\n"); */ } - fclose(infile1); - free(filename); - printf("Finished! Cancatenated alignment written to \"allseqs.fas\"\n\n"); - } - -void clean_exit(int x) - { - printf("Error: out of memory at %d\n", x); - exit(EXIT_FAILURE); + fclose(fp); + + if (verbose) + printf("Read %s\n", filelist[j]); + } + free(buf); +} + +/* ════════════════════════ output helpers ════════════════════════ */ + +static int count_active(void) +{ + int n = 0; + for (int j = 0; j < num_files; j++) if (file_active[j]) n++; + return n; +} + +static int total_aln_len(void) +{ + int len = 0; + for (int j = 0; j < num_files; j++) + if (file_active[j]) len += seqlens[j]; + return len; +} + +/* Return TRUE if species i should appear in the output. */ +static int species_shown(int i) +{ + if (include_count > 0 && + !in_list(species_names[i], include_list, include_count)) + return FALSE; + if (exclude_count > 0 && + in_list(species_names[i], exclude_list, exclude_count)) + return FALSE; + return TRUE; +} + +static int count_shown_species(void) +{ + int n = 0; + for (int i = 0; i < numspecies; i++) if (species_shown(i)) n++; + return n; +} + +/* Longest displayed name (without '>') among shown species. */ +static int max_displayed_name_len(void) +{ + int mx = 0; + for (int i = 0; i < numspecies; i++) { + if (!species_shown(i)) continue; + int l = (int)strlen(species_names[i]) - 1; /* skip '>' */ + if (l > mx) mx = l; + } + return mx; +} + +/* ════════════════════════ write partitions ══════════════════════ */ +static void write_partitions(void) +{ + char path[4096]; + snprintf(path, sizeof(path), "%s.partitions.txt", out_prefix); + FILE *fp = fopen(path, "w"); + if (!fp) { + fprintf(stderr, "Error: Cannot create %s\n", path); + exit(EXIT_FAILURE); + } + + int totlen = 0; + char gname[512]; + + if (part_format == PFMT_NEXUS) + fprintf(fp, "#NEXUS\nbegin sets;\n"); + + for (int j = 0; j < num_files; j++) { + if (!file_active[j]) continue; + int start = totlen + 1; + totlen += seqlens[j]; + gene_name(filelist[j], gname, sizeof(gname)); + + switch (part_format) { + case PFMT_GENERIC: + fprintf(fp, "%s\t=\t%d-%d;\n", filelist[j], start, totlen); + break; + case PFMT_RAXML: + fprintf(fp, "DNA, %s = %d-%d\n", gname, start, totlen); + break; + case PFMT_IQTREE: + fprintf(fp, "%s = %d-%d\n", gname, start, totlen); + break; + case PFMT_NEXUS: + fprintf(fp, "\tcharset %s = %d-%d;\n", gname, start, totlen); + break; + } + } + + if (part_format == PFMT_NEXUS) fprintf(fp, "end;\n"); + fclose(fp); + if (!quiet) printf("Partitions written to \"%s\"\n", path); +} + +/* ════════════════════════ write alignment ═══════════════════════ */ +static void write_alignment(void) +{ + const char *ext = aln_format == AFMT_NEXUS ? "nex" : + aln_format == AFMT_PHYLIP ? "phy" : "fas"; + char path[4096]; + snprintf(path, sizeof(path), "%s.%s", out_prefix, ext); + FILE *fp = fopen(path, "w"); + if (!fp) { + fprintf(stderr, "Error: Cannot create %s\n", path); + exit(EXIT_FAILURE); + } + + int nsp = count_shown_species(); + int tot = total_aln_len(); + int namelen = max_displayed_name_len(); + + /* Format-specific headers */ + if (aln_format == AFMT_PHYLIP) { + fprintf(fp, " %d %d\n", nsp, tot); + } else if (aln_format == AFMT_NEXUS) { + fprintf(fp, "#NEXUS\nbegin data;\n"); + fprintf(fp, "\tdimensions ntax=%d nchar=%d;\n", nsp, tot); + fprintf(fp, "\tformat datatype=dna missing=%c gap=-;\n", missing_char); + fprintf(fp, "matrix\n"); + } + + for (int i = 0; i < numspecies; i++) { + if (!species_shown(i)) continue; + + const char *dname = species_names[i] + 1; /* skip '>' */ + + if (aln_format == AFMT_FASTA) { + fprintf(fp, "%s\n", species_names[i]); /* include '>' */ + } else if (aln_format == AFMT_PHYLIP) { + fprintf(fp, "%-*s ", namelen, dname); + } else { /* NEXUS */ + fprintf(fp, "%-*s\t", namelen, dname); + } + + for (int j = 0; j < num_files; j++) { + if (!file_active[j]) continue; + if (sequences[j][i][0] != '\0') { + fprintf(fp, "%s", sequences[j][i]); + } else { + for (int k = 0; k < seqlens[j]; k++) + fputc(missing_char, fp); + } + } + fputc('\n', fp); + } + + if (aln_format == AFMT_NEXUS) fprintf(fp, ";\nend;\n"); + fclose(fp); + if (!quiet) printf("Alignment written to \"%s\"\n", path); +} + +/* ════════════════════════ write occupancy ═══════════════════════ */ +static void write_occupancy(void) +{ + char path[4096]; + snprintf(path, sizeof(path), "%s.occupancy.txt", out_prefix); + FILE *fp = fopen(path, "w"); + if (!fp) { + fprintf(stderr, "Warning: cannot create %s\n", path); + return; + } + + int active = count_active(); + int tot = total_aln_len(); + long total_cells = 0; + long present_cells = 0; + + fprintf(fp, "Species\tPresent_genes\tMissing_chars\tMissing_pct\n"); + + for (int i = 0; i < numspecies; i++) { + if (!species_shown(i)) continue; + int present = 0, miss_chars = 0; + for (int j = 0; j < num_files; j++) { + if (!file_active[j]) continue; + total_cells++; + if (sequences[j][i][0] != '\0') { + present++; + present_cells++; + } else { + miss_chars += seqlens[j]; + } + } + double pct = tot > 0 ? 100.0 * miss_chars / tot : 0.0; + fprintf(fp, "%s\t%d/%d\t%d\t%.1f%%\n", + species_names[i] + 1, present, active, miss_chars, pct); + } + + if (total_cells > 0) + fprintf(fp, + "\nMatrix completeness: %.1f%%" + " (%ld/%ld gene-species cells present)\n", + 100.0 * (double)present_cells / (double)total_cells, + present_cells, total_cells); + + fclose(fp); + if (!quiet) printf("Occupancy summary written to \"%s\"\n", path); +} + +static void write_occupancy_csv(void) +{ + char path[4096]; + snprintf(path, sizeof(path), "%s.occupancy.csv", out_prefix); + FILE *fp = fopen(path, "w"); + if (!fp) { + fprintf(stderr, "Warning: cannot create %s\n", path); + return; + } + + char gname[512]; + fprintf(fp, "species"); + for (int j = 0; j < num_files; j++) { + if (!file_active[j]) continue; + gene_name(filelist[j], gname, sizeof(gname)); + fprintf(fp, ",%s", gname); } + fputc('\n', fp); + + for (int i = 0; i < numspecies; i++) { + if (!species_shown(i)) continue; + fprintf(fp, "%s", species_names[i] + 1); + for (int j = 0; j < num_files; j++) { + if (!file_active[j]) continue; + fprintf(fp, ",%d", sequences[j][i][0] != '\0' ? 1 : 0); + } + fputc('\n', fp); + } + + fclose(fp); + if (!quiet) printf("Occupancy CSV written to \"%s\"\n", path); +} + +/* ════════════════════════ cleanup ════════════════════════════════ */ +static void free_all(void) +{ + if (species_names) { + for (int i = 0; i < numspecies; i++) free(species_names[i]); + free(species_names); + } + if (sequences) { + for (int j = 0; j < num_files; j++) { + if (sequences[j]) { + for (int k = 0; k < numspecies; k++) + free(sequences[j][k]); + free(sequences[j]); + } + } + free(sequences); + } + if (filelist) { + for (int i = 0; i < num_files; i++) free(filelist[i]); + free(filelist); + } + free(seqlens); + free(file_active); + free(file_species_count); + if (include_list) { + for (int i = 0; i < include_count; i++) free(include_list[i]); + free(include_list); + } + if (exclude_list) { + for (int i = 0; i < exclude_count; i++) free(exclude_list[i]); + free(exclude_list); + } +} + +/* ════════════════════════ usage ══════════════════════════════════ */ +static void usage(const char *prog) +{ + printf( + "\ncatsequences - concatenate FASTA alignments for supermatrix" + " phylogenetics\n" + "\nUsage:\n" + " %s [OPTIONS] LISTFILE\n" + " %s [OPTIONS] FILE1.fasta FILE2.fasta ...\n" + "\nOptions:\n" + " -d DELIMS Name delimiter characters" + " (default: \" |.\")\n" + " -o PREFIX Output filename prefix (default: allseqs)\n" + " -m CHAR Missing-data character (default: ?)\n" + " -q Quiet: suppress progress output\n" + " -v Verbose: print per-file progress\n" + " -f FORMAT Partition format:\n" + " generic (default) | raxml | iqtree | nexus\n" + " -F FORMAT Alignment format:\n" + " fasta (default) | nexus | phylip\n" + " -O Write occupancy matrix CSV\n" + " --min-taxa N Minimum fraction of species required per" + " gene [0.0-1.0]\n" + " --on-duplicate M Duplicate-entry handling:\n" + " warn (default) | error | keep-last\n" + " --include FILE Only include species listed in FILE\n" + " --exclude FILE Exclude species listed in FILE\n" + "\nOutputs:\n" + " .fas / .nex / .phy Concatenated alignment\n" + " .partitions.txt Partition information\n" + " .occupancy.txt Per-species occupancy summary\n" + " .occupancy.csv Occupancy matrix (with -O)\n" + "\n", + prog, prog); +} + +/* ═════════════════════════ main ══════════════════════════════════ */ +int main(int argc, char *argv[]) +{ + static const struct option long_opts[] = { + {"min-taxa", required_argument, NULL, 1}, + {"on-duplicate", required_argument, NULL, 2}, + {"include", required_argument, NULL, 3}, + {"exclude", required_argument, NULL, 4}, + {"help", no_argument, NULL, 'h'}, + {NULL, 0, NULL, 0} + }; + + int opt; + while ((opt = getopt_long(argc, argv, "d:o:m:qvf:F:Oh", + long_opts, NULL)) != -1) { + switch (opt) { + case 'd': name_delimiters = optarg; break; + case 'o': out_prefix = optarg; break; + case 'm': missing_char = optarg[0]; break; + case 'q': quiet = TRUE; break; + case 'v': verbose = TRUE; break; + case 'O': write_occ_csv = TRUE; break; + case 'f': + if (strcmp(optarg, "generic") == 0) part_format = PFMT_GENERIC; + else if (strcmp(optarg, "raxml") == 0) part_format = PFMT_RAXML; + else if (strcmp(optarg, "iqtree") == 0) part_format = PFMT_IQTREE; + else if (strcmp(optarg, "nexus") == 0) part_format = PFMT_NEXUS; + else { + fprintf(stderr, "Unknown partition format: %s\n", optarg); + exit(EXIT_FAILURE); + } + break; + case 'F': + if (strcmp(optarg, "fasta") == 0) aln_format = AFMT_FASTA; + else if (strcmp(optarg, "nexus") == 0) aln_format = AFMT_NEXUS; + else if (strcmp(optarg, "phylip") == 0) aln_format = AFMT_PHYLIP; + else { + fprintf(stderr, "Unknown alignment format: %s\n", optarg); + exit(EXIT_FAILURE); + } + break; + case 1: /* --min-taxa */ + min_taxa_frac = atof(optarg); + if (min_taxa_frac < 0.0 || min_taxa_frac > 1.0) { + fprintf(stderr, + "Error: --min-taxa must be in [0.0, 1.0]\n"); + exit(EXIT_FAILURE); + } + break; + case 2: /* --on-duplicate */ + if (strcmp(optarg, "warn") == 0) dup_mode = DUP_WARN; + else if (strcmp(optarg, "error") == 0) dup_mode = DUP_ERROR; + else if (strcmp(optarg, "keep-last") == 0) dup_mode = DUP_KEEPLAST; + else { + fprintf(stderr, "Unknown --on-duplicate mode: %s\n", optarg); + exit(EXIT_FAILURE); + } + break; + case 3: include_file = optarg; break; + case 4: exclude_file = optarg; break; + case 'h': usage(argv[0]); return EXIT_SUCCESS; + default: usage(argv[0]); exit(EXIT_FAILURE); + } + } + + if (optind >= argc) { + usage(argv[0]); + exit(EXIT_FAILURE); + } + + /* Load include / exclude species lists if provided */ + if (include_file) + load_name_list(include_file, &include_list, &include_count); + if (exclude_file) + load_name_list(exclude_file, &exclude_list, &exclude_count); + + /* + * Determine whether the remaining argument(s) are a list file or + * direct FASTA files. A single argument whose first byte is '>' + * is treated as a FASTA file; otherwise it is a list file. + * Multiple arguments are always treated as direct FASTA files. + */ + if (argc - optind == 1) { + FILE *probe = fopen(argv[optind], "r"); + if (!probe) { + fprintf(stderr, "Error: Cannot open %s: %s\n", + argv[optind], strerror(errno)); + exit(EXIT_FAILURE); + } + int first = getc(probe); + fclose(probe); + if (first == '>') { + add_file(argv[optind]); /* single FASTA file */ + } else { + load_file_list(argv[optind]); /* list file */ + } + } else { + /* Multiple positional arguments → direct FASTA files */ + for (int i = optind; i < argc; i++) + add_pattern(argv[i]); + } + + if (num_files == 0) { + fprintf(stderr, "Error: No input files specified\n"); + exit(EXIT_FAILURE); + } + + preflight_check(); + + if (!quiet) + printf("catsequences: processing %d file(s)\n", num_files); + + first_pass(); + + if (!quiet) + printf("Found %d species across %d file(s)\n", numspecies, num_files); + + second_pass(); + + write_partitions(); + write_alignment(); + write_occupancy(); + if (write_occ_csv) write_occupancy_csv(); + + if (!quiet) { + printf("Total concatenated alignment length: %d bp\n", + total_aln_len()); + printf("Done.\n"); + } + + free_all(); + return EXIT_SUCCESS; +} diff --git a/tests/expected/basic_aln.fas b/tests/expected/basic_aln.fas new file mode 100644 index 0000000..623903a --- /dev/null +++ b/tests/expected/basic_aln.fas @@ -0,0 +1,10 @@ +>species1 +123456789abcdefghijklmno??????????? +>species2 +1234-6789abcd???????????pqrstuvwxyz +>species3 +1234-----a-b-e----jklmno??????????? +>species4 +?????????????efghiiiiiiipqrtttttttt +>species5 +????????????????????????dfgknwognso diff --git a/tests/expected/basic_occupancy.txt b/tests/expected/basic_occupancy.txt new file mode 100644 index 0000000..f7d8e69 --- /dev/null +++ b/tests/expected/basic_occupancy.txt @@ -0,0 +1,8 @@ +Species Present_genes Missing_chars Missing_pct +species1 2/3 11 31.4% +species2 2/3 11 31.4% +species3 2/3 11 31.4% +species4 2/3 13 37.1% +species5 1/3 24 68.6% + +Matrix completeness: 60.0% (9/15 gene-species cells present) diff --git a/tests/expected/basic_partitions.txt b/tests/expected/basic_partitions.txt new file mode 100644 index 0000000..5dd6b18 --- /dev/null +++ b/tests/expected/basic_partitions.txt @@ -0,0 +1,3 @@ +test1.txt = 1-13; +test2.txt = 14-24; +test3.txt = 25-35; diff --git a/tests/expected/direct_aln.fas b/tests/expected/direct_aln.fas new file mode 100644 index 0000000..623903a --- /dev/null +++ b/tests/expected/direct_aln.fas @@ -0,0 +1,10 @@ +>species1 +123456789abcdefghijklmno??????????? +>species2 +1234-6789abcd???????????pqrstuvwxyz +>species3 +1234-----a-b-e----jklmno??????????? +>species4 +?????????????efghiiiiiiipqrtttttttt +>species5 +????????????????????????dfgknwognso diff --git a/tests/expected/direct_partitions.txt b/tests/expected/direct_partitions.txt new file mode 100644 index 0000000..5dd6b18 --- /dev/null +++ b/tests/expected/direct_partitions.txt @@ -0,0 +1,3 @@ +test1.txt = 1-13; +test2.txt = 14-24; +test3.txt = 25-35; diff --git a/tests/expected/exclude.txt b/tests/expected/exclude.txt new file mode 100644 index 0000000..fb18fc5 --- /dev/null +++ b/tests/expected/exclude.txt @@ -0,0 +1,2 @@ +species4 +species5 diff --git a/tests/expected/exclude_aln.fas b/tests/expected/exclude_aln.fas new file mode 100644 index 0000000..46eb20f --- /dev/null +++ b/tests/expected/exclude_aln.fas @@ -0,0 +1,6 @@ +>species1 +123456789abcdefghijklmno??????????? +>species2 +1234-6789abcd???????????pqrstuvwxyz +>species3 +1234-----a-b-e----jklmno??????????? diff --git a/tests/expected/include.txt b/tests/expected/include.txt new file mode 100644 index 0000000..ee1c6b3 --- /dev/null +++ b/tests/expected/include.txt @@ -0,0 +1,2 @@ +species1 +species3 diff --git a/tests/expected/include_aln.fas b/tests/expected/include_aln.fas new file mode 100644 index 0000000..69580d1 --- /dev/null +++ b/tests/expected/include_aln.fas @@ -0,0 +1,4 @@ +>species1 +123456789abcdefghijklmno??????????? +>species3 +1234-----a-b-e----jklmno??????????? diff --git a/tests/expected/iqtree_partitions.txt b/tests/expected/iqtree_partitions.txt new file mode 100644 index 0000000..02b3cb3 --- /dev/null +++ b/tests/expected/iqtree_partitions.txt @@ -0,0 +1,3 @@ +test1 = 1-13 +test2 = 14-24 +test3 = 25-35 diff --git a/tests/expected/missing_n.fas b/tests/expected/missing_n.fas new file mode 100644 index 0000000..a761a18 --- /dev/null +++ b/tests/expected/missing_n.fas @@ -0,0 +1,10 @@ +>species1 +123456789abcdefghijklmnoNNNNNNNNNNN +>species2 +1234-6789abcdNNNNNNNNNNNpqrstuvwxyz +>species3 +1234-----a-b-e----jklmnoNNNNNNNNNNN +>species4 +NNNNNNNNNNNNNefghiiiiiiipqrtttttttt +>species5 +NNNNNNNNNNNNNNNNNNNNNNNNdfgknwognso diff --git a/tests/expected/nexus_aln.nex b/tests/expected/nexus_aln.nex new file mode 100644 index 0000000..3cafd10 --- /dev/null +++ b/tests/expected/nexus_aln.nex @@ -0,0 +1,12 @@ +#NEXUS +begin data; + dimensions ntax=5 nchar=35; + format datatype=dna missing=? gap=-; +matrix +species1 123456789abcdefghijklmno??????????? +species2 1234-6789abcd???????????pqrstuvwxyz +species3 1234-----a-b-e----jklmno??????????? +species4 ?????????????efghiiiiiiipqrtttttttt +species5 ????????????????????????dfgknwognso +; +end; diff --git a/tests/expected/nexus_part_partitions.txt b/tests/expected/nexus_part_partitions.txt new file mode 100644 index 0000000..eec18a9 --- /dev/null +++ b/tests/expected/nexus_part_partitions.txt @@ -0,0 +1,6 @@ +#NEXUS +begin sets; + charset test1 = 1-13; + charset test2 = 14-24; + charset test3 = 25-35; +end; diff --git a/tests/expected/occ_matrix.csv b/tests/expected/occ_matrix.csv new file mode 100644 index 0000000..979ee56 --- /dev/null +++ b/tests/expected/occ_matrix.csv @@ -0,0 +1,6 @@ +species,test1,test2,test3 +species1,1,1,0 +species2,1,0,1 +species3,1,1,0 +species4,0,1,1 +species5,0,0,1 diff --git a/tests/expected/phylip_aln.phy b/tests/expected/phylip_aln.phy new file mode 100644 index 0000000..184605c --- /dev/null +++ b/tests/expected/phylip_aln.phy @@ -0,0 +1,6 @@ + 5 35 +species1 123456789abcdefghijklmno??????????? +species2 1234-6789abcd???????????pqrstuvwxyz +species3 1234-----a-b-e----jklmno??????????? +species4 ?????????????efghiiiiiiipqrtttttttt +species5 ????????????????????????dfgknwognso diff --git a/tests/expected/raxml_partitions.txt b/tests/expected/raxml_partitions.txt new file mode 100644 index 0000000..6b15e71 --- /dev/null +++ b/tests/expected/raxml_partitions.txt @@ -0,0 +1,3 @@ +DNA, test1 = 1-13 +DNA, test2 = 14-24 +DNA, test3 = 25-35 diff --git a/tests/run_tests.sh b/tests/run_tests.sh new file mode 100755 index 0000000..2f84b94 --- /dev/null +++ b/tests/run_tests.sh @@ -0,0 +1,162 @@ +#!/usr/bin/env bash +# tests/run_tests.sh - regression test suite for catsequences +# +# Run from the repository root: +# cc catsequences.c -o catsequences -lm +# bash tests/run_tests.sh +# +# Exits 0 on success, non-zero if any test fails. + +set -euo pipefail + +REPO="$(cd "$(dirname "$0")/.." && pwd)" +BIN="${REPO}/catsequences" +EXPECTED="${REPO}/tests/expected" +TMPDIR_RUN="$(mktemp -d)" +trap 'rm -rf "$TMPDIR_RUN"' EXIT + +PASS=0 +FAIL=0 + +# ── helper: compare two files, report result ───────────────────── +check() { + local label="$1" got="$2" want="$3" + if diff -q "$got" "$want" > /dev/null 2>&1; then + echo "PASS: $label" + PASS=$(( PASS + 1 )) + else + echo "FAIL: $label" + diff "$want" "$got" | head -20 + FAIL=$(( FAIL + 1 )) + fi +} + +# ── test 1: basic concatenation via list file ───────────────────── +OUT="${TMPDIR_RUN}/basic" +"$BIN" -q -o "$OUT" "${REPO}/list" +check "basic FASTA output" "${OUT}.fas" "${EXPECTED}/basic_aln.fas" +check "basic partitions" "${OUT}.partitions.txt" "${EXPECTED}/basic_partitions.txt" +check "basic occupancy summary" "${OUT}.occupancy.txt" "${EXPECTED}/basic_occupancy.txt" + +# ── test 2: direct FASTA files as positional arguments ─────────── +# Run from REPO so relative filenames match the expected partition output +OUT="${TMPDIR_RUN}/direct" +( cd "$REPO" && "$BIN" -q -o "$OUT" test1.txt test2.txt test3.txt ) +check "direct-args FASTA" "${OUT}.fas" "${EXPECTED}/direct_aln.fas" +check "direct-args partitions" "${OUT}.partitions.txt" "${EXPECTED}/direct_partitions.txt" + +# ── test 3: partition formats ───────────────────────────────────── +OUT="${TMPDIR_RUN}/raxml" +"$BIN" -q -f raxml -o "$OUT" "${REPO}/list" +check "RAxML partition format" "${OUT}.partitions.txt" "${EXPECTED}/raxml_partitions.txt" + +OUT="${TMPDIR_RUN}/iqtree" +"$BIN" -q -f iqtree -o "$OUT" "${REPO}/list" +check "IQ-TREE partition format" "${OUT}.partitions.txt" "${EXPECTED}/iqtree_partitions.txt" + +OUT="${TMPDIR_RUN}/nexpart" +"$BIN" -q -f nexus -o "$OUT" "${REPO}/list" +check "NEXUS partition format" "${OUT}.partitions.txt" "${EXPECTED}/nexus_part_partitions.txt" + +# ── test 4: alignment output formats ───────────────────────────── +OUT="${TMPDIR_RUN}/nexus" +"$BIN" -q -F nexus -o "$OUT" "${REPO}/list" +check "NEXUS alignment format" "${OUT}.nex" "${EXPECTED}/nexus_aln.nex" + +OUT="${TMPDIR_RUN}/phylip" +"$BIN" -q -F phylip -o "$OUT" "${REPO}/list" +check "PHYLIP alignment format" "${OUT}.phy" "${EXPECTED}/phylip_aln.phy" + +# ── test 5: occupancy CSV ───────────────────────────────────────── +OUT="${TMPDIR_RUN}/occ" +"$BIN" -q -O -o "$OUT" "${REPO}/list" +check "occupancy matrix CSV" "${OUT}.occupancy.csv" "${EXPECTED}/occ_matrix.csv" + +# ── test 6: missing-data character ──────────────────────────────── +OUT="${TMPDIR_RUN}/missn" +"$BIN" -q -m N -o "$OUT" "${REPO}/list" +check "missing char N" "${OUT}.fas" "${EXPECTED}/missing_n.fas" + +# ── test 7: species include filter ─────────────────────────────── +OUT="${TMPDIR_RUN}/inc" +"$BIN" -q --include "${EXPECTED}/include.txt" -o "$OUT" "${REPO}/list" +check "species include filter" "${OUT}.fas" "${EXPECTED}/include_aln.fas" + +# ── test 8: species exclude filter ─────────────────────────────── +OUT="${TMPDIR_RUN}/exc" +"$BIN" -q --exclude "${EXPECTED}/exclude.txt" -o "$OUT" "${REPO}/list" +check "species exclude filter" "${OUT}.fas" "${EXPECTED}/exclude_aln.fas" + +# ── test 9: pre-flight catches all missing files ────────────────── +printf 'test1.txt\nmissing_a.txt\ntest3.txt\nmissing_b.txt\n' \ + > "${TMPDIR_RUN}/badlist.txt" +ERROUT=$( cd "$REPO" && "$BIN" -o "${TMPDIR_RUN}/junk" \ + "${TMPDIR_RUN}/badlist.txt" 2>&1 || true ) +NMISSING=$(echo "$ERROUT" | grep -c "Cannot open" || true) +if [ "$NMISSING" -eq 2 ]; then + echo "PASS: pre-flight reports all missing files" + PASS=$(( PASS + 1 )) +else + echo "FAIL: pre-flight should report 2 missing files, got: $NMISSING" + echo "$ERROUT" + FAIL=$(( FAIL + 1 )) +fi + +# ── test 10: duplicate handling – warn (default) ────────────────── +printf '>sp1\nAAAA\n>sp1\nTTTT\n>sp2\nCCCC\n' > "${TMPDIR_RUN}/dup.fasta" +OUT="${TMPDIR_RUN}/dup_warn" +STDOUT=$("$BIN" -q -o "$OUT" "${TMPDIR_RUN}/dup.fasta" 2>&1) +if echo "$STDOUT" | grep -q "Duplicate"; then + echo "PASS: duplicate warn message printed" + PASS=$(( PASS + 1 )) +else + echo "FAIL: expected 'Duplicate' warning, got: $STDOUT" + FAIL=$(( FAIL + 1 )) +fi + +# ── test 11: duplicate handling – error ─────────────────────────── +if ! "$BIN" -q --on-duplicate error \ + -o "${TMPDIR_RUN}/dup_err" "${TMPDIR_RUN}/dup.fasta" > /dev/null 2>&1; then + echo "PASS: --on-duplicate error exits non-zero" + PASS=$(( PASS + 1 )) +else + echo "FAIL: --on-duplicate error should exit with failure" + FAIL=$(( FAIL + 1 )) +fi + +# ── test 12: duplicate handling – keep-last ─────────────────────── +OUT="${TMPDIR_RUN}/dup_last" +"$BIN" -q --on-duplicate keep-last -o "$OUT" "${TMPDIR_RUN}/dup.fasta" +SP1SEQ=$(grep -A1 '^>sp1$' "${OUT}.fas" | tail -1) +if [ "$SP1SEQ" = "TTTT" ]; then + echo "PASS: --on-duplicate keep-last keeps last sequence" + PASS=$(( PASS + 1 )) +else + echo "FAIL: keep-last expected 'TTTT', got '$SP1SEQ'" + FAIL=$(( FAIL + 1 )) +fi + +# ── test 13: glob expansion in list file ───────────────────────── +# Create a list file using a glob pattern and verify results match basic run +printf '%s/test*.txt\n' "$REPO" > "${TMPDIR_RUN}/globlist.txt" +OUT="${TMPDIR_RUN}/glob" +"$BIN" -q -o "$OUT" "${TMPDIR_RUN}/globlist.txt" +# The glob expands test1.txt test2.txt test3.txt in sorted order – same as basic +check "glob expansion in list file" "${OUT}.fas" "${EXPECTED}/basic_aln.fas" + +# ── test 14: output prefix (-o) and missing char ───────────────── +OUT="${TMPDIR_RUN}/prefix_test/myprefix" +mkdir -p "${TMPDIR_RUN}/prefix_test" +"$BIN" -q -o "$OUT" "${REPO}/list" +if [ -f "${OUT}.fas" ] && [ -f "${OUT}.partitions.txt" ]; then + echo "PASS: -o prefix creates correctly named files" + PASS=$(( PASS + 1 )) +else + echo "FAIL: -o prefix did not create expected files" + FAIL=$(( FAIL + 1 )) +fi + +# ── summary ─────────────────────────────────────────────────────── +echo "" +echo "Results: ${PASS} passed, ${FAIL} failed" +[ "$FAIL" -eq 0 ] From a71bafcb6521dda51054eb2cdeb80dc153ef64bc Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Thu, 30 Apr 2026 10:00:06 +0000 Subject: [PATCH 2/2] fix: address code review feedback (realloc safety, offset tracking, PHYLIP header, diff order) Agent-Logs-Url: https://github.com/ChrisCreevey/catsequences/sessions/d6cfecf1-be96-4edd-82eb-e4900e7f01b7 Co-authored-by: ChrisCreevey <4115315+ChrisCreevey@users.noreply.github.com> --- catsequences.c | 40 ++++++++++++++++++++++++++--------- tests/expected/phylip_aln.phy | 2 +- tests/run_tests.sh | 2 +- 3 files changed, 32 insertions(+), 12 deletions(-) diff --git a/catsequences.c b/catsequences.c index c242591..6799c95 100755 --- a/catsequences.c +++ b/catsequences.c @@ -115,10 +115,11 @@ static int read_line(FILE *fp, char **buf, size_t *cap) break; } if (len + 1 >= *cap) { - *cap *= 2; - char *tmp = realloc(*buf, *cap); + size_t newcap = *cap * 2; + char *tmp = realloc(*buf, newcap); if (!tmp) die("out of memory"); *buf = tmp; + *cap = newcap; } (*buf)[len++] = (char)c; } @@ -126,12 +127,13 @@ static int read_line(FILE *fp, char **buf, size_t *cap) return (int)len; } -/* Return basename of path with the extension stripped. */ +/* Return basename of path with the last extension stripped. + * e.g. "genes/rbcL.fasta" -> "rbcL", "gene.v2.fasta" -> "gene.v2" */ static void gene_name(const char *path, char *out, size_t outsz) { const char *base = strrchr(path, '/'); base = base ? base + 1 : path; - const char *dot = strrchr(base, '.'); + const char *dot = strrchr(base, '.'); /* strip only the last extension */ size_t len = dot ? (size_t)(dot - base) : strlen(base); if (len >= outsz) len = outsz - 1; memcpy(out, base, len); @@ -375,13 +377,21 @@ static void second_pass(void) } } + /* Per-species write offset within the current file — avoids + * repeated strlen calls when appending multi-line sequences. */ + int *offsets = malloc((size_t)numspecies * sizeof(int)); + if (!offsets) die("out of memory"); + char *buf = NULL; size_t cap = 0; for (int j = 0; j < num_files; j++) { + memset(offsets, 0, (size_t)numspecies * sizeof(int)); + FILE *fp = fopen(filelist[j], "r"); if (!fp) { fprintf(stderr, "Error: Cannot reopen %s\n", filelist[j]); + free(offsets); exit(EXIT_FAILURE); } @@ -403,11 +413,12 @@ static void second_pass(void) if (found < 0) { fprintf(stderr, "Error: species '%s' not in master list\n", buf); + free(offsets); exit(EXIT_FAILURE); } /* Handle duplicates */ - if (sequences[j][found][0] != '\0') { + if (offsets[found] > 0) { switch (dup_mode) { case DUP_WARN: printf("Duplicate %s in %s: ignoring extra copy\n", @@ -418,18 +429,26 @@ static void second_pass(void) fprintf(stderr, "Error: duplicate species %s in %s\n", species_names[found] + 1, filelist[j]); + free(offsets); exit(EXIT_FAILURE); case DUP_KEEPLAST: + offsets[found] = 0; sequences[j][found][0] = '\0'; break; } } } else if (found >= 0 && !previous) { - /* Append sequence line, guarding against overflow */ - size_t have = strlen(sequences[j][found]); - size_t room = (size_t)seqlens[j] - have; - strncat(sequences[j][found], buf, room); + /* Append sequence line using tracked offset */ + int have = offsets[found]; + int room = seqlens[j] - have; + if (room > 0) { + int srclen = (int)strlen(buf); + int copy = srclen < room ? srclen : room; + memcpy(sequences[j][found] + have, buf, (size_t)copy); + offsets[found] = have + copy; + sequences[j][found][offsets[found]] = '\0'; + } } } fclose(fp); @@ -437,6 +456,7 @@ static void second_pass(void) if (verbose) printf("Read %s\n", filelist[j]); } + free(offsets); free(buf); } @@ -551,7 +571,7 @@ static void write_alignment(void) /* Format-specific headers */ if (aln_format == AFMT_PHYLIP) { - fprintf(fp, " %d %d\n", nsp, tot); + fprintf(fp, "%d %d\n", nsp, tot); } else if (aln_format == AFMT_NEXUS) { fprintf(fp, "#NEXUS\nbegin data;\n"); fprintf(fp, "\tdimensions ntax=%d nchar=%d;\n", nsp, tot); diff --git a/tests/expected/phylip_aln.phy b/tests/expected/phylip_aln.phy index 184605c..eb3aaa1 100644 --- a/tests/expected/phylip_aln.phy +++ b/tests/expected/phylip_aln.phy @@ -1,4 +1,4 @@ - 5 35 +5 35 species1 123456789abcdefghijklmno??????????? species2 1234-6789abcd???????????pqrstuvwxyz species3 1234-----a-b-e----jklmno??????????? diff --git a/tests/run_tests.sh b/tests/run_tests.sh index 2f84b94..1b9c39b 100755 --- a/tests/run_tests.sh +++ b/tests/run_tests.sh @@ -21,7 +21,7 @@ FAIL=0 # ── helper: compare two files, report result ───────────────────── check() { local label="$1" got="$2" want="$3" - if diff -q "$got" "$want" > /dev/null 2>&1; then + if diff -q "$want" "$got" > /dev/null 2>&1; then echo "PASS: $label" PASS=$(( PASS + 1 )) else