In this case study, we show an example of running DeepSomatic on WGS tumor-only data. We use HCC1395 as an example for this case study.
For this case-study, we use HCC1395 as an example. We run the analysis on chr1
that we hold out during training.
Please see the metrics page for details on runtime and data.
For accurate tumor-only calling, we use the allele-frequency channel that uses
1000 genomes variant calls using DeepVariant to filter out germline variants
during inference. Currently, the default VCF is set to variant calls against
GRCh38 reference. If you want to customize this to your VCF then please do
so by using --population_vcfs
parameter.
Docker will be used to run DeepSomatic and hap.py,
We will be using GRCh38 for this case study.
BASE="${HOME}/deepsomatic-wgs-tumor-only-case-study"
# Set up input and output directory data
INPUT_DIR="${BASE}/input/data"
OUTPUT_DIR="${BASE}/output"
## Create local directory structure
mkdir -p "${INPUT_DIR}"
mkdir -p "${OUTPUT_DIR}"
mkdir -p "${OUTPUT_DIR}/sompy_output"
# Download bam files to input directory
HTTPDIR=https://storage.googleapis.com/deepvariant/deepsomatic-case-studies/deepsomatic-chr1-case-studies
# Download the reference files
curl ${HTTPDIR}/GCA_000001405.15_GRCh38_no_alt_analysis_set.chr1.fna > ${INPUT_DIR}/GCA_000001405.15_GRCh38_no_alt_analysis_set.chr1.fna
curl ${HTTPDIR}/GCA_000001405.15_GRCh38_no_alt_analysis_set.chr1.fna.fai > ${INPUT_DIR}/GCA_000001405.15_GRCh38_no_alt_analysis_set.chr1.fna.fai
# Download the bam file
curl ${HTTPDIR}/HCC1395_wgs.tumor.chr1.bam > ${INPUT_DIR}/HCC1395_wgs.tumor.chr1.bam
curl ${HTTPDIR}/HCC1395_wgs.tumor.chr1.bam.bai > ${INPUT_DIR}/HCC1395_wgs.tumor.chr1.bam.bai
# Download truth VCF
DATA_HTTP_DIR=https://storage.googleapis.com/deepvariant/deepsomatic-case-studies/SEQC2-S1395-truth
wget -P ${INPUT_DIR} "${DATA_HTTP_DIR}"/High-Confidence_Regions_v1.2.bed
wget -P ${INPUT_DIR} "${DATA_HTTP_DIR}"/high-confidence_sINDEL_sSNV_in_HC_regions_v1.2.1.merged.vcf.gz
wget -P ${INPUT_DIR} "${DATA_HTTP_DIR}"/high-confidence_sINDEL_sSNV_in_HC_regions_v1.2.1.merged.vcf.gz.tbi
DeepVariant pipeline consists of 3 steps: make_examples_somatic
, call_variants
, and
postprocess_variants
. You can run DeepSomatic with one command using the
run_deepvariant
script.
BIN_VERSION="1.8.0"
sudo docker pull google/deepsomatic:"${BIN_VERSION}"
sudo docker run \
-v ${INPUT_DIR}:${INPUT_DIR} \
-v ${OUTPUT_DIR}:${OUTPUT_DIR} \
google/deepsomatic:"${BIN_VERSION}" \
run_deepsomatic \
--model_type=WGS_TUMOR_ONLY \
--ref=${INPUT_DIR}/GCA_000001405.15_GRCh38_no_alt_analysis_set.chr1.fna \
--reads_tumor=${INPUT_DIR}/HCC1395_wgs.tumor.chr1.bam \
--output_vcf=${OUTPUT_DIR}/HCC1395_deepsomatic_output.vcf.gz \
--sample_name_tumor="HCC1395Tumor" \
--num_shards=$(nproc) \
--logging_dir=${OUTPUT_DIR}/logs \
--intermediate_results_dir=${OUTPUT_DIR}/intermediate_results_dir \
--use_default_pon_filtering=true \
--regions=chr1
By using --use_default_pon_filtering=true
the somatic variants will be
filtered using the default PON vcf that contains variant calls from dbSNP,
gnomAD and 1000 genomes. If you plan to customize post-filtering, then you
can set this parameter to false
and use custom filtering.
NOTE: If you want to run each of the steps separately, add --dry_run=true
to the command above to figure out what flags you need in each step. Based on
the different model types, different flags are needed in the make_examples
step.
--intermediate_results_dir
flag is optional. By specifying it, the
intermediate outputs of make_examples_somatic
and call_variants
stages can be found in the directory.
sudo docker pull pkrusche/hap.py:latest
# Run hap.py
sudo docker run \
-v ${INPUT_DIR}:${INPUT_DIR} -v ${OUTPUT_DIR}:${OUTPUT_DIR} \
pkrusche/hap.py:latest \
/opt/hap.py/bin/som.py \
-N ${INPUT_DIR}/high-confidence_sINDEL_sSNV_in_HC_regions_v1.2.1.merged.vcf.gz \
${OUTPUT_DIR}/HCC1395_deepsomatic_output.vcf.gz \
-r ${INPUT_DIR}/GCA_000001405.15_GRCh38_no_alt_analysis_set.chr1.fna \
-o ${OUTPUT_DIR}/sompy_output/deepsomatic.chr1.sompy.output \
--feature-table generic \
-R ${INPUT_DIR}/High-Confidence_Regions_v1.2.bed \
-l chr1
The output:
type total.truth total.query tp fp fn unk ambi recall recall_lower recall_upper recall2 precision precision_lower precision_upper na ambiguous fp.region.size fp.rate
0 indels 133 368 86 282 47 0 0 0.646617 0.562921 0.724009 0.646617 0.233696 0.192656 0.278905 0 0 248956422 1.132728
1 SNVs 3440 4813 2937 1876 503 0 0 0.853779 0.841676 0.865287 0.853779 0.610222 0.596381 0.623931 0 0 248956422 7.535455
5 records 3573 5181 3023 2158 550 0 0 0.846068 0.833956 0.857619 0.846068 0.583478 0.570011 0.596852 0 0 248956422 8.668184