Mutant Proteoform Prediction

Call somatic DNA and RNA variants from long-read DNA and RNA sequencing, link them, and translate the assembled transcripts to mutant proteoforms and peptides.

This pipeline is Exacto’s primary use case: identify somatic DNA variants in a case sample (for example a tumor) against a matched control, assemble the case’s transcriptome from long RNA reads, call RNA variants on them, link each RNA variant to the DNA variant that explains it, and translate the transcripts to proteoforms whose every amino acid records the variants it came from.

External tools you’ll need alongside Exacto:

Workflow

%%{init: {'securityLevel': 'loose', 'flowchart': {'rankSpacing': 20, 'nodeSpacing': 40, 'subGraphTitleMargin': {'top': 5, 'bottom': 20}}}}%%
flowchart TD
    subgraph DNA["<span style='white-space:nowrap;color:#414a4c'>DNA variant calling</span>"]
        direction LR
        DNAIN[/"<span style='white-space:nowrap;color:#414a4c'>Tumor + normal DNA BAMs</span>"/] --> CALLDNA(call-somatic-dna-vars)
        CALLDNA --> DNAOUT[/"<span style='white-space:nowrap'>Somatic DNA variants TSV</span>"/]
    end

    subgraph RNA["<span style='white-space:nowrap;color:#414a4c'>Variant-aware transcriptome assembly</span>"]
        direction LR
        RNAIN[/"<span style='white-space:nowrap;color:#414a4c'>Tumor RNA reads BAM</span>"/] --> RECON(cluster, correct, consensus, stitch)
        RECON --> CALLRNA(call-rna-transcript-vars)
        CALLRNA --> RNAOUT[/"<span style='white-space:nowrap;color:#414a4c'>RNA variants + transcript structures</span>"/]
    end

    DNAOUT --> INT(integrate-vars)
    RNAOUT --> INT
    INT --> INTOUT[/"<span style='white-space:nowrap;color:#414a4c'>Integrated variants TSV</span>"/]

    RNAOUT --> TR(translate-transcripts)
    DNAOUT --> TR
    INTOUT --> TR
    TR --> PROT[/"<span style='white-space:nowrap;color:#414a4c'>Proteoforms TSV + FASTA</span>"/]

    PROT --> CPV(call-peptide-vars)
    CPV --> PEPS[/"<span style='white-space:nowrap;color:#414a4c'>Mutant peptides TSV + FASTA</span>"/]

    style DNA fill:#ffffff,stroke:#bbbbbb
    style RNA fill:#ffffff,stroke:#bbbbbb

    click CALLDNA href "../cli/call-somatic-dna-vars.html" "View call-somatic-dna-vars docs" _self
    click RECON href "variant-aware-transcriptome-assembly.html" "View the variant-aware transcriptome assembly pipeline" _self
    click CALLRNA href "../cli/call-rna-transcript-vars.html" "View call-rna-transcript-vars docs" _self
    click INT href "../cli/integrate-vars.html" "View integrate-vars docs" _self
    click TR href "../cli/translate-transcripts.html" "View translate-transcripts docs" _self
    click CPV href "../cli/call-peptide-vars.html" "View call-peptide-vars docs" _self

    classDef linked text-decoration:underline;
    class CALLDNA,RECON,CALLRNA,INT,TR,CPV linked

Step 1. Align the DNA reads

Align the tumor and normal long-read DNA to the reference genome, then sort and index. --cs is required: Exacto reads variants from the cs tag.

# Tumor
minimap2 -ax map-hifi --cs --eqx -Y -L \
    reference_genome.fasta tumor_dna.fastq.gz \
    | samtools sort -o tumor_dna.sorted.bam
samtools index tumor_dna.sorted.bam

# Normal: repeat with normal_dna.fastq.gz -> normal_dna.sorted.bam

Step 2. Call somatic DNA variants

Call the variants of the tumor that the normal does not carry. Use --preset pb for PacBio HiFi reads or --preset ont for Oxford Nanopore reads:

exacto call-somatic-dna-vars \
    --bam-file tumor_dna.sorted.bam \
    --bai-file tumor_dna.sorted.bam.bai \
    --control-bam-files normal_dna.sorted.bam \
    --control-bai-files normal_dna.sorted.bam.bai \
    --fasta-file reference_genome.fasta \
    --preset pb \
    --output-tsv-file tumor_somatic_dna_variants.tsv

