Merge remote branch 'remotes/master/master' into fb-annots

[samtools.git] / bam_plcmd.c

diff --git a/bam_plcmd.c b/bam_plcmd.c
index e45c8c8c81649630449f78d653cbd037e2046408..9c4f273dcbfe2c3576e9c6ea0e2d7661839bbb85 100644 (file)
--- a/bam_plcmd.c
+++ b/bam_plcmd.c
@@ -4,6 +4,8 @@
 #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"
@@ -80,6 +82,7 @@ int bed_overlap(const void *_h, const char *chr, int beg, int end);
 
 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;
@@ -117,6 +120,8 @@ static int mplp_func(void *data, bam1_t *b)
                        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;
@@ -397,11 +402,15 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn)
 }
 
 #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 )
@@ -410,28 +419,33 @@ static int read_file_list(const char *file_list,int *n,char **argv[])
         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 )
@@ -460,8 +474,16 @@ int bam_mpileup(int argc, char *argv[])
        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;
@@ -517,7 +539,7 @@ int bam_mpileup(int argc, char *argv[])
                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");
@@ -529,6 +551,8 @@ int bam_mpileup(int argc, char *argv[])
                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");