Skip to content

Commit

Permalink
Fix missed VCF regions
Browse files Browse the repository at this point in the history
Awaits the merge of samtools/htslib#1624

Resolves #1918
  • Loading branch information
pd3 committed May 25, 2023
1 parent 7040c10 commit 4435525
Show file tree
Hide file tree
Showing 5 changed files with 47 additions and 2 deletions.
11 changes: 9 additions & 2 deletions consensus.c
Original file line number Diff line number Diff line change
Expand Up @@ -324,7 +324,7 @@ static void init_region(args_t *args, char *line)
{
char *ss, *se = line;
while ( *se && !isspace(*se) && *se!=':' ) se++;
int from = 0, to = 0;
hts_pos_t from = 0, to = 0;
char tmp = 0, *tmp_ptr = NULL;
if ( *se )
{
Expand Down Expand Up @@ -356,7 +356,14 @@ static void init_region(args_t *args, char *line)
args->fa_frz_mod = -1;
args->fa_case = -1;
args->vcf_rbuf.n = 0;
bcf_sr_seek(args->files,line,args->fa_ori_pos);

kstring_t str = {0,0,0};
if ( from==0 ) from = 1;
if ( to==0 ) to = HTS_POS_MAX;
ksprintf(&str,"%s:%"PRIhts_pos"-%"PRIhts_pos,line,from,to);
bcf_sr_set_regions(args->files,line,0);
free(str.s);

if ( tmp_ptr ) *tmp_ptr = tmp;
fprintf(args->fp_out,">%s%s\n",args->chr_prefix?args->chr_prefix:"",line);
if ( args->chain_fname )
Expand Down
6 changes: 6 additions & 0 deletions test/consensus.21.fa
Original file line number Diff line number Diff line change
@@ -0,0 +1,6 @@
>17
ACGTACGT
>18
ACGTACGT
>19
ACGTACGT
11 changes: 11 additions & 0 deletions test/consensus.21.vcf
Original file line number Diff line number Diff line change
@@ -0,0 +1,11 @@
##fileformat=VCFv4.2
##reference=file://some/path/human_g1k_v37.fasta
##contig=<ID=17>
##contig=<ID=18>
##contig=<ID=19>
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT a
19 2 . C A . . . GT 0/0
19 3 . G C . . . GT 0/0
19 4 . T C . . . GT 0/1
19 5 . A C . . . GT 1/1
6 changes: 6 additions & 0 deletions test/consensus21.1.out
Original file line number Diff line number Diff line change
@@ -0,0 +1,6 @@
>17
ACGTACGT
>18
ACGTACGT
>19
ACGYCCGT
15 changes: 15 additions & 0 deletions test/test.pl
Original file line number Diff line number Diff line change
Expand Up @@ -811,6 +811,7 @@
run_test(\&test_vcf_consensus,$opts,in=>'consensus.20',out=>'consensus20.2.out',fa=>'consensus.20.fa',args=>'');
run_test(\&test_vcf_consensus,$opts,in=>'consensus.20',out=>'consensus20.3.out',fa=>'consensus.20.fa',args=>'-M . -s b');
run_test(\&test_vcf_consensus,$opts,in=>'consensus.20',out=>'consensus20.4.out',fa=>'consensus.20.fa',args=>'-M . -s a');
run_test(\&test_vcf_consensus,$opts,in=>'consensus.21',out=>'consensus21.1.out',fa=>'consensus.21.fa',args=>'');
run_test(\&test_mpileup,$opts,in=>[qw(mpileup.1 mpileup.2 mpileup.3)],out=>'mpileup/mpileup.1.out',args=>q[-r17:100-150],test_list=>1);
run_test(\&test_mpileup,$opts,in=>[qw(mpileup.1 mpileup.2 mpileup.3)],out=>'mpileup/mpileup.2.out',args=>q[-a DP,DV -r17:100-600]); # test files from samtools mpileup test suite
run_test(\&test_mpileup,$opts,in=>[qw(mpileup.1)],out=>'mpileup/mpileup.3.out',args=>q[-B --ff 0x14 -r17:1050-1060]); # test file converted to vcf from samtools mpileup test suite
Expand Down Expand Up @@ -1174,6 +1175,18 @@ sub bgzip_tabix_vcf
my ($opts,$file) = @_;
bgzip_tabix($opts,file=>$file,suffix=>'vcf',args=>'-p vcf');
}
sub bgzip_index_bcf
{
my ($opts,$file) = @_;
if ( !-e "$$opts{tmp}/$file.bcf" or is_file_newer("$$opts{path}/$file.vcf","$$opts{tmp}/$file.bcf") )
{
cmd("$$opts{bin}/bcftools view -Ob $$opts{path}/$file.vcf -o $$opts{tmp}/$file.bcf");
}
if ( !-e "$$opts{tmp}/$file.bcf.csi" or is_file_newer("$$opts{tmp}/$file.bcf","$$opts{tmp}/$file.bcf.csi") )
{
cmd("$$opts{bin}/bcftools index -f $$opts{tmp}/$file.bcf");
}
}
# The tests --------------------------
Expand Down Expand Up @@ -1755,10 +1768,12 @@ sub test_vcf_consensus
{
my ($opts,%args) = @_;
bgzip_tabix_vcf($opts,$args{in});
bgzip_index_bcf($opts,$args{in});
$args{args} =~ s/{PATH}/$$opts{path}/g;
my $mask = $args{mask} ? "-m $$opts{path}/$args{mask}" : '';
my $chain = $args{chain} ? "-c $$opts{tmp}/$args{chain}" : '';
test_cmd($opts,%args,cmd=>"$$opts{bin}/bcftools consensus $$opts{tmp}/$args{in}.vcf.gz -f $$opts{path}/$args{fa} $args{args} $mask $chain");
test_cmd($opts,%args,cmd=>"$$opts{bin}/bcftools consensus $$opts{tmp}/$args{in}.bcf -f $$opts{path}/$args{fa} $args{args} $mask $chain");
}
sub test_vcf_consensus_chain
{
Expand Down

0 comments on commit 4435525

Please sign in to comment.