X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=bam_plcmd.c;h=a51e7731b0616a8dcf0c4f709f20544f4997ca47;hb=4b0058f506960786fc4c988c512cd462e02f24ed;hp=07f0a4fe3afb9bc13da71aabced9881ba49085ce;hpb=183e4cb41c6136b9bff3dc83fcfa9f023e1d7629;p=samtools.git

diff --git a/bam_plcmd.c b/bam_plcmd.c
index 07f0a4f..a51e773 100644
--- a/bam_plcmd.c
+++ b/bam_plcmd.c
@@ -4,6 +4,7 @@
 #include <ctype.h>
 #include <string.h>
 #include <errno.h>
+#include <sys/stat.h>
 #include "sam.h"
 #include "faidx.h"
 #include "kstring.h"
@@ -72,6 +73,7 @@ static inline void pileup_seq(const bam_pileup1_t *p, int pos, int ref_len, cons
 #define MPLP_IGNORE_RG 0x2000
 #define MPLP_PRINT_POS 0x4000
 #define MPLP_PRINT_MAPQ 0x8000
+#define MPLP_PER_SAMPLE 0x10000
 
 void *bed_read(const char *fn);
 void bed_destroy(void *_h);
@@ -207,6 +209,11 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn)
 		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
@@ -217,11 +224,11 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn)
 			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;
@@ -271,6 +278,7 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn)
 		bca->openQ = conf->openQ, bca->extQ = conf->extQ, bca->tandemQ = conf->tandemQ;
 		bca->min_frac = conf->min_frac;
 		bca->min_support = conf->min_support;
+        bca->per_sample_flt = conf->flag & MPLP_PER_SAMPLE;
 	}
 	if (tid0 >= 0 && conf->fai) { // region is set
 		ref = faidx_fetch_seq(conf->fai, h->target_name[tid0], 0, 0x7fffffff, &ref_len);
@@ -389,8 +397,12 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn)
 static 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 )
@@ -399,28 +411,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 )
@@ -449,7 +466,7 @@ 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:o:e:h:Im:F:EG:6OsV")) >= 0) {
+	while ((c = getopt(argc, argv, "Agf:r:l:M:q:Q:uaRC:BDSd:L:b:P:po:e:h:Im:F:EG:6OsV")) >= 0) {
 		switch (c) {
 		case 'f':
 			mplp.fai = fai_load(optarg);
@@ -459,6 +476,7 @@ int bam_mpileup(int argc, char *argv[])
 		case 'r': mplp.reg = strdup(optarg); break;
 		case 'l': mplp.bed = bed_read(optarg); break;
 		case 'P': mplp.pl_list = strdup(optarg); break;
+		case 'p': mplp.flag |= MPLP_PER_SAMPLE; break;
 		case 'g': mplp.flag |= MPLP_GLF; break;
 		case 'u': mplp.flag |= MPLP_NO_COMP | MPLP_GLF; break;
 		case 'a': mplp.flag |= MPLP_NO_ORPHAN | MPLP_REALN; break;
@@ -532,6 +550,7 @@ int bam_mpileup(int argc, char *argv[])
 		fprintf(stderr, "       -L INT       max per-sample depth for INDEL calling [%d]\n", mplp.max_indel_depth);
 		fprintf(stderr, "       -m INT       minimum gapped reads for indel candidates [%d]\n", mplp.min_support);
 		fprintf(stderr, "       -o INT       Phred-scaled gap open sequencing error probability [%d]\n", mplp.openQ);
+		fprintf(stderr, "       -p           apply -m and -F per-sample to increase sensitivity\n");
 		fprintf(stderr, "       -P STR       comma separated list of platforms for indels [all]\n");
 		fprintf(stderr, "\n");
 		fprintf(stderr, "Notes: Assuming diploid individuals.\n\n");