Germline variants are called the same way with call-germline-dna-vars, from one BAM file.

Step 3. Assemble the tumor transcriptome

Run the Variant-aware Transcriptome Assembly pipeline on the tumor’s long RNA reads, Steps 1 to 9: two passes of clustering with correction in between, one consensus sequence per transcript, and reference stitching of degraded ends. This step uses the files:

  • pass2/tumor_pass2_exacto_rna_consensus.tsv (the consensus sequences and their reads)
  • stitched/tumor_exacto_stitched_transcripts.tsv (the stitched sequences)
  • stitched/tumor_exacto_stitched_transcripts.sorted.bam (the stitched sequences, aligned)

Step 4. Call RNA variants

exacto call-rna-transcript-vars \
    --bam-file stitched/tumor_exacto_stitched_transcripts.sorted.bam \
    --reference-genome-fasta-file reference_genome.fasta \
    --reference-gene-annotation-file gencode.gtf.gz \
    --reference-gene-annotation-source gencode \
    --reference-gene-annotation-assembly hg38 \
    --reference-gene-annotation-version v45 \
    --dna-variants-tsv-files tumor_somatic_dna_variants.tsv \
    --output-dir rna_variants/ \
    --output-prefix tumor

--dna-variants-tsv-files is optional: an RNA variant that matches a DNA variant (same chromosomes, operations, type and sequence, each position within --dna-variant-match-buffer bases) gets that variant’s origin in the origin column.

Step 6. Translate the transcripts to proteoforms

translate-transcripts takes the transcript sequences from the consensus TSV, replaced by the stitched sequences when --stitched-transcripts-tsv-file is given. The call-rna-transcript-vars tables must then come from the stitched alignment, as in Step 4:

exacto translate-transcripts \
    --assembled-transcript-model-alignments-tsv-file rna_variants/tumor_exacto_assembled_transcript_model_alignments.tsv \
    --assembled-transcript-variants-tsv-file rna_variants/tumor_exacto_assembled_transcript_variants.tsv \
    --dna-variants-tsv-file tumor_somatic_dna_variants.tsv \
    --integrated-variants-tsv-file tumor_integrated_variants.tsv \
    --rna-consensus-tsv-file pass2/tumor_pass2_exacto_rna_consensus.tsv \
    --stitched-transcripts-tsv-file stitched/tumor_exacto_stitched_transcripts.tsv \
    --strategy longest_orf \
    --output-dir translation/ \
    --output-prefix tumor

Transcripts assembled outside exacto are translated by passing --assembled-transcripts-support-tsv-file instead of --rna-consensus-tsv-file.

Step 7. Call mutant peptides

Cut each proteoform into peptides of 8 to 11 amino acids (--min-k, --max-k) and keep those that overlap a mutant amino acid and do not occur in the reference proteome:

exacto call-peptide-vars \
    --proteoforms-tsv-file translation/tumor_exacto_proteoforms.tsv \
    --reference-fasta-file reference_proteome.fasta \
    --output-tsv-file tumor_peptide_variants.tsv \
    --output-fasta-file tumor_peptide_variants.fasta
Warning

call-peptide-vars builds a k-mer index of the reference proteome, which takes about 16 GB of memory for a human proteome. Budget for it when running several samples at once.

Outputs

File Produced by Description
tumor_somatic_dna_variants.tsv call-somatic-dna-vars Somatic DNA variants
rna_variants/tumor_exacto_assembled_transcript_variants.tsv call-rna-transcript-vars RNA variants per transcript
rna_variants/tumor_exacto_assembled_transcript_model_alignments.tsv call-rna-transcript-vars Each transcript’s structure against the genome
tumor_integrated_variants.tsv integrate-vars RNA variants linked to the DNA variants that explain them
translation/tumor_exacto_proteoforms.tsv translate-transcripts Proteoforms with their mutant amino acid intervals and source variants
translation/tumor_exacto_proteoforms.fasta translate-transcripts Proteoform sequences
translation/tumor_exacto_proteoform_nucleotides.tsv translate-transcripts Per-nucleotide provenance of each proteoform
tumor_peptide_variants.tsv call-peptide-vars Mutant peptides