Clinical-Genomics / scout

VCF visualization interface
https://clinical-genomics.github.io/scout
BSD 3-Clause "New" or "Revised" License
150 stars 46 forks source link

No reads shown for MT variants in internal IGV #4866

Closed Jakob37 closed 1 week ago

Jakob37 commented 1 week ago

Describe the bug

Previously (before updating 4.81 -> 4.88), when navigating to an MT variant page, it was possible to bring up the IGV viewer track including coverage tracks.

Now coverage tracks are only present for non MT variants.

A working variant:

igv_working

A non working MT-variant:

igv_m_not_working

Accessed from this "IGV viewer" button:

case_view

Additional info

I saw the IGV JS was updated in 4.88.1, so I tried running on that version. Same situation there - no MT tracks in neither 4.88 nor 4.88.1

A recurrent issue is that we here in Lund call the mitochondrial chromosome "M" instead of "MT". Might it be something related to that?

northwestwitch commented 1 week ago

Hi @Jakob37, where are your MT variants read from, same alignment file as the other variants or are they on a separate file (bam/cram) ?

dnil commented 1 week ago

Right, always interesting with reference genomes. It looks like you are getting the IGV view and the correct locus, but the alignment is not showing. Is that correct?

Jakob37 commented 1 week ago

Hi @Jakob37, where are your MT variants read from, same alignment file as the other variants or are they on a separate file (bam/cram) ?

They are stored in a separate bam file, here is how a sample can look in the case yaml (before adding owner):

---
owner:
family: 'hg002'
lims_id: '0000-00'
samples:
  - analysis_type: wgs
    sample_id: 'hg002'
    sample_name: 'hg002'
    mother: '0'
    father: '0'
    capture_kit:
    phenotype: affected
    sex: male
    bam_path: /access/wgs/bam/hg002_dedup.bam
    mt_bam: /access/wgs/bam/hg002_mito.bam
    d4_file: /access/wgs/cov/hg002_coverage.d4

Right, always interesting with reference genomes. It looks like you are getting the IGV view and the correct locus, but the alignment is not showing. Is that correct?

Yes, that is how I understand it!

Jakob37 commented 1 week ago

OK, actually better I double check our MT BAMs before I say anything. I'll check ...

northwestwitch commented 1 week ago

OK, actually better I double check our MT BAMs before I say anything. I'll check ...

We touched that code quite recently so I think it might be a bug

Jakob37 commented 1 week ago

OK! 🐛

Yes, looking at the MT BAM from the screenshots above, there are reads present.

northwestwitch commented 1 week ago

I can look into it!

dnil commented 1 week ago

It's likely not the chromosome naming or reference genome as such: we are using the same MT genome sequence with the good old 1kg hg19 (displayed on hg38) - more likely how we special case to make the tracks for it.

northwestwitch commented 1 week ago

@Jakob37 any specific error if you right click on the igv page and choose "inspect"?

northwestwitch commented 1 week ago

Other question: are you navigating to the MT variant from another variant (not on MT) or from the IGV browser created from the link on the case page?

Jakob37 commented 1 week ago

I'll check!

Jakob37 commented 1 week ago

Some warnings in the console, but these are also there for the non MT page:

igv_tracks

Other question: are you navigating to the MT variant from another variant (not on MT) or from the IGV browser created from the link on the case page?

I am navigating: case -> Clinical SNV and INDELs -> (finding an MT variant in the list) -> IGV viewer button

northwestwitch commented 1 week ago

Ok, so the javascript shows no errors, just warnings. Thanks!

northwestwitch commented 1 week ago

I've tried to reproduce locally by setting the genome build of the demo case to 38 (your case has build 38, right? ) and modifying the VCF file of the SNV variants so that MT chrom is now M. The track shows .. Mmmm or MTTT :)

Jakob37 commented 1 week ago

your case has build 38, right?

Yes, we only use build 38

and modifying the VCF file of the SNV variants so that MT chrom is now M. The track shows .. Mmmm or MTTT :)

;)

Strange! I'll check if I see the same if loading our hg002 case locally.

dnil commented 1 week ago

Also possibly worth trying out Scout 4.88.1 ie IGV.js 3.0.5 instead of 3.0.4. It shouldn't be that, but if you have like very deep mt bams and not subsampled I guess it might possibly be that region load termination thing anyway?

dnil commented 1 week ago

I've tried to reproduce locally by setting the genome build of the demo case to 38 (your case has build 38, right? ) and modifying the VCF file of the SNV variants so that MT chrom is now M. The track shows .. Mmmm or MTTT :)

@northwestwitch did you also try changing the aln to have chrM instead of MT? The demo aln looks like this: 😊

