Training deepvariant on RNA TUMOR ONLY SOMATIC variant calling

[chstubbe]

Hi!

I am trying to train deepvariant on RNA SEQ TUMOR ONLY SOMATIC variant calling.
I was already in contact with @danielecook regarding some guidelines on this, which is greatly appreciated. I am now posting on Github to make it easier to troubleshoot. :)

I have prepared my input data:

vcfs in right output format (these vcfs are filtered on coverage > 10 based on the RNA bam file)
RNA tumor bam (chr split for eval/train/test)
exons.bed file
GRCh38.d1.vd1.fa reference
I use the deepvariant version 1.9.0 sif image.

Is the header of my vcf file too extensive?

Then, I ran make_examples for train/test/eval and created the tfrecords files.
Below are the specific parameters I used:
--mode training
--customized_classes_labeler_classes_list=ref,germline,somatic
--customized_classes_labeler_info_field_name=type
--split_skip_reads=True
--channel_list='BASE_CHANNELS,insert_size'

I then shuffled these records, and obtained pbtxt files.

But when I looked deeper in my records files, I see that all the variants get label 0 (ref).
For example for one training sample and first 10000 examples:
Variants in truth vcf: 5609 germline, 73 somatic
Class label distribution:
ref (class0): 10,000 (100.0%)
TOTAL: 10,000
I also checked the overlap with the exons.bed file and this was: 5636, so this should not be the problem I think.

I am working with TCGA data, so I suppose the quality of the data should be sufficient.
Is there someone that would have any insight on this? Some help would be greatly appreciated!

With kind regards
Charlotte

[danielecook]

Is the header of my vcf file too extensive?

The header looks fine. Aren't all the contigs normally included?

But when I looked deeper in my records files, I see that all the variants get label 0 (ref).
For example for one training sample and first 10000 examples:
Variants in truth vcf: 5609 germline, 73 somatic

@chstubbe sorry, I think I missed mentioning this flag as being required as well:

--labeler_algorithm=customized_classes_labeler

Can you please give that a try and report back?

[chstubbe]

Hi Daniel

The extra flag resolved the classification! Thanks a lot!

I subdivided the vcfs upfront based on chromosomes, the originals indeed contain all contigs.
I will now shuffle the examples, but I had a question for the training script.

So I have to create a config.yaml file. I understand you need to add the paths to the pbtxt files.

• what needs to be in the pbtxt files? The paths to the examples and what else?
• what is required in the config file; I understand these need to be training parameters, but it is not entirely clear what I need to include. You recommended to use a pre-existing checkpoint of a trained model, which was the "models-DeepVariant-1.9.0-checkpoints-wgs-deepvariant.wgs.ckpt". Would it not be better to use a WES model, which I saw also exists? Do I use the parameters of this pre existing model in the config file?
• You recommended to use a RNA PON, I have one available from GTEx, but how do I integrate this?

My apologies for the rookie questions, as I am not very experienced with machine learning! :)
Thanks a lot for your help!

Kind regards
Charlotte

[danielecook]

@chstubbe how are you planning on performing shuffling?

What needs to be in the pbtxt files? The paths to the examples and what else?

Review the DeepVariant training tutorial for more details, but in short, the output should look like this:

# Generated by shuffle_tfrecords_beam.py
# class2: 124571
# class1: 173664
# class0: 43803

name: "HG001"
tfrecord_path: "OUTPUT_BUCKET/training_set.with_label.shuffled-?????-of-?????.tfrecord.gz"
num_examples: 342038
#
# --input_pattern_list=OUTPUT_BUCKET/training_set.with_label.tfrecord-?????-of-00016.gz
# --output_pattern_prefix=OUTPUT_BUCKET/training_set.with_label.shuffled

There are three mandatory fields for training: name, tfrecord_path, and num_examples. These must point to your TFRecord files and specify the total count of examples. This file is relatively straightforward to construct manually if you use an alternative method for shuffling.

**What is required in the config file? I understand these need to be training parameters, but it's not entirely clear what to include.

It depends what type of model you are using. The base config lists most parameters, and these are then overwritten when you specify a certain type of model (e.g. wes). You probably don't need to fiddle with these too much just yet - I would start with the get_deepsomatic_wgs_tumor_only_config or the get_deepsomatic_wes_tumor_only_config.

You recommended using a pre-existing checkpoint (WGS 1.9.0)—would it be better to use a WES model? Do I use the parameters of the pre-existing model in the config?**

Originally, we hadn't released a WES tumor-only checkpoint, but one is now available here: gs://deepvariant/models/DeepSomatic/1.10.0/checkpoints/wes_tumor_only/.

I recommend trying both the WGS tumor-only and WES tumor-only checkpoints if possible to see which performs better for your specific data.

You recommended using an RNA PON; I have one from GTEx, but how do I integrate this?

PON (Panel of Normals) filtering occurs after variants are called. It's not a concern during training. It is integrated into the post-processing step via the --pon_filtering=<vcf> flag.

Important! you want your training dataset to approximate your inference dataset as closely as possible. Our approach is to avoid filtering during training, allowing the model to learn the distinction between sequencing artifacts, germline variants, and somatic variants.

[chstubbe]

Hi Daniel

how are you planning on performing shuffling?

I shuffled using the apache beam script, and indeed my pbtxt files look exactly like you mention! So that is good!

I recommend trying both the WGS tumor-only and WES tumor-only checkpoints if possible to see which performs better for your specific data.

I will try to use both WGS and WES tumor only checkpoints and not edit the config files so far. I will keep you updated if I come across any bumps on the road.

Important! you want your training dataset to approximate your inference dataset as closely as possible. Our approach is to avoid filtering during training, allowing the model to learn the distinction between sequencing artifacts, germline variants, and somatic variants.

