@@ -1366,6 +1366,9 @@ static void filters_set_alt_string(filter_t *flt, bcf1_t *line, token_t *tok)
13661366}
13671367static void filters_set_nmissing (filter_t * flt , bcf1_t * line , token_t * tok )
13681368{
1369+ if ( !tok -> nsamples ) error ("The function %s works with FORMAT fields\n" , tok -> tag );
1370+ assert (tok -> usmpl );
1371+
13691372 bcf_unpack (line , BCF_UN_FMT );
13701373 if ( !line -> n_sample )
13711374 {
@@ -1384,11 +1387,13 @@ static void filters_set_nmissing(filter_t *flt, bcf1_t *line, token_t *tok)
13841387 return ;
13851388 }
13861389
1387- int j ,nmissing = 0 ;
1390+ int j ,nsmpl = 0 , nmissing = 0 ;
13881391 #define BRANCH (type_t , convert , is_vector_end ) { \
13891392 for (i=0; i<line->n_sample; i++) \
13901393 { \
1394+ if ( !tok->usmpl[i] ) continue; \
13911395 uint8_t *ptr = fmt->p + i*fmt->size; \
1396+ nsmpl++; \
13921397 for (j=0; j<fmt->n; j++) \
13931398 { \
13941399 type_t val = convert(&ptr[j * sizeof(type_t)]); \
@@ -1405,7 +1410,7 @@ static void filters_set_nmissing(filter_t *flt, bcf1_t *line, token_t *tok)
14051410 }
14061411 #undef BRANCH
14071412 tok -> nvalues = 1 ;
1408- tok -> values [0 ] = tok -> tag [0 ]== 'N' ? nmissing : (double )nmissing / line -> n_sample ;
1413+ tok -> values [0 ] = tok -> tag [0 ]== 'N' ? nmissing : (nsmpl ? ( double )nmissing / nsmpl : 0 ) ;
14091414}
14101415static int func_npass (filter_t * flt , bcf1_t * line , token_t * rtok , token_t * * stack , int nstack )
14111416{
@@ -1414,16 +1419,17 @@ static int func_npass(filter_t *flt, bcf1_t *line, token_t *rtok, token_t **stac
14141419 if ( !tok -> nsamples ) error ("The function %s works with FORMAT fields\n" , rtok -> tag );
14151420 assert (tok -> usmpl );
14161421
1417- int i , npass = 0 ;
1422+ int i , nsmpl = 0 , npass = 0 ;
14181423 for (i = 0 ; i < tok -> nsamples ; i ++ )
14191424 {
14201425 if ( !tok -> usmpl [i ] ) continue ;
1426+ nsmpl ++ ;
14211427 if ( tok -> pass_samples [i ] ) npass ++ ;
14221428 }
14231429 hts_expand (double ,1 ,rtok -> mvalues ,rtok -> values );
14241430 rtok -> nsamples = 0 ;
14251431 rtok -> nvalues = 1 ;
1426- rtok -> values [0 ] = rtok -> tag [0 ]== 'N' ? npass : (line -> n_sample ? 1.0 * npass /line -> n_sample : 0 );
1432+ rtok -> values [0 ] = rtok -> tag [0 ]== 'N' ? npass : (nsmpl ? 1.0 * npass /nsmpl : 0 );
14271433
14281434 return 1 ;
14291435}
@@ -3144,6 +3150,14 @@ static int filters_init1_ext(filter_t *filter, char *str, int len, token_t *tok)
31443150 filter -> ext [filter -> n_ext - 1 ] = tok -> ht_type ;
31453151 return 0 ;
31463152}
3153+ static void init_usmpl (filter_t * flt , token_t * tok )
3154+ {
3155+ if ( tok -> nsamples ) return ;
3156+ int i ;
3157+ tok -> nsamples = bcf_hdr_nsamples (flt -> hdr );
3158+ tok -> usmpl = (uint8_t * ) malloc (tok -> nsamples );
3159+ for (i = 0 ; i < tok -> nsamples ; i ++ ) tok -> usmpl [i ] = 1 ;
3160+ }
31473161static int filters_init1 (filter_t * filter , char * str , int len , token_t * tok )
31483162{
31493163 tok -> vl_len = BCF_VL_FIXED ;
@@ -3295,6 +3309,7 @@ static int filters_init1(filter_t *filter, char *str, int len, token_t *tok)
32953309 tok -> setter = & filters_set_nmissing ;
32963310 tok -> tag = strdup ("N_MISSING" );
32973311 tok -> ht_type = BCF_HT_INT ;
3312+ init_usmpl (filter ,tok );
32983313 return 0 ;
32993314 }
33003315 else if ( !strncasecmp (str ,"F_MISSING" ,len ) )
@@ -3303,6 +3318,7 @@ static int filters_init1(filter_t *filter, char *str, int len, token_t *tok)
33033318 tok -> setter = & filters_set_nmissing ;
33043319 tok -> tag = strdup ("F_MISSING" );
33053320 tok -> ht_type = BCF_HT_REAL ;
3321+ init_usmpl (filter ,tok );
33063322 return 0 ;
33073323 }
33083324 }
@@ -3335,13 +3351,7 @@ static int filters_init1(filter_t *filter, char *str, int len, token_t *tok)
33353351 if ( tok -> idx == -3 && bcf_hdr_id2length (filter -> hdr ,BCF_HL_FMT ,tok -> hdr_id )!= BCF_VL_R )
33363352 error ("Error: GT subscripts can be used only with Number=R tags\n" );
33373353 }
3338- else if ( is_fmt && !tok -> nsamples )
3339- {
3340- int i ;
3341- tok -> nsamples = bcf_hdr_nsamples (filter -> hdr );
3342- tok -> usmpl = (uint8_t * ) malloc (tok -> nsamples );
3343- for (i = 0 ; i < tok -> nsamples ; i ++ ) tok -> usmpl [i ] = 1 ;
3344- }
3354+ else if ( is_fmt ) init_usmpl (filter ,tok );
33453355
33463356 tok -> hl_type = is_fmt ? BCF_HL_FMT : BCF_HL_INFO ;
33473357 if ( tok -> hdr_id >= 0 ) tok -> vl_len = bcf_hdr_id2length (filter -> hdr ,tok -> hl_type ,tok -> hdr_id );
0 commit comments