More ploidy fixes to #2122

Author: pd3

ploidy.c

@@ -1,5 +1,5 @@
 /*
-    Copyright (C) 2014-2016 Genome Research Ltd.
+    Copyright (C) 2014-2025 Genome Research Ltd.
 
     Author: Petr Danecek <[email protected]>
 
@@ -58,7 +58,7 @@ int ploidy_parse(const char *line, char **chr_beg, char **chr_end, uint32_t *beg
     ploidy_t *ploidy = (ploidy_t*) usr;
     void *sex2id = ploidy->sex2id;
 
-    // Check for special case of default ploidy "* * * <sex> <ploidy>"
+    // Check for special case of default ploidy "* * * SEX PLOIDY"
     int default_ploidy_def = 0;
 
     char *ss = (char*) line;
@@ -112,7 +112,7 @@ int ploidy_parse(const char *line, char **chr_beg, char **chr_end, uint32_t *beg
     // Special case, chr="*" stands for a default value
     if ( default_ploidy_def )
     {
-        ploidy->sex2dflt[ploidy->nsex-1] = sp->ploidy;
+        ploidy->sex2dflt[sp->sex] = sp->ploidy;
         return -1;
     }
 
@@ -266,3 +266,23 @@ int ploidy_min(ploidy_t *ploidy)
     return ploidy->dflt < ploidy->min ? ploidy->dflt : ploidy->min;
 }
 
+char *ploidy_format(ploidy_t *ploidy)
+{
+    kstring_t str = {0,0,0};
+
+    regitr_t *itr = regitr_init(ploidy->idx);
+    while ( regitr_loop(itr) )
+    {
+        int id  = regitr_payload(itr,sex_ploidy_t).sex;
+        int pld = regitr_payload(itr,sex_ploidy_t).ploidy;
+        ksprintf(&str,"%s\t%d\t%d\t%s\t%d\n", itr->seq, itr->beg+1, itr->end+1, ploidy->id2sex[id],pld);
+    }
+    regitr_destroy(itr);
+
+    int i;
+    for (i=0; i<ploidy->nsex; i++)
+        ksprintf(&str,"*\t*\t*\t%s\t%d\n", ploidy->id2sex[i],ploidy->sex2dflt[i]);
+
+    return str.s;
+}
+

ploidy.h

@@ -1,5 +1,5 @@
-/* 
-    Copyright (C) 2014-2015 Genome Research Ltd.
+/*
+    Copyright (C) 2014-2025 Genome Research Ltd.
 
     Author: Petr Danecek <[email protected]>
 
@@ -9,10 +9,10 @@
     to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
     copies of the Software, and to permit persons to whom the Software is
     furnished to do so, subject to the following conditions:
-    
+
     The above copyright notice and this permission notice shall be included in
     all copies or substantial portions of the Software.
-    
+
     THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
     IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
     FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
@@ -81,8 +81,8 @@ void ploidy_destroy(ploidy_t *ploidy);
  *  @param seq: chromosome name
  *  @param pos: 0-based position
  *  @param sex2ploidy:  if not NULL, array will be filled with mapping from sex id to ploidy
- *  @param min: if not NULL, minimum encountered encountered will be set
- *  @param max: if not NULL, maximum encountered encountered will be set
+ *  @param min: if not NULL, minimum encountered ploidy will be set
+ *  @param max: if not NULL, maximum encountered ploidy will be set
  *
  *  Returns 1 if the position is listed in the regions or 0 otherwise.
  */
@@ -125,5 +125,8 @@ regidx_t *ploidy_regions(ploidy_t *ploidy);
 int ploidy_max(ploidy_t *ploidy);
 int ploidy_min(ploidy_t *ploidy);
 
+/** Create a parseable ploidy file for debugging. The string must be free()-ed by the caller */
+char *ploidy_format(ploidy_t *ploidy);
+
 #endif
 

test/call-ploidy.1.1.out

@@ -0,0 +1,16 @@
+##fileformat=VCFv4.2
+##FILTER=<ID=PASS,Description="All filters passed">
+##reference=file://ref.fa
+##contig=<ID=chr1>
+##contig=<ID=chr2>
+##contig=<ID=chr3>
+##FORMAT=<ID=PL,Number=G,Type=Integer,Description="List of Phred-scaled genotype likelihoods">
+##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
+##INFO=<ID=AC,Number=A,Type=Integer,Description="Allele count in genotypes for each ALT allele, in the same order as listed">
+##INFO=<ID=AN,Number=1,Type=Integer,Description="Total number of alleles in called genotypes">
+##INFO=<ID=DP4,Number=4,Type=Integer,Description="Number of high-quality ref-forward , ref-reverse, alt-forward and alt-reverse bases">
+##INFO=<ID=MQ,Number=1,Type=Integer,Description="Average mapping quality">
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	smpl1	smpl2
+chr1	621	.	A	.	24.4632	.	AN=2;DP4=11,0,1,0;MQ=41	GT	.	0/0
+chr2	621	.	A	.	27.9838	.	AN=2;DP4=11,0,1,0;MQ=41	GT	0	0
+chr3	621	.	A	.	24.4632	.	AN=2;DP4=11,0,1,0;MQ=41	GT	0/0	.

