Loading NEWS +12 −0 Original line number Diff line number Diff line Loading @@ -5,6 +5,18 @@ http://pysam.readthedocs.io/en/latest/release.html Release notes ============= Release 0.14.1 ============== This is mostly a bugfix release, though bcftools has now also been upgraded to 1.7.0. * [#621] Add a warning to count_coverage when an alignment has an empty QUAL field * [#635] Speed-up of AlignedSegment.find_intro() * treat border case of all bases in pileup column below quality score * [#634] Fix access to pileup reference_sequence Release 0.14.0 ============== Loading bcftools/bcftools.h +17 −0 Original line number Diff line number Diff line Loading @@ -63,6 +63,23 @@ static inline char gt2iupac(char a, char b) return iupac[(int)a][(int)b]; } static inline int iupac_consistent(char iupac, char nt) { static const char iupac_mask[90] = { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,14,2, 13,0,0,4,11,0,0,12,0,3,15,0,0,0,5,6,8,0,7,9,0,10 }; if ( iupac > 89 ) return 0; if ( nt > 90 ) nt -= 32; // lowercase if ( nt=='A' ) nt = 1; else if ( nt=='C' ) nt = 2; else if ( nt=='G' ) nt = 4; else if ( nt=='T' ) nt = 8; return iupac_mask[(int)iupac] & nt ? 1 : 0; } static inline char nt_to_upper(char nt) { if ( nt < 97 ) return nt; Loading bcftools/consensus.c +9 −4 Original line number Diff line number Diff line Loading @@ -77,6 +77,7 @@ typedef struct rbuf_t vcf_rbuf; bcf1_t **vcf_buf; int nvcf_buf, rid; char *chr; regidx_t *mask; regitr_t *itr; Loading Loading @@ -121,6 +122,8 @@ static void destroy_chain(args_t *args) free(chain->block_lengths); free(chain); chain = NULL; free(args->chr); args->chr = NULL; } static void print_chain(args_t *args) Loading Loading @@ -162,7 +165,7 @@ static void print_chain(args_t *args) score += chain->block_lengths[n]; } score += last_block_size; fprintf(args->fp_chain, "chain %d %s %d + %d %d %s %d + %d %d %d\n", score, bcf_hdr_id2name(args->hdr,args->rid), ref_end_pos, chain->ori_pos, ref_end_pos, bcf_hdr_id2name(args->hdr,args->rid), alt_end_pos, chain->ori_pos, alt_end_pos, ++args->chain_id); fprintf(args->fp_chain, "chain %d %s %d + %d %d %s %d + %d %d %d\n", score, args->chr, ref_end_pos, chain->ori_pos, ref_end_pos, args->chr, alt_end_pos, chain->ori_pos, alt_end_pos, ++args->chain_id); for (n=0; n<chain->num; n++) { fprintf(args->fp_chain, "%d %d %d\n", chain->block_lengths[n], chain->ref_gaps[n], chain->alt_gaps[n]); } Loading Loading @@ -248,6 +251,7 @@ static void destroy_data(args_t *args) if ( args->vcf_buf[i] ) bcf_destroy1(args->vcf_buf[i]); free(args->vcf_buf); free(args->fa_buf.s); free(args->chr); if ( args->mask ) regidx_destroy(args->mask); if ( args->itr ) regitr_destroy(args->itr); if ( args->chain_fname ) Loading Loading @@ -276,6 +280,7 @@ static void init_region(args_t *args, char *line) else to--; } } args->chr = strdup(line); args->rid = bcf_hdr_name2id(args->hdr,line); if ( args->rid<0 ) fprintf(stderr,"Warning: Sequence \"%s\" not in %s\n", line,args->fname); args->fa_buf.l = 0; Loading Loading @@ -380,7 +385,7 @@ static void apply_variant(args_t *args, bcf1_t *rec) if ( regidx_overlap(args->mask, chr,start,end,NULL) ) return; } int i, ialt = 1; int i, ialt = 1; // the alternate allele if ( args->isample >= 0 ) { bcf_unpack(rec, BCF_UN_FMT); Loading Loading @@ -417,6 +422,7 @@ static void apply_variant(args_t *args, bcf1_t *rec) { char ial = rec->d.allele[ialt][0]; char jal = rec->d.allele[jalt][0]; if ( !ialt ) ialt = jalt; // only ialt is used, make sure 0/1 is not ignored rec->d.allele[ialt][0] = gt2iupac(ial,jal); } } Loading Loading @@ -565,11 +571,10 @@ static void apply_variant(args_t *args, bcf1_t *rec) static void mask_region(args_t *args, char *seq, int len) { char *chr = (char*)bcf_hdr_id2name(args->hdr,args->rid); int start = args->fa_src_pos - len; int end = args->fa_src_pos; if ( !regidx_overlap(args->mask, chr,start,end, args->itr) ) return; if ( !regidx_overlap(args->mask, args->chr,start,end, args->itr) ) return; int idx_start, idx_end, i; while ( regitr_overlap(args->itr) ) Loading bcftools/consensus.c.pysam.c +9 −4 Original line number Diff line number Diff line Loading @@ -79,6 +79,7 @@ typedef struct rbuf_t vcf_rbuf; bcf1_t **vcf_buf; int nvcf_buf, rid; char *chr; regidx_t *mask; regitr_t *itr; Loading Loading @@ -123,6 +124,8 @@ static void destroy_chain(args_t *args) free(chain->block_lengths); free(chain); chain = NULL; free(args->chr); args->chr = NULL; } static void print_chain(args_t *args) Loading Loading @@ -164,7 +167,7 @@ static void print_chain(args_t *args) score += chain->block_lengths[n]; } score += last_block_size; fprintf(args->fp_chain, "chain %d %s %d + %d %d %s %d + %d %d %d\n", score, bcf_hdr_id2name(args->hdr,args->rid), ref_end_pos, chain->ori_pos, ref_end_pos, bcf_hdr_id2name(args->hdr,args->rid), alt_end_pos, chain->ori_pos, alt_end_pos, ++args->chain_id); fprintf(args->fp_chain, "chain %d %s %d + %d %d %s %d + %d %d %d\n", score, args->chr, ref_end_pos, chain->ori_pos, ref_end_pos, args->chr, alt_end_pos, chain->ori_pos, alt_end_pos, ++args->chain_id); for (n=0; n<chain->num; n++) { fprintf(args->fp_chain, "%d %d %d\n", chain->block_lengths[n], chain->ref_gaps[n], chain->alt_gaps[n]); } Loading Loading @@ -250,6 +253,7 @@ static void destroy_data(args_t *args) if ( args->vcf_buf[i] ) bcf_destroy1(args->vcf_buf[i]); free(args->vcf_buf); free(args->fa_buf.s); free(args->chr); if ( args->mask ) regidx_destroy(args->mask); if ( args->itr ) regitr_destroy(args->itr); if ( args->chain_fname ) Loading Loading @@ -278,6 +282,7 @@ static void init_region(args_t *args, char *line) else to--; } } args->chr = strdup(line); args->rid = bcf_hdr_name2id(args->hdr,line); if ( args->rid<0 ) fprintf(bcftools_stderr,"Warning: Sequence \"%s\" not in %s\n", line,args->fname); args->fa_buf.l = 0; Loading Loading @@ -382,7 +387,7 @@ static void apply_variant(args_t *args, bcf1_t *rec) if ( regidx_overlap(args->mask, chr,start,end,NULL) ) return; } int i, ialt = 1; int i, ialt = 1; // the alternate allele if ( args->isample >= 0 ) { bcf_unpack(rec, BCF_UN_FMT); Loading Loading @@ -419,6 +424,7 @@ static void apply_variant(args_t *args, bcf1_t *rec) { char ial = rec->d.allele[ialt][0]; char jal = rec->d.allele[jalt][0]; if ( !ialt ) ialt = jalt; // only ialt is used, make sure 0/1 is not ignored rec->d.allele[ialt][0] = gt2iupac(ial,jal); } } Loading Loading @@ -567,11 +573,10 @@ static void apply_variant(args_t *args, bcf1_t *rec) static void mask_region(args_t *args, char *seq, int len) { char *chr = (char*)bcf_hdr_id2name(args->hdr,args->rid); int start = args->fa_src_pos - len; int end = args->fa_src_pos; if ( !regidx_overlap(args->mask, chr,start,end, args->itr) ) return; if ( !regidx_overlap(args->mask, args->chr,start,end, args->itr) ) return; int idx_start, idx_end, i; while ( regitr_overlap(args->itr) ) Loading bcftools/convert.c +29 −1 Original line number Diff line number Diff line /* convert.c -- functions for converting between VCF/BCF and related formats. Copyright (C) 2013-2017 Genome Research Ltd. Copyright (C) 2013-2018 Genome Research Ltd. Author: Petr Danecek <pd3@sanger.ac.uk> Loading Loading @@ -92,6 +92,7 @@ struct _convert_t int ndat; char *undef_info_tag; int allow_undef_tags; uint8_t **subset_samples; }; typedef struct Loading Loading @@ -174,6 +175,24 @@ static inline int32_t bcf_array_ivalue(void *bcf_array, int type, int idx) } return ((int32_t*)bcf_array)[idx]; } static inline void _copy_field(char *src, uint32_t len, int idx, kstring_t *str) { int n = 0, ibeg = 0; while ( src[ibeg] && ibeg<len && n < idx ) { if ( src[ibeg]==',' ) n++; ibeg++; } if ( ibeg==len ) { kputc('.', str); return; } int iend = ibeg; while ( src[iend] && src[iend]!=',' && iend<len ) iend++; if ( iend>ibeg ) kputsn(src+ibeg, iend-ibeg, str); else kputc('.', str); } static void process_info(convert_t *convert, bcf1_t *line, fmt_t *fmt, int isample, kstring_t *str) { if ( fmt->id<0 ) Loading Loading @@ -232,6 +251,7 @@ static void process_info(convert_t *convert, bcf1_t *line, fmt_t *fmt, int isamp case BCF_BT_INT16: BRANCH(int16_t, val==bcf_int16_missing, val==bcf_int16_vector_end, kputw(val, str)); break; case BCF_BT_INT32: BRANCH(int32_t, val==bcf_int32_missing, val==bcf_int32_vector_end, kputw(val, str)); break; case BCF_BT_FLOAT: BRANCH(float, bcf_float_is_missing(val), bcf_float_is_vector_end(val), kputd(val, str)); break; case BCF_BT_CHAR: _copy_field((char*)info->vptr, info->vptr_len, fmt->subscript, str); break; default: fprintf(stderr,"todo: type %d\n", info->type); exit(1); break; } #undef BRANCH Loading Loading @@ -288,6 +308,8 @@ static void process_format(convert_t *convert, bcf1_t *line, fmt_t *fmt, int isa else kputw(ival, str); } else if ( fmt->fmt->type == BCF_BT_CHAR ) _copy_field((char*)(fmt->fmt->p + isample*fmt->fmt->size), fmt->fmt->size, fmt->subscript, str); else error("TODO: %s:%d .. fmt->type=%d\n", __FILE__,__LINE__, fmt->fmt->type); } else Loading Loading @@ -1312,6 +1334,9 @@ int convert_line(convert_t *convert, bcf1_t *line, kstring_t *str) } for (js=0; js<convert->nsamples; js++) { // Skip samples when filtering was requested if ( *convert->subset_samples && !(*convert->subset_samples)[js] ) continue; // Here comes a hack designed for TBCSQ. When running on large files, // such as 1000GP, there are too many empty fields in the output and // it's very very slow. Therefore in case the handler does not add Loading Loading @@ -1362,6 +1387,9 @@ int convert_set_option(convert_t *convert, enum convert_option opt, ...) case allow_undef_tags: convert->allow_undef_tags = va_arg(args, int); break; case subset_samples: convert->subset_samples = va_arg(args, uint8_t**); break; default: ret = -1; } Loading Loading
NEWS +12 −0 Original line number Diff line number Diff line Loading @@ -5,6 +5,18 @@ http://pysam.readthedocs.io/en/latest/release.html Release notes ============= Release 0.14.1 ============== This is mostly a bugfix release, though bcftools has now also been upgraded to 1.7.0. * [#621] Add a warning to count_coverage when an alignment has an empty QUAL field * [#635] Speed-up of AlignedSegment.find_intro() * treat border case of all bases in pileup column below quality score * [#634] Fix access to pileup reference_sequence Release 0.14.0 ============== Loading
bcftools/bcftools.h +17 −0 Original line number Diff line number Diff line Loading @@ -63,6 +63,23 @@ static inline char gt2iupac(char a, char b) return iupac[(int)a][(int)b]; } static inline int iupac_consistent(char iupac, char nt) { static const char iupac_mask[90] = { 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,14,2, 13,0,0,4,11,0,0,12,0,3,15,0,0,0,5,6,8,0,7,9,0,10 }; if ( iupac > 89 ) return 0; if ( nt > 90 ) nt -= 32; // lowercase if ( nt=='A' ) nt = 1; else if ( nt=='C' ) nt = 2; else if ( nt=='G' ) nt = 4; else if ( nt=='T' ) nt = 8; return iupac_mask[(int)iupac] & nt ? 1 : 0; } static inline char nt_to_upper(char nt) { if ( nt < 97 ) return nt; Loading
bcftools/consensus.c +9 −4 Original line number Diff line number Diff line Loading @@ -77,6 +77,7 @@ typedef struct rbuf_t vcf_rbuf; bcf1_t **vcf_buf; int nvcf_buf, rid; char *chr; regidx_t *mask; regitr_t *itr; Loading Loading @@ -121,6 +122,8 @@ static void destroy_chain(args_t *args) free(chain->block_lengths); free(chain); chain = NULL; free(args->chr); args->chr = NULL; } static void print_chain(args_t *args) Loading Loading @@ -162,7 +165,7 @@ static void print_chain(args_t *args) score += chain->block_lengths[n]; } score += last_block_size; fprintf(args->fp_chain, "chain %d %s %d + %d %d %s %d + %d %d %d\n", score, bcf_hdr_id2name(args->hdr,args->rid), ref_end_pos, chain->ori_pos, ref_end_pos, bcf_hdr_id2name(args->hdr,args->rid), alt_end_pos, chain->ori_pos, alt_end_pos, ++args->chain_id); fprintf(args->fp_chain, "chain %d %s %d + %d %d %s %d + %d %d %d\n", score, args->chr, ref_end_pos, chain->ori_pos, ref_end_pos, args->chr, alt_end_pos, chain->ori_pos, alt_end_pos, ++args->chain_id); for (n=0; n<chain->num; n++) { fprintf(args->fp_chain, "%d %d %d\n", chain->block_lengths[n], chain->ref_gaps[n], chain->alt_gaps[n]); } Loading Loading @@ -248,6 +251,7 @@ static void destroy_data(args_t *args) if ( args->vcf_buf[i] ) bcf_destroy1(args->vcf_buf[i]); free(args->vcf_buf); free(args->fa_buf.s); free(args->chr); if ( args->mask ) regidx_destroy(args->mask); if ( args->itr ) regitr_destroy(args->itr); if ( args->chain_fname ) Loading Loading @@ -276,6 +280,7 @@ static void init_region(args_t *args, char *line) else to--; } } args->chr = strdup(line); args->rid = bcf_hdr_name2id(args->hdr,line); if ( args->rid<0 ) fprintf(stderr,"Warning: Sequence \"%s\" not in %s\n", line,args->fname); args->fa_buf.l = 0; Loading Loading @@ -380,7 +385,7 @@ static void apply_variant(args_t *args, bcf1_t *rec) if ( regidx_overlap(args->mask, chr,start,end,NULL) ) return; } int i, ialt = 1; int i, ialt = 1; // the alternate allele if ( args->isample >= 0 ) { bcf_unpack(rec, BCF_UN_FMT); Loading Loading @@ -417,6 +422,7 @@ static void apply_variant(args_t *args, bcf1_t *rec) { char ial = rec->d.allele[ialt][0]; char jal = rec->d.allele[jalt][0]; if ( !ialt ) ialt = jalt; // only ialt is used, make sure 0/1 is not ignored rec->d.allele[ialt][0] = gt2iupac(ial,jal); } } Loading Loading @@ -565,11 +571,10 @@ static void apply_variant(args_t *args, bcf1_t *rec) static void mask_region(args_t *args, char *seq, int len) { char *chr = (char*)bcf_hdr_id2name(args->hdr,args->rid); int start = args->fa_src_pos - len; int end = args->fa_src_pos; if ( !regidx_overlap(args->mask, chr,start,end, args->itr) ) return; if ( !regidx_overlap(args->mask, args->chr,start,end, args->itr) ) return; int idx_start, idx_end, i; while ( regitr_overlap(args->itr) ) Loading
bcftools/consensus.c.pysam.c +9 −4 Original line number Diff line number Diff line Loading @@ -79,6 +79,7 @@ typedef struct rbuf_t vcf_rbuf; bcf1_t **vcf_buf; int nvcf_buf, rid; char *chr; regidx_t *mask; regitr_t *itr; Loading Loading @@ -123,6 +124,8 @@ static void destroy_chain(args_t *args) free(chain->block_lengths); free(chain); chain = NULL; free(args->chr); args->chr = NULL; } static void print_chain(args_t *args) Loading Loading @@ -164,7 +167,7 @@ static void print_chain(args_t *args) score += chain->block_lengths[n]; } score += last_block_size; fprintf(args->fp_chain, "chain %d %s %d + %d %d %s %d + %d %d %d\n", score, bcf_hdr_id2name(args->hdr,args->rid), ref_end_pos, chain->ori_pos, ref_end_pos, bcf_hdr_id2name(args->hdr,args->rid), alt_end_pos, chain->ori_pos, alt_end_pos, ++args->chain_id); fprintf(args->fp_chain, "chain %d %s %d + %d %d %s %d + %d %d %d\n", score, args->chr, ref_end_pos, chain->ori_pos, ref_end_pos, args->chr, alt_end_pos, chain->ori_pos, alt_end_pos, ++args->chain_id); for (n=0; n<chain->num; n++) { fprintf(args->fp_chain, "%d %d %d\n", chain->block_lengths[n], chain->ref_gaps[n], chain->alt_gaps[n]); } Loading Loading @@ -250,6 +253,7 @@ static void destroy_data(args_t *args) if ( args->vcf_buf[i] ) bcf_destroy1(args->vcf_buf[i]); free(args->vcf_buf); free(args->fa_buf.s); free(args->chr); if ( args->mask ) regidx_destroy(args->mask); if ( args->itr ) regitr_destroy(args->itr); if ( args->chain_fname ) Loading Loading @@ -278,6 +282,7 @@ static void init_region(args_t *args, char *line) else to--; } } args->chr = strdup(line); args->rid = bcf_hdr_name2id(args->hdr,line); if ( args->rid<0 ) fprintf(bcftools_stderr,"Warning: Sequence \"%s\" not in %s\n", line,args->fname); args->fa_buf.l = 0; Loading Loading @@ -382,7 +387,7 @@ static void apply_variant(args_t *args, bcf1_t *rec) if ( regidx_overlap(args->mask, chr,start,end,NULL) ) return; } int i, ialt = 1; int i, ialt = 1; // the alternate allele if ( args->isample >= 0 ) { bcf_unpack(rec, BCF_UN_FMT); Loading Loading @@ -419,6 +424,7 @@ static void apply_variant(args_t *args, bcf1_t *rec) { char ial = rec->d.allele[ialt][0]; char jal = rec->d.allele[jalt][0]; if ( !ialt ) ialt = jalt; // only ialt is used, make sure 0/1 is not ignored rec->d.allele[ialt][0] = gt2iupac(ial,jal); } } Loading Loading @@ -567,11 +573,10 @@ static void apply_variant(args_t *args, bcf1_t *rec) static void mask_region(args_t *args, char *seq, int len) { char *chr = (char*)bcf_hdr_id2name(args->hdr,args->rid); int start = args->fa_src_pos - len; int end = args->fa_src_pos; if ( !regidx_overlap(args->mask, chr,start,end, args->itr) ) return; if ( !regidx_overlap(args->mask, args->chr,start,end, args->itr) ) return; int idx_start, idx_end, i; while ( regitr_overlap(args->itr) ) Loading
bcftools/convert.c +29 −1 Original line number Diff line number Diff line /* convert.c -- functions for converting between VCF/BCF and related formats. Copyright (C) 2013-2017 Genome Research Ltd. Copyright (C) 2013-2018 Genome Research Ltd. Author: Petr Danecek <pd3@sanger.ac.uk> Loading Loading @@ -92,6 +92,7 @@ struct _convert_t int ndat; char *undef_info_tag; int allow_undef_tags; uint8_t **subset_samples; }; typedef struct Loading Loading @@ -174,6 +175,24 @@ static inline int32_t bcf_array_ivalue(void *bcf_array, int type, int idx) } return ((int32_t*)bcf_array)[idx]; } static inline void _copy_field(char *src, uint32_t len, int idx, kstring_t *str) { int n = 0, ibeg = 0; while ( src[ibeg] && ibeg<len && n < idx ) { if ( src[ibeg]==',' ) n++; ibeg++; } if ( ibeg==len ) { kputc('.', str); return; } int iend = ibeg; while ( src[iend] && src[iend]!=',' && iend<len ) iend++; if ( iend>ibeg ) kputsn(src+ibeg, iend-ibeg, str); else kputc('.', str); } static void process_info(convert_t *convert, bcf1_t *line, fmt_t *fmt, int isample, kstring_t *str) { if ( fmt->id<0 ) Loading Loading @@ -232,6 +251,7 @@ static void process_info(convert_t *convert, bcf1_t *line, fmt_t *fmt, int isamp case BCF_BT_INT16: BRANCH(int16_t, val==bcf_int16_missing, val==bcf_int16_vector_end, kputw(val, str)); break; case BCF_BT_INT32: BRANCH(int32_t, val==bcf_int32_missing, val==bcf_int32_vector_end, kputw(val, str)); break; case BCF_BT_FLOAT: BRANCH(float, bcf_float_is_missing(val), bcf_float_is_vector_end(val), kputd(val, str)); break; case BCF_BT_CHAR: _copy_field((char*)info->vptr, info->vptr_len, fmt->subscript, str); break; default: fprintf(stderr,"todo: type %d\n", info->type); exit(1); break; } #undef BRANCH Loading Loading @@ -288,6 +308,8 @@ static void process_format(convert_t *convert, bcf1_t *line, fmt_t *fmt, int isa else kputw(ival, str); } else if ( fmt->fmt->type == BCF_BT_CHAR ) _copy_field((char*)(fmt->fmt->p + isample*fmt->fmt->size), fmt->fmt->size, fmt->subscript, str); else error("TODO: %s:%d .. fmt->type=%d\n", __FILE__,__LINE__, fmt->fmt->type); } else Loading Loading @@ -1312,6 +1334,9 @@ int convert_line(convert_t *convert, bcf1_t *line, kstring_t *str) } for (js=0; js<convert->nsamples; js++) { // Skip samples when filtering was requested if ( *convert->subset_samples && !(*convert->subset_samples)[js] ) continue; // Here comes a hack designed for TBCSQ. When running on large files, // such as 1000GP, there are too many empty fields in the output and // it's very very slow. Therefore in case the handler does not add Loading Loading @@ -1362,6 +1387,9 @@ int convert_set_option(convert_t *convert, enum convert_option opt, ...) case allow_undef_tags: convert->allow_undef_tags = va_arg(args, int); break; case subset_samples: convert->subset_samples = va_arg(args, uint8_t**); break; default: ret = -1; } Loading