backup commit

[samtools.git] / bam_plcmd.c

diff --git a/bam_plcmd.c b/bam_plcmd.c
index 07f0a4fe3afb9bc13da71aabced9881ba49085ce..e45c8c8c81649630449f78d653cbd037e2046408 100644 (file)
--- a/bam_plcmd.c
+++ b/bam_plcmd.c
@@ -72,6 +72,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,8 +208,17 @@ 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);
+        if ( !h_tmp ) {
+            fprintf(stderr,"[%s] fail to read the header of %s\n", __func__, fn[i]);
+            exit(1);
+        }
                data[i]->h = i? h : h_tmp; // for i==0, "h" has not been set yet
                bam_smpl_add(sm, fn[i], (conf->flag&MPLP_IGNORE_RG)? 0 : h_tmp->text);
                rghash = bcf_call_add_rg(rghash, h_tmp->text, conf->pl_list);
@@ -217,11 +227,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 +281,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);
@@ -305,7 +316,7 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn)
                        ref16 = bam_nt16_table[_ref0];
                        for (i = 0; i < gplp.n; ++i)
                                bcf_call_glfgen(gplp.n_plp[i], gplp.plp[i], ref16, bca, bcr + i);
-                       bcf_call_combine(gplp.n, bcr, ref16, &bc);
+                       bcf_call_combine(gplp.n, bcr, bca, ref16, &bc);
                        bcf_call2bcf(tid, pos, &bc, b, bcr, conf->fmt_flag, 0, 0);
                        bcf_write(bp, bh, b);
                        bcf_destroy(b);
@@ -313,7 +324,7 @@ static int mpileup(mplp_conf_t *conf, int n, char **fn)
                        if (!(conf->flag&MPLP_NO_INDEL) && total_depth < max_indel_depth && bcf_call_gap_prep(gplp.n, gplp.n_plp, gplp.plp, pos, bca, ref, rghash) >= 0) {
                                for (i = 0; i < gplp.n; ++i)
                                        bcf_call_glfgen(gplp.n_plp[i], gplp.plp[i], -1, bca, bcr + i);
-                               if (bcf_call_combine(gplp.n, bcr, -1, &bc) >= 0) {
+                               if (bcf_call_combine(gplp.n, bcr, bca, -1, &bc) >= 0) {
                                        b = calloc(1, sizeof(bcf1_t));
                                        bcf_call2bcf(tid, pos, &bc, b, bcr, conf->fmt_flag, bca, ref);
                                        bcf_write(bp, bh, b);
@@ -449,7 +460,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 +470,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 +544,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");