Concordance of resistance predictions with Illumina

mbhall88 commented 3 years ago

Where there is discordance with Illumina, what kind of errors do we see? Missing resistance or false resistance? Are any particular SNPs repeatedly wrong? Table showing Very Major Error (missed resistance), Major Error (missed susceptible), PPV (what % of R calls are R) and NPV (what % of S calls are S) for each tool.

iqbal-lab commented 3 years ago

A sub-question is - can you do better at this than mykrobe

mbhall88 commented 3 years ago

Working through the implementation of this, here is my first-pass process.

For each sample, we take the Illumina predictions (true) and the Nanopore predictions (test) and produce a CSV with classifications for each drug.

Classification

TP: Illumina says resistant and Nanopore says resistant
TN: Illumina says susceptible and Nanopore says susceptible FP: Illumina says susceptible and Nanopore says resistant (Major Error) FN: Illumina says resistant and Nanopore says susceptible (Very Major Error)

Example

Example output from this process

drug,classification,true_call,test_call
Ofloxacin,TN,S,S
Moxifloxacin,TN,S,S
Isoniazid,FP,S,R
Kanamycin,TN,S,S
Ethambutol,TN,S,S
Streptomycin,FN,R,S
Ciprofloxacin,TN,S,S
Pyrazinamide,TP,R,R
Rifampicin,TN,S,S
Amikacin,TN,S,S
Capreomycin,TN,S,S

One main question I have is how to deal with minor resistant calls ("r" as opposed to "R")? Currently, I am treating them as "R". Not sure if we should deal with them differently?

iqbal-lab commented 3 years ago

This looks right to me. Mykrobe doesn't make r calls for nanopore, but it does for illumina, so they might be there in the truth. I'd treat them as R

mbhall88 commented 3 years ago

Very early plot for mykrobe concordance classification counts. Columns (subplots) represent the different classifications (i.e. TP, FN etc.), Y-axis is drug, X-axis is count. Reminder: there are 150 samples.

mbhall88 commented 3 years ago

Sensitivity and specificity for mykrobe concordance

iqbal-lab commented 3 years ago

This https://github.com/mbhall88/head_to_head_pipeline/issues/75#issuecomment-820144138

looks to me like the biggest problem for mykrobe is lots of FNs, which i find believable. Mykrobe needs a lot of coverage

mbhall88 commented 3 years ago

Well, interestingly, it seems that coverage has a somewhat inverse effect on FNs (sensitivity)....

mbhall88 commented 3 years ago

Table showing Very Major Error (missed resistance), Major Error (missed susceptible), PPV (what % of R calls are R) and NPV (what % of S calls are S) for each tool.

	drug	tool	NPV	PPV	sensitivity	specificity	ME/FP	VME/FN
0	Amikacin	mykrobe	0.964539	0.777778	0.583333	0.985507	2	5
1	Capreomycin	mykrobe	0.964539	0.777778	0.583333	0.985507	2	5
2	Ciprofloxacin	mykrobe	0.970803	1	0.764706	1	0	4
3	Ethambutol	mykrobe	0.696296	1	0.267857	1	0	41
4	Isoniazid	mykrobe	0.821429	0.984848	0.8125	0.985714	1	15
5	Kanamycin	mykrobe	0.964286	0.8	0.615385	0.985401	2	5
6	Moxifloxacin	mykrobe	0.970803	1	0.764706	1	0	4
7	Ofloxacin	mykrobe	0.970803	1	0.764706	1	0	4
8	Pyrazinamide	mykrobe	0.901515	1	0.580645	1	0	13
9	Rifampicin	mykrobe	0.727273	1	0.653846	1	0	27
10	Streptomycin	mykrobe	0.824	0.96	0.521739	0.990385	1	22

I'm not sure if it makes sense to aggregate these metrics/counts over drugs though?

mbhall88 commented 3 years ago

So the above plots and tables are for mykrobe with default settings (with --ont for nanopore).

I reran mykrobe and changed the expected error rate to 0.08 (from 0.15) and I get the following results

tl;dr: we lost some specificity (i.e. more FPs; very slight decrease), but massively gain sensitivity (i.e. less FNs)

	drug	tool	NPV	PPV	sensitivity	specificity	ME/FP	VME/FN
