#include <ctype.h>
#include <string.h>
#include <errno.h>
+#include <sys/stat.h>
+#include <getopt.h>
#include "sam.h"
#include "faidx.h"
#include "kstring.h"
typedef struct {
int max_mq, min_mq, flag, min_baseQ, capQ_thres, max_depth, max_indel_depth, fmt_flag;
+ int rflag_require, rflag_filter;
int openQ, extQ, tandemQ, min_support; // for indels
double min_frac; // for indels
char *reg, *pl_list;
skip = 1;
continue;
}
+ if (ma->conf->rflag_require && !(ma->conf->rflag_require&b->core.flag)) { skip = 1; continue; }
+ if (ma->conf->rflag_filter && ma->conf->rflag_filter&b->core.flag) { skip = 1; continue; }
if (ma->conf->bed) { // test overlap
skip = !bed_overlap(ma->conf->bed, ma->h->target_name[b->core.tid], b->core.pos, bam_calend(&b->core, bam1_cigar(b)));
if (skip) continue;
bam_header_t *h_tmp;
data[i] = calloc(1, sizeof(mplp_aux_t));
data[i]->fp = strcmp(fn[i], "-") == 0? bam_dopen(fileno(stdin), "r") : bam_open(fn[i], "r");
+ if ( !data[i]->fp )
+ {
+ fprintf(stderr, "[%s] failed to open %s: %s\n", __func__, fn[i], strerror(errno));
+ exit(1);
+ }
data[i]->conf = conf;
h_tmp = bam_header_read(data[i]->fp);
data[i]->h = i? h : h_tmp; // for i==0, "h" has not been set yet
bam_index_t *idx;
idx = bam_index_load(fn[i]);
if (idx == 0) {
- fprintf(stderr, "[%s] fail to load index for %d-th input.\n", __func__, i+1);
+ fprintf(stderr, "[%s] fail to load index for %s\n", __func__, fn[i]);
exit(1);
}
if (bam_parse_region(h_tmp, conf->reg, &tid, &beg, &end) < 0) {
- fprintf(stderr, "[%s] malformatted region or wrong seqname for %d-th input.\n", __func__, i+1);
+ fprintf(stderr, "[%s] malformatted region or wrong seqname for %s\n", __func__, fn[i]);
exit(1);
}
if (i == 0) tid0 = tid, beg0 = beg, end0 = end;
}
#define MAX_PATH_LEN 1024
-static int read_file_list(const char *file_list,int *n,char **argv[])
+int read_file_list(const char *file_list,int *n,char **argv[])
{
char buf[MAX_PATH_LEN];
- int len, nfiles;
- char **files;
+ int len, nfiles = 0;
+ char **files = NULL;
+ struct stat sb;
+
+ *n = 0;
+ *argv = NULL;
FILE *fh = fopen(file_list,"r");
if ( !fh )
return 1;
}
- // Speed is not an issue here, determine the number of files by reading the file twice
- nfiles = 0;
- while ( fgets(buf,MAX_PATH_LEN,fh) ) nfiles++;
-
- if ( fseek(fh, 0L, SEEK_SET) )
- {
- fprintf(stderr,"%s: %s\n", file_list,strerror(errno));
- return 1;
- }
-
files = calloc(nfiles,sizeof(char*));
nfiles = 0;
while ( fgets(buf,MAX_PATH_LEN,fh) )
{
+ // allow empty lines and trailing spaces
len = strlen(buf);
while ( len>0 && isspace(buf[len-1]) ) len--;
if ( !len ) continue;
- files[nfiles] = malloc(sizeof(char)*(len+1));
- strncpy(files[nfiles],buf,len);
- files[nfiles][len] = 0;
+ // check sanity of the file list
+ buf[len] = 0;
+ if (stat(buf, &sb) != 0)
+ {
+ // no such file, check if it is safe to print its name
+ int i, safe_to_print = 1;
+ for (i=0; i<len; i++)
+ if (!isprint(buf[i])) { safe_to_print = 0; break; }
+ if ( safe_to_print )
+ fprintf(stderr,"The file list \"%s\" appears broken, could not locate: %s\n", file_list,buf);
+ else
+ fprintf(stderr,"Does the file \"%s\" really contain a list of files and do all exist?\n", file_list);
+ return 1;
+ }
+
nfiles++;
+ files = realloc(files,nfiles*sizeof(char*));
+ files[nfiles-1] = strdup(buf);
}
fclose(fh);
if ( !nfiles )
mplp.openQ = 40; mplp.extQ = 20; mplp.tandemQ = 100;
mplp.min_frac = 0.002; mplp.min_support = 1;
mplp.flag = MPLP_NO_ORPHAN | MPLP_REALN;
- while ((c = getopt(argc, argv, "Agf:r:l:M:q:Q:uaRC:BDSd:L:b:P:po:e:h:Im:F:EG:6OsV")) >= 0) {
+ static struct option lopts[] =
+ {
+ {"rf",1,0,1}, // require flag
+ {"ff",1,0,2}, // filter flag
+ {0,0,0,0}
+ };
+ while ((c = getopt_long(argc, argv, "Agf:r:l:M:q:Q:uaRC:BDSd:L:b:P:po:e:h:Im:F:EG:6OsV1:2:",lopts,NULL)) >= 0) {
switch (c) {
+ case 1 : mplp.rflag_require = strtol(optarg,0,0); break;
+ case 2 : mplp.rflag_filter = strtol(optarg,0,0); break;
case 'f':
mplp.fai = fai_load(optarg);
if (mplp.fai == 0) return 1;
fprintf(stderr, " -6 assume the quality is in the Illumina-1.3+ encoding\n");
fprintf(stderr, " -A count anomalous read pairs\n");
fprintf(stderr, " -B disable BAQ computation\n");
- fprintf(stderr, " -b FILE list of input BAM files [null]\n");
+ fprintf(stderr, " -b FILE list of input BAM filenames, one per line [null]\n");
fprintf(stderr, " -C INT parameter for adjusting mapQ; 0 to disable [0]\n");
fprintf(stderr, " -d INT max per-BAM depth to avoid excessive memory usage [%d]\n", mplp.max_depth);
fprintf(stderr, " -E recalculate extended BAQ on the fly thus ignoring existing BQs\n");
fprintf(stderr, " -R ignore RG tags\n");
fprintf(stderr, " -q INT skip alignments with mapQ smaller than INT [%d]\n", mplp.min_mq);
fprintf(stderr, " -Q INT skip bases with baseQ/BAQ smaller than INT [%d]\n", mplp.min_baseQ);
+ fprintf(stderr, " --rf INT required flags: skip reads with mask bits unset []\n");
+ fprintf(stderr, " --ff INT filter flags: skip reads with mask bits set []\n");
fprintf(stderr, "\nOutput options:\n\n");
fprintf(stderr, " -D output per-sample DP in BCF (require -g/-u)\n");
fprintf(stderr, " -g generate BCF output (genotype likelihoods)\n");