X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=bam_plcmd.c;h=391767ab134833f5dcc5ea914866fd7f10b97718;hb=9bd8cd7dd42efa091a9361bb06a29d1bead405ce;hp=46f078a99e715a89ccc4f73119192fe8d1ef280c;hpb=4c4b72301b736e08e37b5e0dd1d435477d969898;p=samtools.git

diff --git a/bam_plcmd.c b/bam_plcmd.c
index 46f078a..391767a 100644
--- a/bam_plcmd.c
+++ b/bam_plcmd.c
@@ -17,6 +17,7 @@ KHASH_MAP_INIT_INT64(64, indel_list_t)
 #define BAM_PLF_INDEL_ONLY 0x04
 #define BAM_PLF_GLF        0x08
 #define BAM_PLF_VAR_ONLY   0x10
+#define BAM_PLF_2ND        0x20
 
 typedef struct {
 	bam_header_t *h;
@@ -189,8 +190,10 @@ static int pileup_func(uint32_t tid, uint32_t pos, int n, const bam_pileup1_t *p
 	// when only variant sites are asked for, test if the site is a variant
 	if ((d->format & BAM_PLF_CNS) && (d->format & BAM_PLF_VAR_ONLY)) {
 		if (!(bam_nt16_table[rb] != 15 && cns>>28 != bam_nt16_table[rb])) { // not a SNP
-			if (!(r && (r->gt == 2 || strcmp(r->s[r->gt], "*")))) // not an indel
+			if (!(r && (r->gt == 2 || strcmp(r->s[r->gt], "*")))) { // not an indel
+				if (r) bam_maqindel_ret_destroy(r);
 				return 0;
+			}
 		}
 	}
 	// print the first 3 columns
@@ -246,6 +249,21 @@ static int pileup_func(uint32_t tid, uint32_t pos, int n, const bam_pileup1_t *p
 		if (c > 126) c = 126;
 		putchar(c);
 	}
+	if (d->format & BAM_PLF_2ND) { // print 2nd calls and qualities
+		const unsigned char *q;
+		putchar('\t');
+		for (i = 0; i < n; ++i) {
+			const bam_pileup1_t *p = pu + i;
+			q = bam_aux_get(p->b, "E2");
+			putchar(q? q[p->qpos + 1] : 'N');
+		}
+		putchar('\t');
+		for (i = 0; i < n; ++i) {
+			const bam_pileup1_t *p = pu + i;
+			q = bam_aux_get(p->b, "U2");
+			putchar(q? q[p->qpos + 1] : '!');
+		}
+	}
 	// print mapping quality if -s is flagged on the command line
 	if (d->format & BAM_PLF_SIMPLE) {
 		putchar('\t');
@@ -274,13 +292,13 @@ static int pileup_func(uint32_t tid, uint32_t pos, int n, const bam_pileup1_t *p
 
 int bam_pileup(int argc, char *argv[])
 {
-	int c;
+	int c, is_SAM = 0;
 	char *fn_list = 0, *fn_fa = 0, *fn_pos = 0;
 	pu_data_t *d = (pu_data_t*)calloc(1, sizeof(pu_data_t));
 	d->tid = -1; d->mask = BAM_DEF_MASK;
 	d->c = bam_maqcns_init();
 	d->ido = bam_maqindel_opt_init();
-	while ((c = getopt(argc, argv, "st:f:cT:N:r:l:im:gI:G:vM:")) >= 0) {
+	while ((c = getopt(argc, argv, "st:f:cT:N:r:l:im:gI:G:vM:S2")) >= 0) {
 		switch (c) {
 		case 's': d->format |= BAM_PLF_SIMPLE; break;
 		case 't': fn_list = strdup(optarg); break;
@@ -295,19 +313,24 @@ int bam_pileup(int argc, char *argv[])
 		case 'v': d->format |= BAM_PLF_VAR_ONLY; break;
 		case 'm': d->mask = strtol(optarg, 0, 0); break;
 		case 'g': d->format |= BAM_PLF_GLF; break;
+		case '2': d->format |= BAM_PLF_2ND; break;
 		case 'I': d->ido->q_indel = atoi(optarg); break;
 		case 'G': d->ido->r_indel = atof(optarg); break;
+		case 'S': is_SAM = 1; break;
 		default: fprintf(stderr, "Unrecognizd option '-%c'.\n", c); return 1;
 		}
 	}
+	if (fn_list) is_SAM = 1;
 	if (optind == argc) {
 		fprintf(stderr, "\n");
 		fprintf(stderr, "Usage:  samtools pileup [options] <in.bam>|<in.sam>\n\n");
 		fprintf(stderr, "Option: -s        simple (yet incomplete) pileup format\n");
+		fprintf(stderr, "        -S        the input is in SAM\n");
+		fprintf(stderr, "        -2        output the 2nd best call and quality\n");
 		fprintf(stderr, "        -i        only show lines/consensus with indels\n");
 		fprintf(stderr, "        -m INT    filtering reads with bits in INT [%d]\n", d->mask);
 		fprintf(stderr, "        -M INT    cap mapping quality at INT [%d]\n", d->c->cap_mapQ);
-		fprintf(stderr, "        -t FILE   list of reference sequences (assume the input is in SAM)\n");
+		fprintf(stderr, "        -t FILE   list of reference sequences (force -S)\n");
 		fprintf(stderr, "        -l FILE   list of sites at which pileup is output\n");
 		fprintf(stderr, "        -f FILE   reference sequence in the FASTA format\n\n");
 		fprintf(stderr, "        -c        output the maq consensus sequence\n");
@@ -333,9 +356,11 @@ int bam_pileup(int argc, char *argv[])
 	}
 	if (d->fai == 0 && (d->format & (BAM_PLF_CNS|BAM_PLF_INDEL_ONLY)))
 		fprintf(stderr, "[bam_pileup] indels will not be called when -f is absent.\n");
+	if (fn_fa && is_SAM && fn_list == 0) fn_list = samfaipath(fn_fa);
+
 	{
 		samfile_t *fp;
-		fp = fn_list? samopen(argv[optind], "r", fn_list) : samopen(argv[optind], "rb", 0);
+		fp = is_SAM? samopen(argv[optind], "r", fn_list) : samopen(argv[optind], "rb", 0);
 		if (fp == 0 || fp->header == 0) {
 			fprintf(stderr, "[bam_pileup] fail to read the header: non-exisiting file or wrong format.\n");
 			return 1;