0	Amikacin	mykrobe	1	0.857143	1	0.985507	2	0
1	Capreomycin	mykrobe	1	0.857143	1	0.985507	2	0
2	Ciprofloxacin	mykrobe	0.992537	1	0.941176	1	0	1
3	Ethambutol	mykrobe	0.979167	1	0.964286	1	0	2
4	Isoniazid	mykrobe	1	0.91954	1	0.9	7	0
5	Kanamycin	mykrobe	1	0.866667	1	0.985401	2	0
6	Moxifloxacin	mykrobe	0.992537	1	0.941176	1	0	1
7	Ofloxacin	mykrobe	1	1	1	1	0	0
8	Pyrazinamide	mykrobe	0.991667	1	0.967742	1	0	1
9	Rifampicin	mykrobe	1	1	1	1	0	0
10	Streptomycin	mykrobe	0.980769	0.956522	0.956522	0.980769	2	2

Again, somewhat counter-intuitively, coverage seems to have a slight detrimental effect on performance. I'm guessing this is related to https://github.com/Mykrobe-tools/mykrobe/issues/121 This is probably most clearly illustrated with

Apologies for all the different plots. I'm still trying to figure out the best way to summarise this info/relationship.

iqbal-lab commented 3 years ago

Wow, that is a huge improvement in sensitivity by changing the error rate, congratulations and thanks for trying that

iqbal-lab commented 3 years ago

I agree best not to aggregate across drugs

mbhall88 commented 3 years ago

First glimpse at `drprg` concordance with mykrobe Illumina

This is the first results for concordance of nanopore data to Illumina-Mykrobe results.

Caveats:

Not yet handling any of the situations where a frameshift causes resistance - i.e. katG and pncA
not doing any de novo variant calling
using filtering from clustering, which is more aimed at increasing precision rather than recall

	drug	tool	NPV	PPV	sensitivity	specificity	ME/FP	VME/FN	TP	TN
0	Amikacin	drprg	0.92	0	0	1	0	12	0	138
1	Amikacin	mykrobe	1	0.857143	1	0.985507	2	0	12	136
2	Capreomycin	drprg	0.92	0	0	1	0	12	0	138
3	Capreomycin	mykrobe	1	0.857143	1	0.985507	2	0	12	136
4	Ciprofloxacin	drprg	0.886667	0	0	1	0	17	0	133
5	Ciprofloxacin	mykrobe	0.992537	1	0.941176	1	0	1	16	133
6	Ethambutol	drprg	0.740157	1	0.410714	1	0	33	23	94
7	Ethambutol	mykrobe	0.979167	1	0.964286	1	0	2	54	94
8	Isoniazid	drprg	0.907895	0.986486	0.9125	0.985714	1	7	73	69
9	Isoniazid	mykrobe	1	0.91954	1	0.9	7	0	80	63
10	Kanamycin	drprg	0.919463	1	0.0769231	1	0	12	1	137
11	Kanamycin	mykrobe	1	0.866667	1	0.985401	2	0	13	135
12	Moxifloxacin	drprg	0.886667	0	0	1	0	17	0	133
13	Moxifloxacin	mykrobe	0.992537	1	0.941176	1	0	1	16	133
14	Ofloxacin	drprg	0.886667	0	0	1	0	17	0	133
15	Ofloxacin	mykrobe	1	1	1	1	0	0	17	133
16	Pyrazinamide	drprg	0.793333	0	0	1	0	31	0	119
17	Pyrazinamide	mykrobe	0.991667	1	0.967742	1	0	1	30	119
18	Rifampicin	drprg	0.48	0	0	1	0	78	0	72
19	Rifampicin	mykrobe	1	1	1	1	0	0	78	72
20	Streptomycin	drprg	0.87395	1	0.673913	1	0	15	31	104
21	Streptomycin	mykrobe	0.980769	0.956522	0.956522	0.980769	2	2	44	102

mbhall88 commented 3 years ago

