Error running example with hg19 genome

[pprieto]

Hi,

I am getting an error at postprocess variants.

The command executed is:

/opt/deepvariant/bin/postprocess_variants     --ref "hg19.fa.gz"     --infile call_variants_output.tfrecord     --outfile "NA12878_S1.chr20.10_10p1mb.bam.vcf"

And the ouput is:

  2018-03-06 11:34:21.456020: I deepvariant/postprocess_variants.cc:87] Read from: call_variants_output.tfrecord
  2018-03-06 11:34:21.457925: I deepvariant/postprocess_variants.cc:96] Done reading: call_variants_output.tfrecord. #entries in single_site_calls = 289
  2018-03-06 11:34:21.457943: I deepvariant/postprocess_variants.cc:100] Total #entries in single_site_calls = 289
  2018-03-06 11:34:21.457949: I deepvariant/postprocess_variants.cc:102] Start SortSingleSiteCalls
  2018-03-06 11:34:21.457957: F deepvariant/core/utils.cc:84] Check failed: pos_in_fasta != contig_name_to_pos_in_fasta.end() Reference name chr20 not in contig info.

Any idea why I cannot change the genome to run the example?

[depristo]

Are you saying you've swapped hg19.fa.gz for our hs37d5.fa.gz fasta? That won't work because the chromosome names are different - 20 in hs37d5 and it sounds like chr20 in hg19.fa.gz. DeepVariant requires the BAM file and the reference genome to have consistent contigs.

[pprieto]

I do not think this is the case. The example BAM file, NA12878_S1.chr20.10_10p1mb.bam, contains reads mapped to chromosome 20 which show chr20 in RNAME for all the reads. This means the fasta file should contain chromosomes named as 'chr'.

Morever, the example fasta, ucsc.hg19.chr20.unittest.fasta, comes from hg19 and it works. Any idea why using the entire genome would make the run fail?

[depristo]

The only real requirement is that the fasta and the BAM need to have a large number of consistent contigs (meaning same length in bp and with the same name) and DeepVariant will process the dataset. If you do samtools view -H on the BAM and cat the hg19.fa.gz.fai file, do these have the same set of contigs with the same lengths?

[pprieto]

Ok, might be doing something wrong. How do you create your .gz.fai file? Is it coming from indexing a fasta file compressed with bgzip or how?

[depristo]

I already have a copy of hs37d5.fa on my machine, so to create a gzipped, indexed version I do:

cp hs37d5.fa /tmp
bgzip /tmp/hs37d5.fa
samtools faidx /tmp/hs37d5.fa.gz
head /tmp/hs37d5.fa.gz.fai
1       249250621       52      60      61
2       243199373       253404903       60      61
3       198022430       500657651       60      61
4       191154276       701980507       60      61
5       180915260       896320740       60      61
6       171115067       1080251307      60      61
7       159138663       1254218344      60      61
8       146364022       1416009371      60      61
9       141213431       1564812846      60      61
10      135534747       1708379889      60      61

[pprieto]

Ok, I think that solves the issue. When indexing over the gz only indexed the first chromosome for some reason in my machine