Analysis of held-out cancer type classification results

[jjc2718]

In my previous PR (#20), I did an analysis of mutation prediction results from stratified cross-validation experiments; i.e. the training set was composed of the same cancer types in the same proportions as the test set. In this PR, I'm adding a similar analysis, but I'm now holding out single cancer types. I did two types of experiments:

1. Train on 75% of data from a single cancer type, test on the other 25% of data from that cancer type
2. Train on the same 75% of data + data from all other cancer types in TCGA with sufficiently many mutations in the given gene, and test on the same 25% of data as before

The plots show p-values from t-tests comparing cross-validation results from each of these setups and the negative control (shuffled labels), and between the two setups (more documentation in the analysis notebooks).

I'm planning to do some follow-up on these results in future PRs.

[ajlee21]

The code all looks good to me! I like all the the markdown comments where you are explaining your logic and what you expect to see.

Most of my comments are about points I got a little confused. I think I'm missing more of the big picture behind what you're doing. I'm sure once you answer them the plots will make sense to me:)

[ajlee21]

Some comments for this notebook overall, couldn't put in line comments since it appears these are not new changes so I apologize that this might be a little out of scope and you can address them in a separate PR if you'd like.

1. I tend to find it helpful to have a description of the experiment I am performing at the top of the notebook just to orient myself.

2. Just wanted to clarify the term stratified here. So you're saying that you're training set includes say 10 samples of cancer type A, 10 of cancer type B, 10 of cancer type C (total of 30 samples each with 1/3 protion). And so your test set contains maybe 9 total samples with 3 samples of cancer type A, 3 of cancer type B, 3 of cancer type C. So the proportions are the same in the test and training..?

3. Trying to understand your dfs. So for the first 3 rows top50_df you have auroc and aupr that tells you how well gene info from training set (including mutation burden, etc) predict mutation status of TP53 (binary i assume) on training/test/validation sets where the labels used to train the model were shuffled. I assume that means you have the same training dataset with multiple labels for mutation status of gene X, Y, Z. So you'd train your model on its ability to predict mutation of gene X, then train your model on its ability to predict mutation of gene Y,...So you have multiple models here?

[jjc2718]

These are all good questions! See answers below:

1. I tend to find it helpful to have a description of the experiment I am performing at the top of the notebook just to orient myself.

Good idea! I'll add this to the top here.

2. Just wanted to clarify the term stratified here. So you're saying that you're training set includes say 10 samples of cancer type A, 10 of cancer type B, 10 of cancer type C (total of 30 samples each with 1/3 protion). And so your test set contains maybe 9 total samples with 3 samples of cancer type A, 3 of cancer type B, 3 of cancer type C. So the proportions are the same in the test and training..?

Yep, exactly (not always exactly the same proportion between train/test sets, but +/- 1 sample I think - this is implemented in StratifiedKFold from scikit-learn).

3. Trying to understand your dfs. So for the first 3 rows top50_df you have auroc and aupr that tells you how well gene info from training set (including mutation burden, etc) predict mutation status of TP53 (binary i assume) on training/test/validation sets where the labels used to train the model were shuffled. I assume that means you have the same training dataset with multiple labels for mutation status of gene X, Y, Z. So you'd train your model on its ability to predict mutation of gene X, then train your model on its ability to predict mutation of gene Y,...So you have multiple models here?

Right - so each row in top50_df and vogelstein_df is one model, trained on the (binary, mutated or not mutated) mutation status of one gene on one cross-validation fold, either on the true labels or the shuffled labels.

We have mutation information for (almost) every gene in the genome - "top_50" and "vogelstein" are two different ways to select the genes to train models on. If we just train models on every gene in the genome, our statistical power to detect true relationships between mutation and gene expression won't be very good (and also it will take forever), so we want to start with sets of known cancer genes to improve power. In each case, we train one model for each gene/true-shuffled/cross-validation fold combination.

Let me know if that doesn't answer your question.

[ajlee21]

This makes sense. Thank you!

As a followup. So you're including those cancer genes instead of all genes. Do you expect most of those genes to have a mutated status = 1 (I guess this'll probably depend on the cancer type)? Would it make sense to include genes that are not mutated in cancers as a control?

[jjc2718]