test/call-ploidy.1.ped

@@ -0,0 +1,2 @@
+Fam1	smpl1	0	0	1
+Fam1	smpl2	0	0	2

test/call-ploidy.1.txt

@@ -0,0 +1,6 @@
+chr1 1 16750 1 0
+chr2 1 16750 1 1
+* * * 1 2
+chr1 1 16750 2 2
+chr2 1 16750 2 1
+* * * 2 0

test/call-ploidy.1.vcf

@@ -0,0 +1,12 @@
+##fileformat=VCFv4.2
+##reference=file://ref.fa
+##contig=<ID=chr1>
+##contig=<ID=chr2>
+##contig=<ID=chr3>
+##INFO=<ID=I16,Number=16,Type=Float,Description="Auxiliary tag used for calling, see description of bcf_callret1_t in bam2bcf.h">
+##INFO=<ID=QS,Number=R,Type=Float,Description="Auxiliary tag used for calling">
+##FORMAT=<ID=PL,Number=G,Type=Integer,Description="List of Phred-scaled genotype likelihoods">
+#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	smpl1	smpl2
+chr1	621	.	A	C,<*>	0	.	I16=11,0,1,0,407,15059,37,1369,440,17600,60,3600,252,6032,25,625;QS=8,1,0	PL	0,3,37,3,37,37	0,3,37,3,37,37
+chr2	621	.	A	C,<*>	0	.	I16=11,0,1,0,407,15059,37,1369,440,17600,60,3600,252,6032,25,625;QS=8,1,0	PL	0,3,37,3,37,37	0,3,37,3,37,37
+chr3	621	.	A	C,<*>	0	.	I16=11,0,1,0,407,15059,37,1369,440,17600,60,3600,252,6032,25,625;QS=8,1,0	PL	0,3,37,3,37,37	0,3,37,3,37,37

test/test.pl

@@ -430,6 +430,7 @@
 run_test(\&test_vcf_head2,$opts,in=>'mpileup.2',out=>'head.1.out',args=>'-s0');
 run_test(\&test_vcf_head2,$opts,in=>'mpileup.2',out=>'head.2.out',args=>'-s1');
 run_test(\&test_vcf_head2,$opts,in=>'mpileup.2',out=>'head.3.out',args=>'-s2 -h2');
+run_test(\&test_vcf_call,$opts,in=>'call-ploidy.1',out=>'call-ploidy.1.1.out',args=>'-m --ploidy-file {PATH}/call-ploidy.1.txt -S {PATH}/call-ploidy.1.ped');
 run_test(\&test_vcf_call,$opts,in=>'mpileup',out=>'mpileup.1.out',args=>'-mv');
 run_test(\&test_vcf_call,$opts,in=>'mpileup',out=>'mpileup.2.out',args=>'-mg0');
 run_test(\&test_vcf_call,$opts,in=>'mpileup',out=>'mpileup.3.out',args=>'-mv -S {PATH}/mpileup.3.samples');

vcfcall.c

@@ -239,7 +239,6 @@ static char **parse_ped_samples(args_t *args, call_t *call, char **vals, int nva
             else if ( !strcmp(sex,"2") && ploidy_sex2id(args->ploidy,fsex)>=0 ) sex = fsex;
             else error("[E::%s] The sex \"%s\" has not been declared in --ploidy/--ploidy-file\n",__func__,sex);
         }
-
         lines = add_sample(name2idx, lines, &nlines, &mlines, col_ends[0]+1, sex, &j);
         if ( strcmp(col_ends[1]+1,"0") && strcmp(col_ends[2]+1,"0") )   // father and mother
         {
@@ -325,10 +324,14 @@ static void set_samples(args_t *args, const char *fn, int is_file)
         while ( *se && !isspace(*se) ) se++;
         if ( se==ss ) { *xptr = x; error("Could not parse: \"%s\"\n", lines[i]); }
 
-        if ( ss[1]==0 && (ss[0]=='0' || ss[0]=='1' || ss[0]=='2') )
-            args->sample2sex[nsmpl] = -1*(ss[0]-'0');
+        char *sex = ss;
+        if ( ploidy_sex2id(args->ploidy,sex)<0 )
+        {
+            if ( sex[1]==0 && (sex[0]=='0' || sex[0]=='1' || sex[0]=='2') ) args->sample2sex[nsmpl] = -1*(sex[0]-'0');
+            else error("[E::%s] The sex \"%s\" has not been declared in --ploidy/--ploidy-file\n",__func__,sex);
+        }
         else
-            args->sample2sex[nsmpl] = ploidy_add_sex(args->ploidy, ss);
+            args->sample2sex[nsmpl] = ploidy_add_sex(args->ploidy,sex);
 
         if ( ismpl!=nsmpl ) map_needed = 1;
         args->samples_map[nsmpl] = ismpl;