Change in the --unify-chr-names option

pd3 · pd3 · commit 2d84e07efdc1 · 2025-03-19T14:09:04.000+01:00
Automatic sequence name modification is not attempted anymore, the prefixes to trim must be given explictly. For example, if run with --unify-chr-names chr,Chromosome,- the program will trim the "chr" prefix in the VCF, "Chromosome" in the GFF, leaving the fasta unchanged Resolves #2378
diff --git a/NEWS b/NEWS
@@ -37,6 +37,11 @@ Changes affecting specific commands:
     - Add experimental support for alternative genetic code tables, accessible via
       a new option `-C, --genetic-code` (#2368)
 
+    - Change in the `--unify-chr-names` option, no automatic sequence name modification
+      is attempted anymore, the prefixes to trim must be given explictly. For example,
+      if run with `--unify-chr-names chr,Chromosome,-`, the program will trim the "chr"
+      prefix in the VCF, "Chromosome" in the GFF, leaving the fasta unchanged (#2378)
+
 * bcftools +fill-tags
 
     - Thanks to the extension of filtering expressions with Fisher's exact test, the plugin
diff --git a/csq.c b/csq.c
diff --git a/doc/bcftools.txt b/doc/bcftools.txt
@@ -1503,10 +1503,12 @@ output VCF and are ignored for the prediction analysis.
 *--targets-overlap* '0'|'1'|'2'::
     see *<<common_options,Common Options>>*
 
-*--unify-chr-names* '0'|'1'::
-    Automatically detect and unify chromosome naming conventions in the GFF, fasta
-    and VCF, such as "chrX" vs "X". The chromosome names in the output VCF will match
-    that of the input VCF. The default is to attempt the automatic translation.
+*--unify-chr-names* '0'|'LIST'::
+    unify chromosome naming by stripping a prefix in VCF, GFF, fasta, respectively.
+    For example, if the LIST is "chr,Chromosome,-", the program will trim the "chr" prefix
+    in the VCF, "Chromosome" in the GFF, and the fasta will be left unchanged.
+    The chromosome names in the output VCF will match that of the input VCF. The default is
+    not to modify the sequence names.
 
 *-v, --verbose* 'INT'::
     verbosity level (0-2)
diff --git a/gff.c b/gff.c
@@ -1,6 +1,6 @@
 /* The MIT License
 
-   Copyright (c) 2023-2024 Genome Research Ltd.
+   Copyright (c) 2023-2025 Genome Research Ltd.
 
    Author: Petr Danecek <pd3@sanger.ac.uk>
 
@@ -88,7 +88,6 @@ typedef struct
     kh_int2tscript_t *id2tr;
 
     // sequences
-    void *seq2int;  // str2int hash
     char **seq;
     int nseq, mseq;
 
@@ -111,13 +110,16 @@ struct gff_t_
     // index iterator
     regidx_t *idx_cds, *idx_utr, *idx_exon, *idx_tscript;
 
+    // str2int hash with parsed sequence names
+    void *seq2int;
+
     // temporary structures, deleted after initializtion
     aux_t init;
 
     // mapping between transcript id (eg. Zm00001d027245_T001) and a numeric idx
     id_tbl_t tscript_ids;
 
-    int strip_chr_names, verbosity;
+    int verbosity;
     int force;      // force run under various conditions. Currently only to skip out-of-phase transcripts
 
     struct {
@@ -158,12 +160,6 @@ int gff_set(gff_t *gff, gff_opt_t key, ...)
             va_end(args);
             return 0;
 
-        case strip_chr_names:
-            va_start(args, key);
-            gff->strip_chr_names = va_arg(args,int);
-            va_end(args);
-            return 0;
-
         case verbosity:
             va_start(args, key);
             gff->verbosity = va_arg(args,int);
@@ -216,12 +212,12 @@ static inline int feature_set_seq(gff_t *gff, char *chr_beg, char *chr_end)
     char tmp = chr_end[1];
     chr_end[1] = 0;
     int iseq;
-    if ( khash_str2int_get(aux->seq2int, chr_beg, &iseq)!=0 )
+    if ( khash_str2int_get(gff->seq2int, chr_beg, &iseq)!=0 )
     {
         char *new_chr = strdup(chr_beg);
         hts_expand(char*, aux->nseq+1, aux->mseq, aux->seq);
         aux->seq[aux->nseq] = new_chr;
-        iseq = khash_str2int_inc(aux->seq2int, aux->seq[aux->nseq]);
+        iseq = khash_str2int_inc(gff->seq2int, aux->seq[aux->nseq]);
         aux->nseq++;
         assert( aux->nseq < 1<<29 );  // see gf_gene_t.iseq and ftr_t.iseq
     }
@@ -239,7 +235,6 @@ static inline void gff_parse_chr(gff_t *gff, const char *line, char **chr_beg, c
     char *se = (char*) line;
     while ( *se && *se!='\t' ) se++;
     if ( !*se ) error("[%s:%d %s] Could not parse the line: %s\n",__FILE__,__LINE__,__FUNCTION__,line);
-    if ( gff->strip_chr_names && !strncasecmp("chr",line,3) ) line += 3;
     *chr_beg = (char*) line;
     *chr_end = se-1;
 }
@@ -974,7 +969,7 @@ int gff_parse(gff_t *gff)
     if ( gff->verbosity > 0 ) fprintf(stderr,"Parsing %s ...\n", gff->fname);
 
     aux_t *aux = &gff->init;
-    aux->seq2int   = khash_str2int_init();   // chrom's numeric id
+    gff->seq2int   = khash_str2int_init();   // chrom's numeric id
     aux->gid2gene  = kh_init(int2gene);      // gene id to gf_gene_t, for idx_gene
     aux->id2tr     = kh_init(int2tscript);   // transcript id to tscript_t
     gff->idx_tscript = regidx_init(NULL, NULL, regidx_free_tscript, sizeof(gf_tscript_t*), NULL);
@@ -1085,7 +1080,6 @@ int gff_parse(gff_t *gff)
 
     free(aux->seq);
     free(aux->ftr);
-    khash_str2int_destroy_free(aux->seq2int);
     // keeping only to destroy the genes at the end: kh_destroy(int2gene,aux->gid2gene);
     kh_destroy(int2tscript,aux->id2tr);
     gff_id_destroy(&aux->gene_ids);
@@ -1119,7 +1113,12 @@ void gff_destroy(gff_t *gff)
     regidx_destroy(gff->idx_exon);
     regidx_destroy(gff->idx_tscript);
 
+    khash_str2int_destroy_free(gff->seq2int);
     gff_id_destroy(&gff->tscript_ids);
     free(gff);
 }
+int gff_has_seq(gff_t *gff, const char *seq)
+{
+    return khash_str2int_has_key(gff->seq2int, seq);
+}
 
diff --git a/gff.h b/gff.h
@@ -1,6 +1,6 @@
 /* The MIT License
 
-   Copyright (c) 2023-2024 Genome Research Ltd.
+   Copyright (c) 2023-2025 Genome Research Ltd.
 
    Author: Petr Danecek <pd3@sanger.ac.uk>
 
@@ -289,7 +289,6 @@ typedef enum
 {
     // write options
     verbosity,          // int, 0-2
-    strip_chr_names,    // int, 0 to leave as is, 1 to strip 'chr' prefix
     force_out_of_phase, // int, 1 to proceed even CDS exon out of expected phase
     dump_fname,         // const char*, dump the parsed GFF into this file, for debugging purposes
 
@@ -314,4 +313,6 @@ void *gff_get(gff_t *gff, gff_opt_t key);
 const char *gff_id2string(gff_t *gff, id_type_t type, int id);
 const char *gf_type2gff_string(int type);
 
+int gff_has_seq(gff_t *gff, const char *chr);
+
 #endif
diff --git a/test/csq.yychr.fa b/test/csq.yychr.fa
@@ -0,0 +1,2 @@
+>Chromosome1
+TTAAGGCTGTTTTTTTATTAATGTCATCGTCCATCCTGCAGGGTTGAACTTGAAAGAATA
diff --git a/test/csq.yychr.gff b/test/csq.yychr.gff
@@ -0,0 +1,3 @@
+Chromosome1	.	gene	1	60	.	+	.	ID=gene:gtrV;Name=gtrV;gene_biotype=protein_coding
+Chromosome1	.	transcript	1	60	.	+	.	ID=transcript:gtrV;Parent=gene:gtrV;gene_biotype=protein_coding
+Chromosome1	.	CDS	1	60	.	+	.	Parent=transcript:gtrV
diff --git a/test/csq.yychr.out b/test/csq.yychr.out
@@ -0,0 +1,2 @@
+.
+missense|gtrV|gtrV|protein_coding|+|1L>1I|1T>A
diff --git a/test/csq.yychr.vcf b/test/csq.yychr.vcf
@@ -0,0 +1,6 @@
+##fileformat=VCFv4.2
+##reference=dummy.fa
+##contig=<ID=Chromosome1,length=1254>
+##FORMAT=<ID=GT,Type=String,Number=1,Description="Genotype">
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	smpl
+Chromosome1	1	.	T	A	.	.	.	GT	0|1
diff --git a/test/test.pl b/test/test.pl
@@ -1015,17 +1015,18 @@
 run_test(\&test_mpileup,$opts,in=>[qw(annot-NMBZ.2)],ref=>'annot-NMBZ.2.fa',out=>'mpileup/annot-NMBZ.2.1.out',args=>q[-a -AD,INFO/NMBZ -r chr6:75]);
 run_test(\&test_mpileup,$opts,in=>[qw(annot-NMBZ.3.1 annot-NMBZ.3.2)],ref=>'annot-NMBZ.3.fa',out=>'mpileup/annot-NMBZ.3.1.out',args=>q[-a -AD,INFO/NMBZ -r chr16:75]);
 run_test(\&test_csq,$opts,in=>'csq',out=>'csq.1.out',cmd=>'-f {PATH}/csq.fa -g {PATH}/csq.gff3');
-run_test(\&test_csq,$opts,in=>'csq',out=>'csq.1.out',cmd=>'-f {PATH}/csq.fa -g {PATH}/csq.chr.gff3');
+run_test(\&test_csq,$opts,in=>'csq',out=>'csq.1.out',cmd=>'-f {PATH}/csq.fa -g {PATH}/csq.chr.gff3 --unify-chr-names -,chr,-');
 run_test(\&test_csq,$opts,in=>'csq.2',out=>'csq.2.out',cmd=>'-f {PATH}/csq.fa -g {PATH}/csq.2.gff',tbcsq=>1);
 run_test(\&test_csq,$opts,in=>'csq.2',out=>'csq.3.out',cmd=>'-f {PATH}/csq.fa -g {PATH}/csq.2.gff --ncsq 64',tbcsq=>1);
-run_test(\&test_csq,$opts,in=>'csq.nchr',out=>'csq.chr.out',cmd=>'-f {PATH}/csq.nchr.fa -g {PATH}/csq.nchr.gff',tbcsq=>1);
-run_test(\&test_csq,$opts,in=>'csq.nchr',out=>'csq.chr.out',cmd=>'-f {PATH}/csq.ychr.fa -g {PATH}/csq.nchr.gff',tbcsq=>1);
-run_test(\&test_csq,$opts,in=>'csq.nchr',out=>'csq.chr.out',cmd=>'-f {PATH}/csq.nchr.fa -g {PATH}/csq.ychr.gff',tbcsq=>1);
-run_test(\&test_csq,$opts,in=>'csq.nchr',out=>'csq.chr.out',cmd=>'-f {PATH}/csq.ychr.fa -g {PATH}/csq.ychr.gff',tbcsq=>1);
-run_test(\&test_csq,$opts,in=>'csq.ychr',out=>'csq.chr.out',cmd=>'-f {PATH}/csq.ychr.fa -g {PATH}/csq.ychr.gff',tbcsq=>1);
-run_test(\&test_csq,$opts,in=>'csq.ychr',out=>'csq.chr.out',cmd=>'-f {PATH}/csq.ychr.fa -g {PATH}/csq.nchr.gff',tbcsq=>1);
-run_test(\&test_csq,$opts,in=>'csq.ychr',out=>'csq.chr.out',cmd=>'-f {PATH}/csq.nchr.fa -g {PATH}/csq.ychr.gff',tbcsq=>1);
-run_test(\&test_csq,$opts,in=>'csq.ychr',out=>'csq.chr.out',cmd=>'-f {PATH}/csq.nchr.fa -g {PATH}/csq.nchr.gff',tbcsq=>1);
+run_test(\&test_csq,$opts,in=>'csq.nchr',out=>'csq.chr.out',cmd=>'-g {PATH}/csq.nchr.gff -f {PATH}/csq.nchr.fa --unify-chr-names -,-,-',tbcsq=>1);
+run_test(\&test_csq,$opts,in=>'csq.nchr',out=>'csq.chr.out',cmd=>'-g {PATH}/csq.nchr.gff -f {PATH}/csq.ychr.fa --unify-chr-names -,-,chr',tbcsq=>1);
+run_test(\&test_csq,$opts,in=>'csq.nchr',out=>'csq.chr.out',cmd=>'-g {PATH}/csq.ychr.gff -f {PATH}/csq.nchr.fa --unify-chr-names -,chr,-',tbcsq=>1);
+run_test(\&test_csq,$opts,in=>'csq.nchr',out=>'csq.chr.out',cmd=>'-g {PATH}/csq.ychr.gff -f {PATH}/csq.ychr.fa --unify-chr-names -,chr,chr',tbcsq=>1);
+run_test(\&test_csq,$opts,in=>'csq.ychr',out=>'csq.chr.out',cmd=>'-g {PATH}/csq.ychr.gff -f {PATH}/csq.ychr.fa --unify-chr-names chr,chr,chr',tbcsq=>1);
+run_test(\&test_csq,$opts,in=>'csq.ychr',out=>'csq.chr.out',cmd=>'-g {PATH}/csq.nchr.gff -f {PATH}/csq.ychr.fa --unify-chr-names chr,-,chr',tbcsq=>1);
+run_test(\&test_csq,$opts,in=>'csq.ychr',out=>'csq.chr.out',cmd=>'-g {PATH}/csq.ychr.gff -f {PATH}/csq.nchr.fa --unify-chr-names chr,chr,-',tbcsq=>1);
+run_test(\&test_csq,$opts,in=>'csq.ychr',out=>'csq.chr.out',cmd=>'-g {PATH}/csq.nchr.gff -f {PATH}/csq.nchr.fa --unify-chr-names chr,-,-',tbcsq=>1);
+run_test(\&test_csq,$opts,in=>'csq.yychr',out=>'csq.yychr.out',cmd=>'-f {PATH}/csq.nchr.fa -g {PATH}/csq.ychr.gff --unify-chr-names Chromosome,chr,-',tbcsq=>1);
 run_test(\&test_csq_real,$opts,in=>'csq');
 run_test(\&test_roh,$opts,in=>'roh.1',out=>'roh.1.1.out',args=>q[-Or -G30 --AF-dflt 0.4]);
 run_test(\&test_roh,$opts,in=>'roh.1',out=>'roh.1.1.out',args=>q[-Or -G30 --AF-file {PATH}/roh.1.tab.gz]);
diff --git a/vcfnorm.c b/vcfnorm.c
@@ -409,7 +409,6 @@ static int is_left_align(args_t *args, bcf1_t *line)
     if ( args->right_align ) return 0;
     if ( !args->gff ) return 1;
     const char *chr = bcf_seqname(args->hdr,line);
-    if ( !strncasecmp("chr",chr,3) ) chr += 3;  // strip 'chr' prefix, that's what we requested the GFF reader to do
     if ( !regidx_overlap(args->idx_tscript,chr,line->pos,line->pos+line->rlen, args->itr_tscript) ) return 1;
 
     // if there are two conflicting overlapping transcripts, go with the default left-alignment
@@ -2248,7 +2247,6 @@ static void init_data(args_t *args)
     {
         args->gff = gff_init(args->gff_fname);
         gff_set(args->gff,verbosity,args->gff_verbosity);
-        gff_set(args->gff,strip_chr_names,1);
         gff_parse(args->gff);
         args->idx_tscript = gff_get(args->gff,idx_tscript);
         args->itr_tscript = regitr_init(NULL);

Original file line number	Diff line number	Diff line change
`@@ -0,0 +1,2 @@`
	`1`	`+>Chromosome1`
	`2`	`+TTAAGGCTGTTTTTTTATTAATGTCATCGTCCATCCTGCAGGGTTGAACTTGAAAGAATA`
Original file line number	Diff line number	Diff line change
`@@ -0,0 +1,3 @@`
	`1`	`+Chromosome1 . gene 1 60 . + . ID=gene:gtrV;Name=gtrV;gene_biotype=protein_coding`
	`2`	`+Chromosome1 . transcript 1 60 . + . ID=transcript:gtrV;Parent=gene:gtrV;gene_biotype=protein_coding`
	`3`	`+Chromosome1 . CDS 1 60 . + . Parent=transcript:gtrV`
Original file line number	Diff line number	Diff line change
`@@ -0,0 +1,2 @@`
	`1`	`+.`
	`2`	`+missense\|gtrV\|gtrV\|protein_coding\|+\|1L>1I\|1T>A`