Hello,
I tried to run vargeno (last commit, 00ee0f0) on the following data: reference, VCFs (provided by 1000genomes), and sample sequenced from NA12878 individual.
I ran some tests on different chromosomes but I always got an output VCF file that contains only the header. I tried to figure out the reasons of this and these are my hypotheses.
I think the main problem is in the index step. Indeed, in my various tests, I obtain
...
[BloomFilter constructBfFromVCF] bit vector: 0/1120000000
...
From here, I started digging in the code and I saw that there could be some problems on how you index the input reference. In more details:
- when you read the VCF file and extract the chromosome name from each line, you add "chr" at its start:
|
if(chr_name[0] != 'c') chr_name = "chr" + chr_name; |
but you don't do the same when you read the headers from the input FASTA file:
For instance, this is a problem if the headers in the FASTA file contain only the chromosome number. I think this problem affects also the way you store information in the .chrlens file of your index:
|
fprintf(chrlens, "%s %lu\n", ref.seqs[i].name, ref.seqs[i].size); |
- when a header in the FASTA file contains the unique identifier for the sequence and also additional information (such as: ">22 dna:chromosome chromosome:GRCh37:22:1:51304566:1"), you consider all the line as unique identifier:
but when you parse the VCF file, you consider as chromosome name only the unique identifier since you get the chromosome number from the first column of the VCF (this should not affect the .chrlens file since you use a different function to read the reference FASTA when you write the .chrlens file)
I tried to solve these two problems by changing some lines in your code but I'm not sure if what I've done is right and enough (I'll anyway open a pull request: it could be a good starting point for you). With my fixes, now the output is not empty anymore.
Moreover, I think that this behaviour occurs also if the input VCF contains the field GT specified in the header and the GT columns (as in the VCFs provided by the 1000genomes project). If I run vargeno index, I obtain:
...
SNP Dictionary
Total k-mers: 0
Unambig k-mers: 0
Ambig unique k-mers: 0
Ambig total k-mers: 0
...
Currently, I solved this problem by removing out from the VCFs the line in the header and the ~2500 columns of the samples. Maybe you could find a better solution to this or maybe you can update the readme accordingly.
Thanks in advance!
Best,
Luca