[Debian-med-packaging] Bug#653042: Bug#653042: segfault in mpileup when alignment is outside of the ref

[Charles Plessy]

Dear SAM tools developers,

here is a patch proposed by a Debian user on our bug tracker, that prevents
samtools mpileup to crash when a read is aligned outside of the reference.

Please consider including it to your next release, and do not hesitate to send
us comments if there are potential problems with this patch.

You can refer to the email quoted below or to the following URL for more
information.

  http://bugs.debian.org/653042

Have a nice day,

-- Charles Plessy, Debian Med team, Illkirch-Graffenstaden, France.


Le Thu, Dec 22, 2011 at 05:42:44PM -0800, Don Armstrong a écrit :
> Package: samtools
> Version: 0.1.18-1
> Severity: important
> Tag: patch
> 
> If a read is aligned outside of the reference (possibly due to
> reference mismatch or truncation), samtools mpileup will segfault. It
> should instead warn, and continue on. The attached patch fixes this.

> Description: Fix segfault if position is outside of the reference sequence length
> Origin: Don Armstrong <don at donarmstrong.com>
> 
> Index: samtools-0.1.18/bam_plcmd.c
> ===================================================================
> --- samtools-0.1.18.orig/bam_plcmd.c	2011-07-06 20:41:26.000000000 -0700
> +++ samtools-0.1.18/bam_plcmd.c	2011-12-22 14:44:05.000000000 -0800
> @@ -93,6 +93,7 @@
>  	bam_iter_t iter;
>  	bam_header_t *h;
>  	int ref_id;
> +        int ref_len;
>  	char *ref;
>  	const mplp_conf_t *conf;
>  } mplp_aux_t;
> @@ -134,6 +135,11 @@
>  				qual[i] = qual[i] > 31? qual[i] - 31 : 0;
>  		}
>  		has_ref = (ma->ref && ma->ref_id == b->core.tid)? 1 : 0;
> +		if (has_ref && (ma->ref_len <= b->core.pos)) { // exclude reads outside of the reference sequence
> +		  fprintf(stderr,"[%s] Skipping because %d is outside of %d [ref:%s]\n",__func__,b->core.pos,ma->ref_len,ma->ref_id);
> +		  skip = 1;
> +		  continue;
> +		}
>  		skip = 0;
>  		if (has_ref && (ma->conf->flag&MPLP_REALN)) bam_prob_realn_core(b, ma->ref, (ma->conf->flag & MPLP_EXT_BAQ)? 3 : 1);
>  		if (has_ref && ma->conf->capQ_thres > 10) {
> @@ -277,7 +283,11 @@
>  	if (tid0 >= 0 && conf->fai) { // region is set
>  		ref = faidx_fetch_seq(conf->fai, h->target_name[tid0], 0, 0x7fffffff, &ref_len);
>  		ref_tid = tid0;
> -		for (i = 0; i < n; ++i) data[i]->ref = ref, data[i]->ref_id = tid0;
> +		for (i = 0; i < n; ++i) {
> +		  data[i]->ref = ref;
> +                  data[i]->ref_id = tid0;
> +		  data[i]->ref_len = ref_len;
> +                }
>  	} else ref_tid = -1, ref = 0;
>  	iter = bam_mplp_init(n, mplp_func, (void**)data);
>  	max_depth = conf->max_depth;
> @@ -295,7 +305,11 @@
>  		if (tid != ref_tid) {
>  			free(ref); ref = 0;
>  			if (conf->fai) ref = faidx_fetch_seq(conf->fai, h->target_name[tid], 0, 0x7fffffff, &ref_len);
> -			for (i = 0; i < n; ++i) data[i]->ref = ref, data[i]->ref_id = tid;
> +			for (i = 0; i < n; ++i) {
> +			  data[i]->ref = ref;
> +			  data[i]->ref_id = tid;
> +			  data[i]->ref_len = ref_len;
> +			}
>  			ref_tid = tid;
>  		}
>  		if (conf->flag & MPLP_GLF) {