ST-E00269:73:HLNFFCCXX:2:1103:25499:65669       163     MT      1       60      13S138M =       229     376     AGACATCACGATGGATCACAGGTCTATCACCCTATTAACCACTCACGGGAGCTCTCCATGCATTTGGTATTTTCGTCTGGGGGGTATGCACGCGATAGCATTGCGAGACGCTGGAGCCGGAGCACCCTATGTCGCAGTATCTGTCTTTGAT      +????????I?I???II?I????I??II?I???IIII+???I?I?????????I?I??I???III???III?III?I???5????II???II?I?I5????I??I???II???'???????+?+????II??II????55????I?II??I MD:Z:138        PG:Z:MarkDuplicates          RG:Z:ADM1464A1.171015_HLNFFCCXX-11_XXXXXX.lane2 NM:i:0  AS:i:138        XS:i:86
ST-E00269:73:HLNFFCCXX:2:2113:11505:16041       99      MT      1       60      3S148M  =       142     292     ATGGATCACAGGTCTATCACCCTATTAACCACTCACGGGAGCTCTCCATGCATTTGGTATTTTCGTCTGGGGGGTATGCACGCGATAGCATTGCGAGACGCTGGAGCCGGAGCACCCTATGTCGCAGTATCTGTCTTTGATTCCTGCCTCA      ?????II?I????I?III?I???IIIIII??I?I?I???????I?I??I???III???IIIIII????????????I???I????II???II??????I???????????????I???II??I?????III???I?II??III???????? MD:Z:148        PG:Z:MarkDuplicates          RG:Z:ADM1464A1.171015_HLNFFCCXX-21_XXXXXX.lane2 NM:i:0  AS:i:148        XS:i:83
Jakob37 commented 1 week ago

Also possibly worth trying out Scout 4.88.1 ie IGV.js 3.0.5 instead of 3.0.4. It shouldn't be that, but if you have like very deep mt bams and not subsampled I guess it might possibly be that region load termination thing anyway?

Hmm, yes, I am testing things out with v4.88.1. For some variants I get BAM file(s) missing. While for others I get the very same situation - i.e. the button is there but no reads.

Some debug printing of the values when the button is there of the data sent to the template:

scout_dev  | variant.chromosome: M
scout_dev  | case.mt_bams True

if you have like very deep mt bams and not subsampled I guess it might possibly be that region load termination thing anyway?

I haven't been in contact with that before. Is that something inside the JS-part?

(I'll continue after 🥗 )

dnil commented 1 week ago

Oh, so your chromosome is a plain M? Not chrM? Here is the entry from the hg38_alias.tab from IGV:

chrM    MT  J01415.2    NC_012920.1 chrMT

I guess it could be that then!

Jakob37 commented 1 week ago

Oh, so your chromosome is a plain M?

Yes, that is correct! Not sure about the origin story, but now it is there by inertia.

Is the solution to rename our MTs? Or can IGV be configured somehow? (Guess something must have changed, as the issue seems to have appeared now)

northwestwitch commented 1 week ago

Is the solution to rename our MTs? Or can IGV be configured somehow? (Guess something must have changed, as the issue seems to have appeared now)

No no, no change of VCF files. If the problem is indeed the alias tab files, why does the demo with only M works locally? Are you also providing an index for the bam files? because the track is not added if you are missing the index, check the code here: https://github.com/Clinical-Genomics/scout/blob/0bc179bdc8e26ffa65ad463eb4452a25b6fc81e1/scout/server/blueprints/alignviewers/controllers.py#L315

dnil commented 1 week ago

Tentatively it is not the aln chromosome name, at least I could not reproduce it with the demo.

EDIT: just the GUI fooling me. I should have used the CLI also to change the case. It reproduces with M in the aln for one of the individuals. Screenshot 2024-09-18 at 12 43 58'

Screenshot 2024-09-18 at 12 55 29 Screenshot 2024-09-18 at 12 55 14
dnil commented 1 week ago

No no, no change of VCF files. If the problem is indeed the alias tab files, why does the demo with only M works locally?

Did you really also try with M in the aln though? 😉 But you are right, I probably botched my attempt. It should not work. Let me recheck.

dnil commented 1 week ago

Oh, so your chromosome is a plain M?

Yes, that is correct! Not sure about the origin story, but now it is there by inertia.

Is the solution to rename our MTs? Or can IGV be configured somehow? (Guess something must have changed, as the issue seems to have appeared now)

Hm, consider asking IGV to include it in the alias file? That would be an instant fix. Some backstory and official reference claims might help, but it might be enough just that it seems very very unlikely someone would ever name some other human genome contig in their hg38 reference M. Alternatively, copy their file, add an "M" and change the URL to somewhere you control, like a static served by scout. You can find the URL to change here: https://github.com/Clinical-Genomics/scout/blob/0bc179bdc8e26ffa65ad463eb4452a25b6fc81e1/scout/constants/igv_tracks.py#L11 Or revert to IGV.js <3 I suppose, but you don't want that.

dnil commented 1 week ago

If we decide to approach IGV.js about it it might help to note that either MT or M in the bam just works anyway on desktop IGV. Which I guess has an additional fallback for the chr naming, much like IGV.js did before.

Screenshot 2024-09-18 at 12 41 05
northwestwitch commented 1 week ago

Did you really also try with M in the aln though?

Right, I did not modify the alignment. It makes sense. We could try to open an issue on igv.js and see if they fix it, otherwise we can include the alias tab files in scout and serve them from the app (bit more cumbersome)

northwestwitch commented 1 week ago

Here --> https://github.com/igvteam/igv.js/issues/1896

Jakob37 commented 1 week ago

Nice job hunting this down 💪 Hope they are fine with updating it.