I've spent all of today manually looking at (drprg) FNs for a few samples and the overwhelming reason for most of them is the FRS filter. If we were to remove that filter, about 95% of the FNs I looked at would become TPs as we made the correct calls, but the FRS filter was applied so we ignored the R call when producing the output JSON.

The reason the FRS filter is not suited to this situation is probably best illustrated with the following (drprg/pandora) VCF record

rrs     1501    bca790cf        AC      AA,AG,AT,CC,GC,TC       .       frs     SVTYPE=PH_SNPs;GRAPHTYPE=SIMPLE;VARID=rrs_A1401X,rrs_C1402X;PREDICT=R,S GT:MEAN_FWD_COVG:MEAN_REV_COVG:MED_FWD_COVG:MED_REV_COVG:SUM_FWD_COVG:SUM_REV_COVG:GAPS:LIKELIHOOD:GT_CONF    5:7,11,11,11,5,26,11:10,15,15,15,7,37,15:0,11,11,11,0,27,
11:0,15,15,15,0,36,15:23,23,23,23,23,79,23:31,31,31,31,31,112,31:0.666667,0.5,0.5,0.5,0.75,0,0.5:-899.432,-836.433,-836.433,-836.433,-935.704,-615.668,-836.433:220.765

This shows us that this variant covers rrs_A1401X and rrs_C1402X and is deemed 'R' for rrs_A1401X (which agrees with mykrobe Illumina result). The alt alleles that actually describe rrs_A1401X are 4, 5, and 6. So whilst we would expect/hope that only one of those would get coverage, in reality we get coverage on all three, which ends up skewing the FRS quite considerably.

I will test out removing the FRS filter and see what impact that has on the results.

iqbal-lab commented 3 years ago

So, we actually genotype correctly, but frs is triggered? I see we genotype as 5, that's right, right? Is the problem that our implementation of frs is wrong, our conception of it is wrong for complex sites, or something else?

iqbal-lab commented 3 years ago

Sorry, I get it now. Looks to me like noise on wrong alleles is triggering frs?

iqbal-lab commented 3 years ago

I guess quite a lot of random snp errors in genes will look like X alleles

iqbal-lab commented 3 years ago

What I mean is, is it a problem that some if the wrong DNA alleles happen to trigger the same, correct, amino allele?

mbhall88 commented 3 years ago

So, we actually genotype correctly, but frs is triggered? I see we genotype as 5, that's right, right?

Yes and yes

Is the problem that our implementation of frs is wrong, our conception of it is wrong for complex sites, or something else?

What I mean is, is it a problem that some if the wrong DNA alleles happen to trigger the same, correct, amino allele?

I think it is a combination of

our conception is wrong for complex sites
complex sites are slightly different to our mental model. What I mean by this is if we have two codons together in a pandora site, we don't get an allele for the every possible combination of all codons (the above variant illustrates this) which will also cause even more problems if both are present at the same time
it does worry me a little that pandora assigns so much coverage to "incorrect" alleles

iqbal-lab commented 3 years ago

But at the DNA level there's nothing wrong with the model? And frs acts at that level? I'm confused. I can see how we might need to drop frs threshold die to nanopore error rate, but I don't get what cod9ns have to do with it.

mbhall88 commented 3 years ago

But at the DNA level there's nothing wrong with the model?

Everything I've talked about so far is at the DNA level... Not sure I get what you're asking...

I can see how we might need to drop frs threshold die to nanopore error rate

I'm not convinced it's the ONT error rate that is the reason for needing to drop FRS. The site in this example is pretty "clean" in the bcftools VCF but we get (quite a bit of) coverage on all of the alleles in pandora.

I don't get what cod9ns have to do with it.

Just substitute "codons" with "SNPs" then - it doesn't really matter - I'm just trying to say two panel variants in the same pandora VCF record.

mbhall88 commented 3 years ago

Results after removing FRS filter

So that gets rid of a lot of FNs. I'll manually look into the remaining FNs in the coming days, but it looks like they're probably just the frame shifts we didn't add in. I'll also start digging into the FPs.

One small win is we have slightly more sensitivity/recall for Streptomycin than mykrobe.

	drug	tool	NPV	PPV	sensitivity	specificity	ME/FP	VME/FN	TP	TN