Okay I understand! Right now I am using the approach of leaving out samples from the same dataset and I will then test the full pipeline on these samples for benchmarking purposes. (n=45, left out samples=8)

Thank you!

[chstubbe]

Hi Daniel

I have a couple of questions:

1. I want to try the deepsomatic wes tumor only setting. I cannot find this setting in the config file. Is there a config file that has the settings for this?

2. Below are the number of examples I have for 45 DLBCL samples. Do you think these numbers are sufficient? Do you recommend using more samples from other tumor types, or using another method to split up the data? Do you think it is better to restrict to one tumor type or to use a pan-cancer approach?

3. A lot of parameters can be adjusted in the config file. Do you have any recommendations for our specific RNA dataset? These are the current ones:

def get_deepsomatic_wgs_tumor_only_config(config: ml_collections.ConfigDict):
  get_wgs_config(config)
  config.best_checkpoint_metric = 'tune/f1_homalt'
  config.class_weights = '1,5,50'
  # config.train_dataset_pbtxt = '/placer/prod/home/brain-genomics/danielecook/deepsomatic/ds_wgs_tumor_only/ds_wgs_tumor_only_train.dataset_config.pbtxt'
  # config.tune_dataset_pbtxt = '/placer/prod/home/brain-genomics/danielecook/deepsomatic/ds_wgs_tumor_only/ds_wgs_tumor_only_tune.dataset_config.pbtxt'
  config.train_dataset_pbtxt = '/scratch/gent/vo/000/gvo00027/projects/CST/CST2501/DEEPSOMATIC_TRAINING/PBTXT/training_dataset_filtered_shuffled_somatic.pbtxt'
  config.tune_dataset_pbtxt = '/scratch/gent/vo/000/gvo00027/projects/CST/CST2501/DEEPSOMATIC_TRAINING/PBTXT/evaluation_dataset_filtered_shuffled_somatic.pbtxt'
  config.init_checkpoint = '/scratch/gent/vo/000/gvo00027/projects/CST/CST2501/DEEPSOMATIC_TRAINING/RESOURCES/wgs_tumor_only/models-DeepSomatic-1.10.0-checkpoints-wgs_tumor_only-deepsomatic.wgs_tumor_only.ckpt'
  config.batch_size = 1024
  config.num_epochs = 5
  config.learning_rate = 0.0001
  config.tune_every_steps = 1000
  config.num_validation_examples = 0 

Thanks a lot in advance for your help and insights!

kind regards
Charlotte

[danielecook]

1. I want to try the deepsomatic wes tumor only setting. I cannot find this setting in the config file. Is there a config file that has the settings for this?

The dv_config.py has a config called deepsomatic_wes_tumor_only; During training you would specify your config as --config path/to/dv_config.py:deepsomatic_wes_tumor_only.

deepvariant/deepvariant/dv_config.py

Lines 534 to 535 in 45f2627

	elif config_name == 'deepsomatic_wes_tumor_only':
	get_deepsomatic_wes_tumor_only_config(config)

2. Below are the number of examples I have for 45 DLBCL samples. Do you think these numbers are sufficient? Do you recommend using more samples from other tumor types, or using another method to split up the data? Do you think it is better to restrict to one tumor type or to use a pan-cancer approach

The total number of training examples is low (1326). You can increase this by using downsampling, or mixing and matching different T-N pairs (for example, if you have replicates of a tumor, you could use different combinations of tumor with normal (e.g. T1-N1, T2-N1, T2-N1, T2-N2). For comparison, a recent training run we performed using WES tumor only looks like this:

# Classes:
#   class 0:    10482181
#   class 1:    39901763
#   class 2:    314495

This amounts to 0.620% somatic variants; You have a lower rate 0.048% of somatic variants, likely due to a lower TMB? [(1326/(1326+248213+2477907))*100]; It looks like you are observing approx. ~30 somatic variants per sample?

It will be difficult to train a model with such low numbers. I would look into incorporating additional cancer types (perhaps aiming for a more generalizable model).

Have you trained and evaluated any models yet? I would be curious how they perform.

3. A lot of parameters can be adjusted in the config file. Do you have any recommendations for our specific RNA dataset? These are the current ones:

• I would change config.class_weights = '1,5,50' to config.class_weights = '1,1,50'. You have plenty of germline samples to train on.

[chstubbe]

Hi Daniel

Thank you for the follow up!

So, I trained the DeepSomatic model on 37 DLBCL samples and ran the pipeline with the WGS tumor only checkpoint. I need to note that I did not use all default parameters in the config file, like you can see in my previous post.

Here are the performance results that I obtained with the som.py script:

Overall: 
Precision: 0.018 (1,8 %)
Recall: 0.1377 (13,78 %)
F1 Score: 0.0320 (3,2%)

For comparison, these were the results for running 1 DLBCL sample through the “untrained” DeepSomatic pipeline: 

Performance: 
True Positives (TP):          15
False Positives (FP):         1434
False Negatives (FN):         76
Precision:                    0.0104 (1.04%)
Recall (Sensitivity):         0.1648 (16.48%)
F1 Score:                     0.0195
 
It’s clear that results do not really improve with this training. 

Currently, I am obtaining more samples from other tumor types (lung cancer (n=83) and skin cancer (n=103)). I am also running with the wes_tumor_only checkpoint. I hope with this results can further improve.

I have sent you an email with the specific steps I used for training, in case you would like to check any pitfalls or think I can adjust anything. Are these results what you would expect? Do you think training with more samples will improve the results drastically?

Looking forward to hearing your insights!

Kind regards
Charlotte