Missing HPV33 - potential false negative

[seedgeorge]

Description of the bug

I think there is an issue with the slimmed-down kraken2 database for virusbreakend distributed with Oncoanalyser. I expected to find HPV33 in one of my samples, but when I ran the pipeline I ended up with an empty file. I replaced the database with the one suggested by virusbreakend itself, and the tool was able to find my expected virus.

Command used and terminal output

Relevant files

No response

System information

No response

[rhassaine]

Thanks for making the issue George - We had some back and forth on slack regarding this issue and I don't know if you ever got around to running the pipeline with one of the 2 hmf reference genomes and if your virus still isn't being picked up?

[seedgeorge]

I haven't, but here is a more full accounting of my setup:

• Custom genome with viral references appended to the hartwig preferred hg38.
• Indexed and prepared using Oncoanalyser in prepare-references mode.
• Virusbreakend run with the --viralreferences parameter, which adds all reads in the listed contigs to the unmapped reads normally considered for VB.
• I've checked that the reads mapped to HPV33 (there are many) are being extracted by VB and processed along with the others.

This worked fine for the majority of my cases, correctly retrieving HPV16, but failed to pickup HPV33. When I swapped out the OA database for the one supplied by VB itself, HPV is correctly detected and reported by downstream tools.

[rhassaine]

Given the variability in your setup, I think it best to proceed by eliminating potential culprits one by one. Given that you've played with both the reference genome (and forced into the pipeline, in the appropriate fashion) and adjusted the args used in virusbreakend (which is totally fine and should work), we can't be certain of which one of these is going wrong. I think there are 3 steps/things that would be beneficial for us to understand what's going on:

1. I would first try using the included human reference genome without any additional contigs and the existing virusbreakend database and check if you find what you're looking for.