0	Amikacin	drprg	1	0.857143	1	0.985507	2	0	12	136
1	Amikacin	mykrobe	1	0.857143	1	0.985507	2	0	12	136
2	Capreomycin	drprg	1	0.857143	1	0.985507	2	0	12	136
3	Capreomycin	mykrobe	1	0.857143	1	0.985507	2	0	12	136
4	Ciprofloxacin	drprg	0.886667	0	0	1	0	17	0	133
5	Ciprofloxacin	mykrobe	0.992537	1	0.941176	1	0	1	16	133
6	Ethambutol	drprg	0.979167	1	0.964286	1	0	2	54	94
7	Ethambutol	mykrobe	0.979167	1	0.964286	1	0	2	54	94
8	Isoniazid	drprg	0.969697	0.928571	0.975	0.914286	6	2	78	64
9	Isoniazid	mykrobe	1	0.91954	1	0.9	7	0	80	63
10	Kanamycin	drprg	1	0.866667	1	0.985401	2	0	13	135
11	Kanamycin	mykrobe	1	0.866667	1	0.985401	2	0	13	135
12	Moxifloxacin	drprg	0.886667	0	0	1	0	17	0	133
13	Moxifloxacin	mykrobe	0.992537	1	0.941176	1	0	1	16	133
14	Ofloxacin	drprg	0.892617	1	0.0588235	1	0	16	1	133
15	Ofloxacin	mykrobe	1	1	1	1	0	0	17	133
16	Pyrazinamide	drprg	0.933333	0.766667	0.741935	0.941176	7	8	23	112
17	Pyrazinamide	mykrobe	0.991667	1	0.967742	1	0	1	30	119
18	Rifampicin	drprg	1	1	1	1	0	0	78	72
19	Rifampicin	mykrobe	1	1	1	1	0	0	78	72
20	Streptomycin	drprg	0.990291	0.957447	0.978261	0.980769	2	1	45	102
21	Streptomycin	mykrobe	0.980769	0.956522	0.956522	0.980769	2	2	44	102

mbhall88 commented 3 years ago

I have look at every single FP and FN to get an idea of where we are making mistakes (considering mykrobe as the truth).

I won't like the results of every investigation here, but will try to summarise general themes

FNs

7 / 29 are not in the panel - i.e frame-shift mutations in katG or pncA
13 / 29 are due to a natural polymorphism missing from the PRG. The resistance variant we are missing is gyrA_D94X. We get 0 coverage on this VCF position because there is a very common SNP in the population gyrA_S95T (known to be silent) that these samples have and is not in our PRG. Because this silent SNP is directly after the variant of interest, we basically don't get any minimizer hits. Indeed, when I look at samples that don't have this silent SNP, they do have coverage on this VCF position.
2 / 29 relate to Ethambutol variant embB_M306[V,L]. These two are called minor resistance by mykrobe Illumina and susceptible by mykrobe nanopore. They were filtered out in the COMPASS and bcftools VCFs for HET/FRS reasons. We had low-ish coverage on the ALT. All in all this was a tough one.
2 / 29 was Pyrazinamide variant pncA_H137P and pncA_Q141P. This site in the PRG has a huge number of alleles due to lots of variants close together. There was coverage on all of the alleles so we didn't really have a chance. Interestingly, we get this variant correct in the pandora clustering VCF (i.e with the population PRG). The sample that missed the Q141P variant was a minor resistance call by mykrobe
3 / 29 fluoroquinolone variant gyrA_A90V. The drprg VCF has it as a null call. The weird thing about this one is that the pandora consensus VCF from the drprg run actually had the correct GT for these 3. This is one of the next things I am looking into as I was under the assumption the consensus should be REF for all positions as I am specifying a VCF ref
1 / 29 Streptomycin variant rpsL_K43R. Minor resistance call by mykrobe Illumina and susceptible by mykrobe ONT. It seems like a fairly clear ref call in the pandora VCF. COMPASS called HET and bcftools was filtered due to low QUAL. So it is hard to say whether this is a real FN without pheno data, which we don't have for Strep, for this sample.

FPs

