"Options:\n",
" -r, --refseq <file.fa> The reference sequence, required when indels are present.\n",
" -s, --snps-only Ignore indels.\n",
+ " -t, --column-title <string> The column title.\n",
" -h, -?, --help This help message.\n",
"\n";
}
while (my $arg=shift(@ARGV))
{
if ( $arg eq '-r' || $arg eq '--refseq' ) { $opts{refseq}=shift(@ARGV); next; }
+ if ( $arg eq '-t' || $arg eq '--column-title' ) { $opts{title}=shift(@ARGV); next; }
if ( $arg eq '-s' || $arg eq '--snps-only' ) { $opts{snps_only}=1; next; }
if ( $arg eq '-?' || $arg eq '-h' || $arg eq '--help' ) { error(); }
);
if ( !exists($iupac{$base}) )
{
- if ( $ref eq $base ) { return ('.','0|0'); }
- return ($base,'1|1');
+ if ( $base ne 'A' && $base ne 'C' && $base ne 'G' && $base ne 'T' ) { error("FIXME: what is this [$base]?\n"); }
+ if ( $ref eq $base ) { return ('.','0/0'); }
+ return ($base,'1/1');
}
my $gt = $iupac{$base};
- if ( $$gt[0] eq $ref ) { return ($$gt[1],'0|1'); }
- elsif ( $$gt[1] eq $ref ) { return ($$gt[0],'0|1'); }
- return ("$$gt[0],$$gt[1]",'1|2');
+ if ( $$gt[0] eq $ref ) { return ($$gt[1],'0/1'); }
+ elsif ( $$gt[1] eq $ref ) { return ($$gt[0],'0/1'); }
+ return ("$$gt[0],$$gt[1]",'1/2');
}
my ($prev_chr,$prev_pos,$prev_ref);
my $refseq;
my $ignore_indels = $$opts{snps_only} ? 1 : 0;
+ my $title = exists($$opts{title}) ? $$opts{title} : 'data';
+
+ print $fh_out
+ qq[##fileformat=VCFv3.3\n],
+ qq[##INFO=DP,1,Integer,"Total Depth"\n],
+ qq[##FORMAT=GT,1,String,"Genotype"\n],
+ qq[##FORMAT=GQ,1,Integer,"Genotype Quality"\n],
+ qq[##FORMAT=DP,1,Integer,"Read Depth"\n],
+ qq[#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\t$title\n]
+ ;
while (my $line=<$fh_in>)
{
chomp($line);
- my ($chr,$pos,$ref,$cons,$cons_qual,$snp_qual,$rms_qual,$depth,@items) = split(/\t/,$line);
+ my (@items) = split(/\t/,$line);
+ if ( scalar @items<8 )
+ {
+ error("\nToo few columns, does not look like output of 'samtools pileup -c': $line\n");
+ }
+ my ($chr,$pos,$ref,$cons,$cons_qual,$snp_qual,$rms_qual,$depth) = @items;
my ($alt,$gt);
if ( $ref eq '*' )
if ( !$alt1 )
{
$alt=$alt2;
- $gt='0|1';
+ $gt='0/1';
}
elsif ( !$alt2 )
{
$alt=$alt1;
- $gt='0|1';
+ $gt='0/1';
}
else
{
$alt="$alt1,$alt2";
- $gt='1|2';
+ $gt='1/2';
}
}
else
($alt,$gt) = iupac_to_gtype($ref,$cons);
}
- print $fh_out "$chr\t$pos\t.\t$ref\t$alt\t$snp_qual\t0\t.\tGT:GQ:DP\t$gt:$cons_qual:$depth\n";
+ print $fh_out "$chr\t$pos\t.\t$ref\t$alt\t$snp_qual\t0\tDP=$depth\tGT:GQ:DP\t$gt:$cons_qual:$depth\n";
$prev_ref = $ref;
$prev_pos = $pos;