Do you expect most of those genes to have a mutated status = 1 (I guess this'll probably depend on the cancer type)?

Yeah, we definitely expect many of the cancer-related genes to be frequently mutated in at least a few cancer types (and we are filtering out gene/cancer type combos where less than 5% of samples are mutated, like in my last exploratory data analysis notebook).

Would it make sense to include genes that are not mutated in cancers as a control?

Yes, this is definitely what we're trying to do, but it's hard to choose control genes well. Ideally we'd include some genes that aren't drivers, but this isn't really documented anywhere (and absence of evidence for a gene being a driver in some cancer type isn't the same as evidence of absence; some drivers are just rarely mutated or haven't been studied in depth).

In the past we've used TTN as an example of a gene that isn't thought to be a driver of any cancer type, but remember that lots of genes are mutated in cancer, even those that aren't actually driving the cancer to form. TTN is a large gene that is frequently mutated as a passenger (just by chance) in many cancers, so its mutation status correlates with mutation burden (and thus with cancer type). So (we think) it turns out that TTN mutation status can actually be predicted to some degree from gene expression, because gene expression -> cancer type -> mutation burden -> probability of TTN being mutated.

I guess we could pick smaller genes to reduce the chances of passenger mutations, but I'd have to think about what the best way is to do this.

[ajlee21]

So is the takeaway from these plots that those genes that were found to be DE in cancer vs normal (y-axis), were also found to be most mutated (x-axis), which we'd expect because these mutations will likely change the expression of these genes?

[jjc2718]

It might help if I explain the results_dfs (with the statistical testing results), then maybe that will make it a bit clearer what the plots are showing.

The question I wanted to answer in this notebook is: for which genes can we train a model to predict mutation from gene expression that outperforms the negative control? So for each gene, we ran 8 total cross-validation replicates (4 folds x 2 random seeds), for the experimental case and the control (shuffled) case.

That then gives us 2 distributions of results (in this case we're using AUPR), and we can compare these using a t-test. The plot is a bit like a volcano plot from a DE analysis: on the x-axis it shows the AUPR difference between the true labels and the shuffled labels (positive = better model performance for true labels), and on the y-axis it shows the p-value for the t-test comparing the two distributions. So points in the upper right (better performance for true labels, and low p-values) are the genes we're interested in, showing that we can build effective classifiers on this dataset for these genes.

Does that make more sense? It's not quite a standard use of a volcano plot, but you can interpret it in a similar way.

[jjc2718]

Oh, and as far as the takeaway - in this notebook, there were two for me:

• For the top 50 analysis, we mostly reproduced the results from BioBombe which also used this gene set (some of the less significant hits weren't found in BioBombe, but we should have better statistical power here so it makes sense that we see more results)
• For the Vogelstein analysis, it was surprising/interesting that we saw lots more significant hits than we did for the top 50 analysis! On some level it's not shocking (if a gene is mutated frequently that doesn't necessarily make it a driver, and conversely drivers aren't always frequently mutated across all samples) but seeing visual confirmation of this was neat.

(I'll add this interpretation to the notebook - I wrote this out in some slides I discussed with Casey, but forgot to add it here)

[ajlee21]

I don't think I'm completely following this part. So for a single-cancer type analysis, you're training on genome data from say cancer type A and you're trying to predict mutation status for BRCA in cancer type A? Or predict in a different cancer type B?

For pancancer is it that you train on genomic data from all cancer types except cancer A and predict mutation status of BRCA on cancer type A?

[jjc2718]

So for a single-cancer type analysis, you're training on genome data from say cancer type A and you're trying to predict mutation status for BRCA in cancer type A? Or predict in a different cancer type B?

We're always training and testing on the same cancer type (so train on cancer type A and test on cancer type A).

For pancancer is it that you train on genomic data from all cancer types except cancer A and predict mutation status of BRCA on cancer type A?

In the "pancancer" experiments, we're training on cancer type A + all other cancers, and testing on cancer type A. We haven't tried training on all other cancers (without A) and testing on A; that might be a good next step! That experiment could help us better understand where the signal is coming from (i.e. are the pancancer models relying heavily on the relevant cancer type, or are they actually benefiting from the additional data).

[ajlee21]

I am a little lost on what these plots are telling you, but I think that is because I got a little lost on what these different df are from above.

[ajlee21]

Looking back now, these plots make a lot more sense now. For your last plot where you're identifying cancer genes that were better predicted using your pan-cancer model vs your single cancer model, were you able to identify the certain effects (e.g. in relatively rare cancer types) you mentioned? I guess I'm not sure what you mean by these effects? Do you mean genes that are rarely mutated except in a single cancer type that has a low sample count?

[jjc2718]

were you able to identify the certain effects (e.g. in relatively rare cancer types) you mentioned?

Ehh, kinda. I was hoping to identify more examples that clearly benefit from adding the pan-cancer data (in the top right of the last plot), but there aren't too many unfortunately.

I have some ideas for how to interpret the last result, but I haven't talked with Casey about this plot yet, so he might have ideas on what it means or how to follow up.

I guess I'm not sure what you mean by these effects? Do you mean genes that are rarely mutated except in a single cancer type that has a low sample count?

I was anticipating that for some genes, there would be a strong signal in one cancer type, but this signal would get washed out when the pancancer data was added (i.e. the gene has strong cancer-specific effects that aren't preserved in other cancer types). These genes would be in the top left. We don't really see much of that, though, which in retrospect isn't that surprising since we have a model covariate for cancer type (so the model can just choose to put lots of weight on samples from the relevant cancer type). I should definitely check if this is actually the case, which I think would be a neat follow-up (#23).

I was also anticipating that for some genes, there wouldn't be a strong signal in one cancer type, but the signal would be preserved across cancer types so adding more data would lead to a better model. These genes would be in the top right. This somewhat seems to be the case for a few genes (in particular, NF1 and FBXW7 seem to show up frequently) but I need to think more about what this means and why these genes.