There were 13 in total. Quite a few of them seem to solvable by raising the GT CONF filter to 15 for ONT and 30 for Illumina. However, 5 / 13 have a phenotype of 'R', meaning it is actually a FN for mykrobe and a TP for us.

mbhall88 commented 3 years ago

Adding in the GT CONF filter is a bit of a mixed bag

	drug	tool	NPV	PPV	sensitivity	specificity	ME/FP	VME/FN	TP	TN
0	Amikacin	drprg	1	0.857143	1	0.985507	2	0	12	136
1	Amikacin	mykrobe	1	0.857143	1	0.985507	2	0	12	136
2	Capreomycin	drprg	1	0.857143	1	0.985507	2	0	12	136
3	Capreomycin	mykrobe	1	0.857143	1	0.985507	2	0	12	136
4	Ciprofloxacin	drprg	0.886667	0	0	1	0	17	0	133
5	Ciprofloxacin	mykrobe	0.992537	1	0.941176	1	0	1	16	133
6	Ethambutol	drprg	0.979167	1	0.964286	1	0	2	54	94
7	Ethambutol	mykrobe	0.979167	1	0.964286	1	0	2	54	94
8	Isoniazid	drprg	0.985714	0.9875	0.9875	0.985714	1	1	79	69
9	Isoniazid	mykrobe	1	0.91954	1	0.9	7	0	80	63
10	Kanamycin	drprg	1	0.866667	1	0.985401	2	0	13	135
11	Kanamycin	mykrobe	1	0.866667	1	0.985401	2	0	13	135
12	Moxifloxacin	drprg	0.886667	0	0	1	0	17	0	133
13	Moxifloxacin	mykrobe	0.992537	1	0.941176	1	0	1	16	133
14	Ofloxacin	drprg	0.892617	1	0.0588235	1	0	16	1	133
15	Ofloxacin	mykrobe	1	1	1	1	0	0	17	133
16	Pyrazinamide	drprg	0.883721	0.761905	0.516129	0.957983	5	15	16	114
17	Pyrazinamide	mykrobe	0.991667	1	0.967742	1	0	1	30	119
18	Rifampicin	drprg	0.9	1	0.897436	1	0	8	70	72
19	Rifampicin	mykrobe	1	1	1	1	0	0	78	72
20	Streptomycin	drprg	0.990291	0.957447	0.978261	0.980769	2	1	45	102
21	Streptomycin	mykrobe	0.980769	0.956522	0.956522	0.980769	2	2	44	102

we improve on FPs for Isoniazid, but do worse for FNs in Pyrazinamide and Rifampicin

iqbal-lab commented 3 years ago

So you can recover 20/29 FNs out of the box by fixing catalog and background mutations , is good news. Still absorbing FP info

mbhall88 commented 3 years ago

Table 2 Comparison of Nanopore drug resistance predictions with Illumina predictions

For this comparison, we assume the mykrobe resistance prediction from Illumina data is correct and evaluate the Nanopore prediction accordingly.

FN=false negative; R=number of resistant samples; FP=false positive; S=number of susceptible samples; FNR=false negative rate; FPR=false positive rate; PPV=positive predictive value; NPV=negative predictive value; CI=Wilson score confidence interval