2. I would try the above scenario with the adjusted virusbreakend arg, --viralreferences (if the above scenario doesn't work and/or for the sake of understanding).

3. Sticking to only a single variable component, use your custom ref genome that includes your viral component --> nothing should be reported, given that the HPV33 reads should be mapping to the reference genome entirely.

The database could be easily adjusted with your virus if this proves unfruitful but if there's a bigger issue at hand here (whether in OA, virusbreakend, virusinterpreter or one or more of the aforementioned), we'll be able to catch it and have a permanent fix.

I think this would be a good 1st step to troubleshooting this! Don't hesitate to let me know what your thoughts are :)

[luan-n-nguyen]

Hi George,

I suspect that your custom ref genome is somehow clashing with the virusbreakend database which itself contains a fasta file with viral genomes. I think this suggestion by Rayan would be good to try:

1. I would first try using the included human reference genome without any additional contigs and the existing virusbreakend database and check if you find what you're looking for.

The slim virusbreakend database only removes reference sequences for bacteria, archaea and vectors. All viral sequences remain. We're able to detect HPV-33 our own Hartwig cohort.

I can confirm that the virusbreakend database shipped with OA contains the ref sequences for Human papillomavirus type 33:

lnguyen@dev-aus-nextflow:/opt/resources/oncoanalyser/ref_data/2.0.0--3/37/hmftools/misc/virusbreakend/library$ grep "papillomavirus type 33" viral/library.fna
lnguyen@dev-aus-nextflow:/opt/resources/oncoanalyser/ref_data/2.0.0--3/37/hmftools/misc/virusbreakend/library$ grep "papillomavirus type 33" added/rzW5MDJa1g.fna
>LR862077.1 Human papillomavirus type 33 isolate LNS9453833_HPV33, complete genome
>LR861872.1 Human papillomavirus type 33 isolate LNS3510974_HPV33, complete genome
>EU918766.1 Human papillomavirus type 33 isolate LZcc12-33, complete genome
>HQ537707.1 Human papillomavirus type 33 isolate Qv34189, complete genome
>HQ537706.1 Human papillomavirus type 33 isolate BF131, complete genome
>HQ537705.1 Human papillomavirus type 33 isolate Qv23819, complete genome
>HQ537704.1 Human papillomavirus type 33 isolate Z066, complete genome
>HQ537703.1 Human papillomavirus type 33 isolate IN241890, complete genome
>HQ537702.1 Human papillomavirus type 33 isolate RW136, complete genome
>HQ537701.1 Human papillomavirus type 33 isolate Qv32431, complete genome
>HQ537700.1 Human papillomavirus type 33 isolate INJP06456, complete genome
>HQ537699.1 Human papillomavirus type 33 isolate INCC0137, complete genome
>HQ537698.1 Human papillomavirus type 33 isolate Qv34060, complete genome
>HQ537697.1 Human papillomavirus type 33 isolate AS099, complete genome
>HQ537696.1 Human papillomavirus type 33 isolate Qv06895, complete genome
>HQ537695.1 Human papillomavirus type 33 isolate Qv18347, complete genome
>HQ537694.1 Human papillomavirus type 33 isolate RW926, complete genome
>HQ537693.1 Human papillomavirus type 33 isolate Z84, complete genome
>HQ537692.1 Human papillomavirus type 33 isolate RW702, complete genome
>HQ537691.1 Human papillomavirus type 33 isolate BF266, complete genome
>HQ537690.1 Human papillomavirus type 33 isolate Qv35834, complete genome
>HQ537689.1 Human papillomavirus type 33 isolate Qv22751, complete genome
>HQ537688.1 Human papillomavirus type 33 isolate Qv32494, complete genome

Below, doing the same for the default virusbreakend database from here. Note that differences in contigs present are due to different accession dates for the source databases (kraken, NCBI).

lnguyen@dev-aus-nextflow:/data/downloads/virusbreakenddb_20210401/library$ grep "papillomavirus type 33" viral/library.fna
lnguyen@dev-aus-nextflow:/data/downloads/virusbreakenddb_20210401/library$ grep "papillomavirus type 33" added/bjqKvtb8yE.fna
>LR862077.1 Human papillomavirus type 33 isolate LNS9453833_HPV33, complete genome
>LR861872.1 Human papillomavirus type 33 isolate LNS3510974_HPV33, complete genome
>KU298894.1 Human papillomavirus type 33 isolate 67B.33, complete genome
>KU298893.1 Human papillomavirus type 33 isolate 65C.33, complete genome
>KU298892.1 Human papillomavirus type 33 isolate 65B.33, complete genome
>KU298891.1 Human papillomavirus type 33 isolate 65A.33, complete genome
>KF436865.1 Human papillomavirus type 33 isolate BF375, complete genome
>EU918766.1 Human papillomavirus type 33 isolate LZcc12-33, complete genome
>M12732.1 Human papillomavirus type 33, complete genome
>HQ537707.1 Human papillomavirus type 33 isolate Qv34189, complete genome
>HQ537706.1 Human papillomavirus type 33 isolate BF131, complete genome
>HQ537705.1 Human papillomavirus type 33 isolate Qv23819, complete genome
>HQ537704.1 Human papillomavirus type 33 isolate Z066, complete genome
>HQ537703.1 Human papillomavirus type 33 isolate IN241890, complete genome
>HQ537702.1 Human papillomavirus type 33 isolate RW136, complete genome
>HQ537701.1 Human papillomavirus type 33 isolate Qv32431, complete genome
>HQ537700.1 Human papillomavirus type 33 isolate INJP06456, complete genome
>HQ537699.1 Human papillomavirus type 33 isolate INCC0137, complete genome
>HQ537698.1 Human papillomavirus type 33 isolate Qv34060, complete genome
>HQ537697.1 Human papillomavirus type 33 isolate AS099, complete genome
>HQ537696.1 Human papillomavirus type 33 isolate Qv06895, complete genome
>HQ537695.1 Human papillomavirus type 33 isolate Qv18347, complete genome
>HQ537694.1 Human papillomavirus type 33 isolate RW926, complete genome
>HQ537693.1 Human papillomavirus type 33 isolate Z84, complete genome
>HQ537692.1 Human papillomavirus type 33 isolate RW702, complete genome
>HQ537691.1 Human papillomavirus type 33 isolate BF266, complete genome
>HQ537690.1 Human papillomavirus type 33 isolate Qv35834, complete genome
>HQ537689.1 Human papillomavirus type 33 isolate Qv22751, complete genome
>HQ537688.1 Human papillomavirus type 33 isolate Qv32494, complete genome

[seedgeorge]

Reasonable thoughts. I'm running a 'vanilla' OA run now, using all default genome and database, for comparison with my setup.

Note, when I run:
grep M12732 seqid2taxid.map
On the two database seqid files, I receive conflicting results, hence my initial suspicions. M12732 is an HPV reference.

[seedgeorge]

Update:

• with a vanilla run, standard reference genome with no additional viral contigs in the alignment, OA-provided kraken2 database: no HPV33
• with custom reference genome with additional viral contigs in the alignment, and these provided to virusbreakend as a reference file, OA-provided kraken2 database: no HPV33
• with custom reference genome with additional viral contigs in the alignment, and these provided to virusbreakend as a reference file, virusbreakend-provided kraken2 database: HPV33 detected

From my perspective, this is solved but I appreciate you may want to dig more.

[ohofmann]

We'll have a look at some of our samples with confirmed HPV33 presence as well and report back.

[luan-n-nguyen]

Thanks George, i'll take a deeper look at why this contig is missing from the current OA virusbreakend db:
>M12732.1 Human papillomavirus type 33, complete genome