In depth analysis of variant calling on matching FF and FFPE metastatic prostate tumor samples
We aimed at testing to what extent WGS of FF- and FFPE-derived material results in the identification of the same somatic SNVs. Hereto, we used one FF and one FFPE sample of the same metastatic prostate tumor from patient UZ001 (see Materials and Methods). We opted for a metastatic sample as this is in general more homogenous than primary samples, reducing differences in variant calls due to ITH and sampling variation. A matching blood sample, taken from the same patient, was used as a reference to identify germline mutations.
Quality checking showed that in the matching FF and FFPE samples, a comparable number of variants was called (see Suppl. Table 4, Additional File 1). No bias towards specific mutational patterns was observed, indicating the absence of prominent artefacts in the FFPE sample (see Suppl. Table 4, Additional File 1). Both samples showed a similar average coverage (see Suppl. Table 5, Additional File 1). To ensure that comparing the somatic variants called from the matching FF and FFPE sample would be independent of the specificities of the variant caller, we used four somatic variant callers, i.e. Strelka2 (15), Mutect2 (16), VarScan2 (18) and Shimmer (17).
Comparison of variants detected in respectively the FF and FFPE sample
At first, we compared the extent to which each caller identifies the same variants in the matching FF and FFPE sample. All variant callers except VarScan2 were run under default settings (see Materials and Methods). In general, more calls were reported in the FFPE than in the FF sample, except for VarScan2 (see Suppl. Fig. 5, Additional File 2). Table 5 shows the overlap between somatic calls in the matching FF and FFPE samples in terms of sensitivity and precision using the calls from the same variant caller on the FF sample as gold standard. Overall, variant callers have an average sensitivity and precision of respectively 50.97% and 52.41%, meaning that 50.97% of the variants from the FF sample are also detected in the FFPE sample and 52.41% of the FFPE somatic variants are detected in the FF sample. We also report the F1-score, which is the harmonic mean of sensitivity and precision. On average, the variant callers achieved an F1-score of 50.05%. Among the four variant callers considered in this study, Strelka2 achieved the highest F1-score.
Table 5. Performance measures of calls considering the FF sample as gold standard for each variant caller.
|
FF
|
FFPE
|
Overlap
|
Sensitivity
|
Precision
|
F1-score
|
Strelka2
|
6292
|
6761
|
4225
|
0.6715
|
0.6249
|
0.6474
|
Mutect2
|
10460
|
11815
|
5755
|
0.5502
|
0.4871
|
0.5167
|
VarScan2
|
4067
|
1760
|
883
|
0.2171
|
0.5017
|
0.3031
|
Shimmer
|
8109
|
10080
|
4865
|
0.6000
|
0.4826
|
0.5349
|
To get an idea of how good the sensitivity and precision from Table 5 are, we compared the calls from different variant callers on the same sample. In general, the difference between variant callers in the same sample is larger than the difference between the FF and the FFPE sample observed for the same caller (see Suppl. Fig. 5 and 6, Additional File 2). This implies that the choice of variant caller is at least as important as whether or not the sample was FF or FFPE. Overall, the overlap between variant callers on the same sample is low, especially for the FFPE sample (see Suppl. Fig. 6, Additional File 2). Shimmer agrees the least with the other variant callers, whereas Strelka2 and Mutect2 tend to mutually agree in both the FF and the FFPE samples (see Suppl. Table 6 and 7, Additional File 1).
To further confirm that no bias towards specific mutational patterns were present, for each caller, we measured the cosine similarity between the mutational profiles between the FF and the FFPE sample. These profiles represent the combination of the 96 different mutation types in each sample. Suppl. Table 8 (Additional File 1) illustrates that the mutational profiles of variants reported by each caller were consistent between the FF and the FFPE sample except for Shimmer.
While the choice of variant caller turns out to play a pivotal role in the calls that are obtained, the basic analysis above does not account for the properties of the called variants. Indeed, we will demonstrate that the variants considered in Table 5 represent both true variants and artefacts of the used variant caller, consistently made in the FF and the FFPE sample. To assess the relevance of the calls made by each of the callers in the overlap between the FF and the FFPE sample, we first assessed to what extent each of the callers tends to reconstruct the same (sub)clonal subpopulations in both the FF and the FFPE samples and secondly whether high confidence calls are consistent between both sample types.
Assessing whether variants called in the FF and FFPE samples represent the same (sub)clonal populations
Given that the FF and the FFPE sample should have a similar subclonal structure, as they derive from the same metastatic tumor, we judged the relevance of the variants called by either method by assessing whether they corresponded in the FF and the FFPE to the same (sub)clonal populations. To visualize the different populations in each sample, we plotted for all detected variants the coverage as a function of the VAF and the corresponding histogram representing the distribution of the VAFs (26). Fig. 1 shows the results for Strelka2 in both the FF and the FFPE samples. Similar figures for the three other variant callers can be found in the Supplementary Figures, Additional File 2. Note that the VAF estimate of the major distributions slightly differs between the different tools because of small discrepancies in their definition of the VAF (see Suppl. Fig. 7, Additional File 2).
According to Strelka2 (Fig. 1), two distinct clusters of variants can be observed in the FF sample, denoted by the blue bars in the histogram: a small cluster at a VAF of 0.05 and a larger cluster around a VAF of 0.25. In the FFPE sample, only one cluster can be observed around a VAF of 0.15. The calls that are common to the FF and the FFPE sample are shown in orange, demonstrating that the peak at a VAF of 0.25 in FF does indeed shift to 0.15 in FFPE. Almost all variants from the largest cluster in the FF sample, i.e. at a VAF of 0.25, are identified in the FFPE sample but most of the FF calls belonging to the smaller peak at 0.05 are not detectable in the FFPE sample. This shows that most variants from the major population identified in the FF sample are recovered in the FFPE sample.
Figure 1. (Sub)clonal populations detection in the FF (left) and the FFPE (right) sample using Strelka2. The upper panel shows coverage as a function of the VAF (26), where a higher variance in the coverage can be observed for FFPE. The lower panel shows the distribution of the VAFs. The blue distribution denotes all calls made in a given sample, while the orange distribution shows only the calls common to FF and FFPE.
We observed similar results for Mutect2 with the peak representing the major cluster at VAF 0.25 in the FF sample being shifted to a lower VAF in the FFPE sample (see Suppl. Fig. 8, Additional File 2). For VarScan2, only a small proportion of the variants detected at an average VAF of 0.25 in the FF is recovered in the FFPE sample (see Suppl. Fig. 9, Additional File 2). Shimmer could barely detect the cluster of variants at VAF 0.25 in the FF and completely misses the major cluster in the FFPE sample (see Suppl. Fig. 10, Additional File 2). A deeper analysis shows that Shimmer reports many somatic variants in the tumor sample that have a VAF above zero in the normal sample (putative germline calls). This behavior is not expected and not observed for any of the other variant callers, explaining the poor overlap between Shimmer and the other variant callers (see Suppl. Tables 6 and 7, Additional File 1). Only keeping for Shimmer the variants with zero VAF in the normal sample filtered those putative germline calls and allows to better recover a cluster of variants at VAF 0.25 in the FF, but still not in the FFPE sample (see Suppl. Fig. 11, Additional File 2).
In addition to the cluster at VAF 0.25, VarScan2 also reports a peak of variants at an average VAF of 0.5 in the FF sample. These are likely germline variants as more than 65% of variants with VAF above 0.5 are present in dbSNP database (see Suppl. Table 9, Additional File 1). Except for Mutect2, which has a similar number of germline calls above and below a VAF of 0.5, other callers reported up to 8 times more calls from dbSNP at a VAF above 0.5. This suggests that calls reported with VAF above 0.5 are more likely to be false positives. Indeed, most of the calls with VAF above 0.5 detected consistently in both the FF and the FFPE samples tend to be residual germline calls.
Hence, the cluster of variants detected at VAF 0.25 likely represents the major population of variants in the metastatic sample. All variant callers except Shimmer could at least partially recover these variants in the FFPE sample albeit at lower VAF (see Suppl. Fig. 12, Additional File 2). Imposing a zero VAF criterion in the normal sample helped Shimmer to recover the major cluster of variants that was also recovered by the other callers, but only in the FF sample. In addition, only a minority of calls (813) are retained of which most (604) were also reported by other callers. This indicates that many of the calls uniquely made by Shimmer without imposing the zero VAF criteria and reported in Table 5 are likely spurious calls, despite being consistently detected in both the FF and the FFPE samples. For the remainder of the analysis, variants reported by Shimmer with a positive VAF in the normal sample were filtered.
Comparing for each of the different callers the significance levels between the FF and FFPE sample.
To assess the relevance of the calls in the overlap between the FF and the FFPE sample, we also compared, for each caller, the distributions of the significance scores of the somatic variant calls detected in the FFPE sample that were also called in the FF sample, hereby assuming that the FF sample is less prone to artefacts and constitutes the reference as to what should be detected. Table 6 shows that calls reported in both samples typically received lower significant scores in the FFPE than in the FF sample, which can be explained by the discrepancies in VAF observed in Fig. 1.
Table 6. Average significance score for somatic variants reported in both samples for each variant caller. For Strelka2 and Mutect2, a higher Somatic EVS and TLOD means a higher confidence in the calls, while for VarScan2 and Shimmer a lower value implies a higher confidence.
|
Strelka2
|
Mutect2
|
VarScan2
|
Shimmer
|
FF
|
17.33
|
51.88
|
0.0024
|
0.0160
|
FFPE
|
14.34
|
26.00
|
0.0038
|
0.0166
|
In addition, the boxplots in Fig. 2 show that the FFPE calls that were also made in the FF sample, received a relatively higher significance score than FFPE calls not made in the FF sample. This shows that the most reliable variants in the FFPE sample generally correspond to those detected in the FF sample.
Figure 2. Boxplots comparing the significance level of calls from FFPE reported or not in FF sample. For Strelka2 and Mutect2, a higher Somatic EVS and TLOD means a higher confidence in the calls, while for Varscan2 and Shimmer a lower value implies a higher confidence.
For each caller, we calculated the correlation between the significance levels reported in the FF and in the FFPE sample of somatic variants called in both samples. Suppl. Table 10 (Additional File 1) shows that the significance levels of Strelka2 and Mutect2 are significantly correlated in both the FF and the FFPE samples. In addition, the significance levels between the callers themselves seems to be consistent between the FF and the FFPE sample. For Shimmer and VarScan2, this consistency between both samples cannot be observed. Furthermore, all callers assign a similar relative rank to common variants in the FF sample but not necessarily in the FFPE sample (i.e. the same variant is ranked high by all callers in the FF sample but in the FFPE sample the ranks for that variant vary depending on the caller). This indicates that these variant callers are less performant in the FFPE sample. A possible explanation for this could be that the underlying hypergeometric testing procedure (used in both VarScan2 and Shimmer) cannot cope with the lower VAFs present in the FFPE sample.
Relation between the (sub)clonal structure and the significance level
To investigate the relation between the (sub)clonal structure and the significance of the calls, we map for each variant caller the 25% highest confidence calls on the coverage versus VAF plots. The upper panel of Fig. 3 shows how, for Strelka2, these most significant calls are located around a VAF of 0.25 in the FF and 0.15 in the FFPE sample, and hence make up the aforementioned major subpopulation that was detected in both the FF and the FFPE samples.
Figure 3. Coverage versus VAF for somatic variants reported by Strelka2, comparing FF (left) against FFPE (right). This plot is identical to the upper panel of Fig. 1 but with a color used to indicate the most significant calls. The upper panel shows the 25% highest confidence calls in orange and the lower confidence in blue. The lower panel shows which calls are also found by other callers where blue= unique calls, orange = calls reported by 2 callers, green = calls reported by 3 callers, red = calls reported by 4 callers.
For Mutect2, most significant calls also belonged to this cluster of variants (see Suppl. Fig. 13, Additional File 2). For VarScan2 and Shimmer (see Suppl. Fig. 14-16, Additional File 2), a similar effect is observed, although not as pronounced as in Strelka2 and Mutect2. Subsequently, many of the highly significant calls belong to the major subpopulation at a VAF of 0.25 in the FF and 0.15 in the FFPE sample; such that the analysis of the significance levels and the clonal subpopulations points at the existence of a highly confident subset of variants, present in both the FF and the FFPE samples.
In addition, the lower panel in Fig. 3 shows that many variants called by other callers are also located in the variant cluster (see also the lower panels of Suppl. Fig. 13-16, Additional File 2). This is in line with the observation that for each variant caller (except Shimmer), calls made by any of the other three variant callers in general received a higher significance especially in the FF sample (Suppl. Fig. 17 and 18, Additional File 2). Using two criteria, based the (sub)clonal structure (Fig. 1) and significance score of the variants (Fig. 2), we could see that calls common to the FF and the FFPE sample tend to belong to the major subpopulation and are highly significant. Importantly, the lower panel of Fig. 3 and Suppl. Fig. 17 and 18 (Additional File 2) show that variants called by more than one caller tend to satisfy these two criteria as well. In the next section, we investigated how well these calls, reported by more than one variant caller, overlap between the FF and the FFPE sample.
Optimizing for each variant caller its significance threshold by optimizing the overlap in variants between the FF and FFPE sample.
To assure that only biologically relevant calls are considered, we first identified a highly reliable subset of calls in the FF sample. This subset, referred to as the ground truth, can then be used to more qualitatively assess how well each caller can recover these highly reliable calls in the FFPE sample. From our analysis above, an ideal ground truth would consist of all highly reliable calls made in the reference FF sample that also represent the major population in the metastatic sample. We have shown that calls made by at least two callers tend to satisfy these criteria. Therefore, the ground truth is defined as the union of all calls that were detected by at least two callers in the FF sample.
Our previous analysis also shows that the subpopulation represented by the ground truth is also present in the FFPE sample albeit at lower VAF (median VAF of 0.1083 in FFPE instead of 0.2323 in FF). Hence, fully recovering the subpopulation of highly significant variants from the FFPE sample will require the identification of a significance threshold in the FFPE sample. Because calls in the FF sample were reported with a higher significance scores, the threshold in the FFPE sample will typically be lower than the threshold that would be necessary in the FF sample to capture the ground truth (see Suppl. Fig. 17 and 18, Additional File 2). In addition, previous analysis shows that it is feasible to set a threshold in the FFPE sample as calls common to the FF and the FFPE sample tend to have a higher significance level in the FFPE sample (see Fig. 2). However, the lower VAF in the FFPE sample complicates identifying a threshold that distinguishes the true calls from the noise (the significance distribution of noisy and true calls starts overlapping, see Suppl. Fig. 18, Additional File 2). We determined, for each variant caller, a threshold in the FFPE sample that optimized the F1-score with the FF-derived ground truth (i.e. that optimizes the tradeoff between precision and sensitivity in recovering the ground truth). Table 3 shows for each caller their sensitivity and precision in recovering the ground truth. This table also shows that Strelka2 and Mutect2 are performing the best in recovering from the FFPE sample the calls that belong to the ground truth. Table 7 now quantitatively illustrates how Shimmer underperforms on the ground truth despite calling consistently the same mutations in the FF and the FFPE sample, which was shown in Table 5. This suggests that most of the calls made by Shimmer in the overlap between the FF and the FFPE sample are likely spurious. The same is to some extent true for VarScan2 because of the high number of residual germline calls.
Table 7. Optimized F1-scores of calls made by each variant caller considering FF sample as gold standard. Optimized F1-scores of calls made by each variant caller considering FF sample as gold standard.
Caller - threshold
|
FF (gold std.)
|
FFPE
|
Overlap
|
Sensitivity
|
Precision
|
F1-score
|
Strelka2 – 9.5
|
4656
|
4559
|
3316
|
0.7122
|
0.7274
|
0.7197
|
Mutect2 – 13
|
4656
|
5658
|
3418
|
0.7341
|
0.6041
|
0.6628
|
VarScan2 – 0.00995
|
4656
|
1755
|
425
|
0.1045
|
0.2422
|
0.1460
|
Shimmer – 0.0495
|
4656
|
262
|
16
|
0.0197
|
0.0611
|
0.0298
|
Table 7 gives an estimate of the expected overlap between the FF and the FFPE sample after optimizing the thresholds in the FFPE sample using the ground truth based on the variants detected in the FF sample. However, often only FFPE samples are available, such that the stringency thresholds cannot be optimized based on observations in the FF sample. Therefore, rather than optimizing the threshold for each caller separately, we assessed to what extent combining the output of different variant callers in the FPPE sample allows recovering the ground truth from the FF sample.
Table 8. Strategies to retrieve the ground truth calls from the FF in the FFPE sample.
Reported by… (in FFPE)
|
FF (gold std.)
|
FFPE
|
Overlap
|
Sensitivity
|
Precision
|
F1-score
|
at least 1 caller
|
4656
|
16020
|
4155
|
0.8924
|
0.2594
|
0.4019
|
at least 2 callers
|
4656
|
4232
|
3684
|
0.7912
|
0.8705
|
0.8290
|
at least 3 callers
|
4656
|
340
|
325
|
0.0698
|
0.9559
|
0.1301
|
all 4 callers
|
4656
|
8
|
8
|
0.0017
|
1
|
0.0034
|
Table 8 shows how taking the intersection of the four variant callers maximizes the precision but comes at the expense of losing almost all sensitivity. As discussed above, Shimmer can barely retrieve the cluster representing the major subpopulation in the FFPE sample (see Suppl. Fig. 15 and 16, Additional File 2) and most of the calls retrieved by Shimmer in the FFPE sample are unique (even after correcting for the so-called somatic calls with high VAF in the normal sample). Because the intersection seems too strict and limits the sensitivity, we considered calls reported by at least three callers. It results in a high precision but still relatively low sensitivity in recovering the ground truth (Table 8) because VarScan2 reports less calls in the FFPE than in the FF sample and loses many true positive calls. Using the calls returned by at least two of the four callers in the FFPE sample drastically increases sensitivity while only slightly decreasing precision. The performance here is even better than the performance of the best caller (Strelka2) that was obtained after optimizing its stringency thresholds based on the ground truth. This indicates that there is some complementarity in the calls made by different callers.
Note that when considering only variants called by at least two variant callers in both samples, the cosine similarity between the mutational profiles increased (0.995). Moreover, we decomposed the mutational profiles of the variants called in both the FF and the FFPE samples into COSMIC signatures (20). Some of those signatures are representative for sequencing artefacts (see Suppl. Table 1, Additional File 1). When considering the union of calls made by the four callers in the FF and the FFPE sample, the contribution of the COSMIC signatures representative for sequencing artefacts to the total mutational profile was of 15.21% and 17.68% respectively. When considering only variants called by at least two callers those proportions decreased significantly (8.22% for the FF and 9.36% for the FFPE sample) indicating that considerably less spurious calls related to sequencing artefacts were made.
Validation of the ‘at least two variant callers strategy’ on samples from the 100,000 Genomes Project
Comparison of variants detected in respectively the FF and FFPE sample
Nine paired FF-FFPE samples from the same tumors obtained from the 100,000 Genomes Project England have been analyzed in order to validate our approach. As for patient UZ001, variant callers reported systematically more calls in the FFPE than in the FF samples (except for Shimmer in GeL365 patient). Either Strelka2 or Mutect2 achieved the highest F1-scores. Overall, the highest F1-score per variant caller varied between 10.12% (GeL032) to 77.30% (GeL300) with an average of 44.10% (see Table 9). On average, the variant callers were less consistent in the FFPE samples than in the FF samples and the difference between callers in the FFPE sample was larger than the difference between the FF and the FFPE samples using the same variant caller. This is consistent with our findings on the metastatic prostate sample. The difference between the variant callers in the FF sample was larger than the difference of each caller between the FF and the FFPE samples for GeL004, GeL300, GeL365 and UZ001 (see Suppl. Table 11, Additional File 1). In addition, for two tumor samples, the VAF was lower in the FFPE sample than in the FF sample (GeL007, GeL028), for one of them the opposite holds true (GeL008). For the rest of the tumor samples, the VAF was consistent between the FF and the FFPE sample (see Suppl. Table 12, Additional File 1). In general, Mutect2 and Strelka2 reported the most similar mutational profiles between the FF and the FFPE samples (see Suppl. Table 13, Additional File 1). No bias towards specific mutational patterns was observed in the FFPE samples.
Table 9. Summary table of the consistency (F1-score) between samples.
Patient ID
|
Strelka2
|
Mutect2
|
VarScan2
|
Shimmer
|
Maximum
|
At least 2
|
Improvement
|
UZ001
|
0.6474
|
0.5167
|
0.3031
|
0.5349
|
0.6474
|
0.8290
|
0.1816
|
GeL007
|
0.1533
|
0.1563
|
0.0297
|
0.0061
|
0.1563
|
0.2317
|
0.0754
|
GeL008
|
0.2344
|
0.0823
|
0.0168
|
0.0023
|
0.2344
|
0.3895
|
0.1551
|
GeL024
|
0.2298
|
0.1756
|
0.0210
|
0.0107
|
0.2298
|
0.4203
|
0.1904
|
GeL028
|
0.3586
|
0.2485
|
0.0449
|
0.0335
|
0.3586
|
0.5326
|
0.1740
|
GeL004
|
0.5656
|
0.5226
|
0.1502
|
0.2200
|
0.5656
|
0.7149
|
0.1493
|
GeL032
|
0.0646
|
0.1012
|
0.0034
|
0.0004
|
0.1012
|
0.0647
|
-0.0365
|
GeL065
|
0.6308
|
0.5309
|
0.1604
|
0.0304
|
0.6308
|
0.6345
|
0.0037
|
GeL300
|
0.7724
|
0.7730
|
0.2358
|
0.5837
|
0.7730
|
0.8716
|
0.0986
|
GeL365
|
0.6532
|
0.6155
|
0.2473
|
0.1530
|
0.6532
|
0.8051
|
0.1519
|
Average
|
0.4310
|
0.3753
|
0.1212
|
0.1575
|
0.4350
|
0.5494
|
0.1144
|
Validation of ‘the at least two strategy’ to call variants
To validate our ‘at least two’ variant calling approach, we ran all four variant callers on the matched FF and FFPE samples for each patient. Then, we combined the variants using the at least two strategy that gave the best performance for the UZ001 patient. We calculated the F1-score between the variants called in the FF and the FFPE samples by each caller and compared them with the F1-score obtained when considering only variants reported by at least two callers in the FF and the FFPE sample. The results are shown in Table 9. For most samples, it can be appreciated that our variant calling approach significantly improves upon the best single caller score under default conditions. The only exception being patient GeL032, which receives a very low F1-score for any caller. Overall, the F1-score improves on average by 0.1144, suggesting that our variant calling strategy allows detecting robustly somatic SNVs in the FFPE samples, offering both a high sensitivity and precision. More detailed information about the sensitivity and precision of all variant callers on each sample are available in Suppl. Table 14 (Additional File 1). Note that selecting variants called by at least two callers also increases the correspondence of the mutational profiles between the FF and the FFPE samples (higher cosine similarity scores, see Suppl. Table 15, Additional File 1).
Similarly to the analysis we conducted for UZ001 patient, the mutational profiles of the 100,000 Genomes Project England samples have been decomposed into COSMIC signatures thereby assessing the contribution of artefact signatures. In general, the contribution of artefact signatures to the FF and the FFPE mutational profiles was of the same order of magnitude for each of the paired FF-FFPE samples. All variant callers except VarScan2 reported a slightly higher proportion of artefact signatures for the FFPE sample. Interestingly, Table 10 shows that when considering only variants called by at least two variant callers, the proportion of signatures representative of artefacts (see Suppl. Table 1, Additional File 1) is significantly lower. This further confirms that our strategy successfully removes sequencing artefacts.
Table 10. Artefact contribution to the mutational profiles of the FF and the FFPE sample per caller.
Patient ID
|
Strelka2
|
Mutect2
|
VarScan2
|
Shimmer
|
At least 2
|
FF
|
FFPE
|
FF
|
FFPE
|
FF
|
FFPE
|
FF
|
FFPE
|
FF
|
FFPE
|
UZ001
|
0.1419
|
0.0847
|
0.1024
|
0.0757
|
0.2010
|
0.0981
|
0.3252
|
0.6759
|
0.0822
|
0.0936
|
GeL007
|
0.1626
|
0.5381
|
0.1285
|
0.1824
|
0.3843
|
0.1981
|
0.1555
|
0.1863
|
0.0714
|
0.1132
|
GeL008
|
0.2133
|
0.1377
|
0.1302
|
0.2501
|
0.4345
|
0.1925
|
0.2466
|
0.1631
|
0.0796
|
0.1704
|
GeL024
|
0.4953
|
0.3111
|
0.1244
|
0.1261
|
0.5331
|
0.2191
|
0.1882
|
0.0556
|
0.1050
|
0.0899
|
GeL028
|
0.1344
|
0.1670
|
0.1041
|
0.0938
|
0.4947
|
0.1503
|
0.1152
|
0.0509
|
0.0666
|
0.0859
|
GeL004
|
0.2067
|
0.2896
|
0.1099
|
0.1611
|
0.3565
|
0.0932
|
0.1940
|
0.3090
|
0.0590
|
0.0775
|
GeL032
|
0.1164
|
0.1168
|
0.0833
|
0.2664
|
0.4955
|
0.1528
|
0.1474
|
0.3310
|
0.0747
|
0.0506
|
GeL065
|
0.1047
|
0.1686
|
0.0580
|
0.1127
|
0.1180
|
0.1425
|
0.1662
|
0.1458
|
0.1574
|
0.0816
|
GeL300
|
0.0822
|
0.2095
|
0.0781
|
0.0827
|
0.0583
|
0.0758
|
0.1022
|
0.1120
|
0.0544
|
0.0558
|
GeL365
|
0.0886
|
0.1701
|
0.0662
|
0.1192
|
0.0862
|
0.1051
|
0.0894
|
0.1362
|
0.0331
|
0.0436
|
Mean
|
0.1746
|
0.2200
|
0.0985
|
0.1201
|
0.3162
|
0.1402
|
0.1730
|
0.2030
|
0.0783
|
0.0817
|
Assessing the extent to which intra-tumor heterogeneity (ITH) explains discrepancies in variants called between the FF and FFPE samples
Table 10 shows that the ‘at least two’ variant calling strategy allows to reliably detect variants in most samples. By applying this strategy on both the FF and the FFPE samples, we can investigate why certain variants are unique to each sample type. Based on Fig. 1, the most plausible explanation would be ITH that results in sampling different subclonal variants. If ITH plays a key role in the difference between the FF and the FFPE samples, then we would expect that the subclonal variants show a lower overlap than the clonal ones. However, to properly define clonal and subclonal variants, it is necessary to correct the obtained VAFs for copy number (CN) status and sample purity.
To perform CN correction, the CN variant calling pipeline from GATK was used to identify regions that are amplified (+), deleted (-) or neutral (0) (see Materials and Methods). Although there is a limited number of publications demonstrating the feasibility of CN calling on FFPE, our results clearly showed a distinction between CN alterations called in the FF and the FFPE samples (see Suppl. Table 2, Additional File 1). The most plausible explanation for this observation is the presence of short CN altered regions that are only observed in the FFPE samples (see Suppl. Fig. 3, Additional File 2). It is highly improbable that these regions are effectively CN altered and therefore variants lying in these regions would be falsely corrected based on their CN status. To avoid introducing such additional bias, only variants called by at least two variants callers mapping to segments diploid in both the FF and the FFPE samples were kept (bona fide variants, see Materials and Methods and Suppl. Table 3, Additional File 1). Additionally, we also discarded variants located on the mitochondrial DNA and the sex chromosomes.
To estimate sample purity, we compared the pathologist purity estimates to those of two tools, TPES (21) and FACETS (22). In general, the purity estimates for the same sample showed little consistency (see Suppl. Table 16, Additional File 1). Therefore, correcting for purity risks introducing yet another bias. Therefore, we aimed at identifying clonal variants for each sample individually, using a Kernel Density Estimation (KDE) approach (see Materials and Methods). For each sample, only considering bona fide diploid variants, the clonal variants are identified based on the VAF distribution. The KDE approach then allows finding the clonal population, regardless of the sample purity. Fig. 4 illustrates the KDE approach on samples from UZ001 patient and demonstrates that the agreement between the FF and the FFPE sample is higher when considering only clonal variants.
Figure 4. Evaluating the overlap between clonal variants from FF and FFPE for UZ001. By first defining the clonal and subclonal variants (see text), it is possible to calculate the overlap between clonal variants only. Compared to Fig.1, only variants using the at least two approach are shown, which are known to lie in diploid regions in both FF and FFPE and can be considered clonal. Overlap: red (green) refers to variants in the FFPE (FF) sample also found in the FF (FFPE) sample. Clearly, these additional filtering steps lead to an appreciable improvement in overlap between the FF and the FFPE sample.
Table 11 shows the overlap between clonal variants in the FF and the FFPE samples where variants were obtained using the ‘at least two’ variant calling strategy. The results can be compared to subclonal variants (defined as all bona fide diploid variants that are not clonal). Note that, as expected, clonal variants are in general easier to retrieve than their subclonal counterparts are. There is one example (GeL007) were the opposite holds, but further investigation showed that many clonal variants from the FF sample appeared as subclonal variants in the FFPE sample (see Suppl. Fig. 19, Additional File 2). These results demonstrate that some of the discrepancies between the FF and the FFPE samples from Table 10 can be attributed to ITH. Note that for some patients (e.g. patient UZ001, see Fig. 4), we completely lose the clonal structure in the FFPE sample, such that no subclonal variants can be defined in this sample. For each patient, we also explicitly verified how many subclonal variants in the FFPE sample were classified as clonal in the FF sample and vice versa (see Suppl. Table 17, Additional File 1). For most patients, the clonal structure in the FFPE sample agrees well with the FF sample, the most notable exception being samples from GeL032 patient.
Table 11. Number of clonal and subclonal variants detected in diploid regions common to FF and FFPE. For both clonal and subclonal variants the F1-score between the FF and the FFPE variants was calculated. Variants were called using the at least two variant calling strategy.
Clonality
|
Clonal
|
Subclonal
|
At least 2
|
Patient ID
|
FF
|
FFPE
|
Overlap
|
F1-score
|
FF
|
FFPE
|
Overlap
|
F1-score
|
F1-score
|
UZ001
|
3229
|
3246
|
2795
|
0.86
|
372
|
0
|
0
|
0.00
|
0.83
|
GeL007
|
662
|
25
|
1
|
0.04
|
1862
|
1183
|
223
|
0.19
|
0.23
|
GeL08
|
1109
|
2417
|
996
|
0.41
|
425
|
0
|
0
|
0.00
|
0.39
|
GeL024
|
1016
|
296
|
244
|
0.82
|
275
|
1413
|
4
|
0.00
|
0.42
|
GeL028
|
1037
|
1585
|
908
|
0.57
|
721
|
0
|
0
|
0.00
|
0.53
|
GeL004
|
1074
|
1074
|
965
|
0.90
|
710
|
630
|
297
|
0.47
|
0.71
|
GeL032
|
3249
|
202
|
126
|
0.62
|
506
|
2103
|
4
|
0.00
|
0.06
|
GeL065
|
1434
|
1362
|
1287
|
0.94
|
587
|
552
|
321
|
0.58
|
0.63
|
GeL300
|
9660
|
9821
|
8890
|
0.91
|
354
|
678
|
16
|
0.02
|
0.87
|
GeL365
|
1099
|
1013
|
962
|
0.95
|
1331
|
1309
|
937
|
0.72
|
0.80
|