Drug	FN(R)	FP(S)	FNR(95% CI)	FPR(95% CI)	PPV(95% CI)	NPV(95% CI)
Amikacin	0(12)	2(138)	0.0% (0.0-24.2%)	1.4% (0.4-5.1%)	85.7% (60.1-96.0%)	100.0% (97.3-100.0%)
Capreomycin	0(12)	2(138)	0.0% (0.0-24.2%)	1.4% (0.4-5.1%)	85.7% (60.1-96.0%)	100.0% (97.3-100.0%)
Ciprofloxacin	1(17)	0(133)	5.9% (1.0-27.0%)	0.0% (0.0-2.8%)	100.0% (80.6-100.0%)	99.3% (95.9-99.9%)
Ethambutol	2(56)	0(94)	3.6% (1.0-12.1%)	0.0% (0.0-3.9%)	100.0% (93.4-100.0%)	97.9% (92.7-99.4%)
Isoniazid	0(80)	7(70)	0.0% (0.0-4.6%)	10.0% (4.9-19.2%)	92.0% (84.3-96.0%)	100.0% (94.3-100.0%)
Kanamycin	0(13)	2(137)	0.0% (0.0-22.8%)	1.5% (0.4-5.2%)	86.7% (62.1-96.3%)	100.0% (97.2-100.0%)
Moxifloxacin	1(17)	0(133)	5.9% (1.0-27.0%)	0.0% (0.0-2.8%)	100.0% (80.6-100.0%)	99.3% (95.9-99.9%)
Ofloxacin	0(17)	0(133)	0.0% (0.0-18.4%)	0.0% (0.0-2.8%)	100.0% (81.6-100.0%)	100.0% (97.2-100.0%)
Pyrazinamide	1(31)	0(119)	3.2% (0.6-16.2%)	0.0% (0.0-3.1%)	100.0% (88.6-100.0%)	99.2% (95.4-99.9%)
Rifampicin	0(78)	0(72)	0.0% (0.0-4.7%)	0.0% (0.0-5.1%)	100.0% (95.3-100.0%)	100.0% (94.9-100.0%)
Streptomycin	2(46)	2(104)	4.3% (1.2-14.5%)	1.9% (0.5-6.7%)	95.7% (85.5-98.8%)	98.1% (93.3-99.5%)

One thing I have changed from the mykrobe papers is the use of FNR/FPR instead of VME/ME. I find VME and ME confusing. Happy to change it back though if this is an issue

mbhall88 commented 3 years ago

Final table that will be used in the paper (pending major issues)

Drug	FN(R)	FP(S)	FNR(95% CI)	FPR(95% CI)	PPV(95% CI)	NPV(95% CI)
Amikacin	0(12)	2(138)	0.0% (0.0-24.2%)	1.4% (0.4-5.1%)	85.7% (60.1-96.0%)	100.0% (97.3-100.0%)
Capreomycin	0(12)	2(138)	0.0% (0.0-24.2%)	1.4% (0.4-5.1%)	85.7% (60.1-96.0%)	100.0% (97.3-100.0%)
Ciprofloxacin	1(17)	0(133)	5.9% (1.0-27.0%)	0.0% (0.0-2.8%)	100.0% (80.6-100.0%)	99.3% (95.9-99.9%)
Ethambutol	3(57)	0(93)	5.3% (1.8-14.4%)	0.0% (0.0-4.0%)	100.0% (93.4-100.0%)	96.9% (91.2-98.9%)
Isoniazid	0(80)	7(70)	0.0% (0.0-4.6%)	10.0% (4.9-19.2%)	92.0% (84.3-96.0%)	100.0% (94.3-100.0%)
Kanamycin	0(13)	2(137)	0.0% (0.0-22.8%)	1.5% (0.4-5.2%)	86.7% (62.1-96.3%)	100.0% (97.2-100.0%)
Moxifloxacin	1(17)	0(133)	5.9% (1.0-27.0%)	0.0% (0.0-2.8%)	100.0% (80.6-100.0%)	99.3% (95.9-99.9%)
Ofloxacin	0(17)	0(133)	0.0% (0.0-18.4%)	0.0% (0.0-2.8%)	100.0% (81.6-100.0%)	100.0% (97.2-100.0%)
Pyrazinamide	1(31)	0(119)	3.2% (0.6-16.2%)	0.0% (0.0-3.1%)	100.0% (88.6-100.0%)	99.2% (95.4-99.9%)
Rifampicin	1(79)	0(71)	1.3% (0.2-6.8%)	0.0% (0.0-5.1%)	100.0% (95.3-100.0%)	98.6% (92.5-99.8%)
Streptomycin	2(46)	2(104)	4.3% (1.2-14.5%)	1.9% (0.5-6.7%)	95.7% (85.5-98.8%)	98.1% (93.3-99.5%)

mbhall88 commented 2 years ago

Based on #80, we will not use mykrobe minor resistance calls. @iqbal-lab should we consider minor resistant calls susceptible or ignore drug predictions for this samples (or drugs)?

iqbal-lab commented 2 years ago

Susceptible

mbhall88 / head_to_head_pipeline