How to correctly call a VCF file using damaged DNA? (IonTorrent & FFPE)

Question

EDIT: I am updating this question to make it more specific to my issue.

For context - original question prior to edit:

How do I obtain a deamination metric when doing the variant calling using the IonTorrent variant caller, and secondly, how do I correct my called variants for deamination to ensure that these don't provide false positives for the downstream analysis?

Background

I'm calling variants using the IonTorrent TorrentSuite on DNA which has been sequenced from formalin-fixed paraffin-embedded (FFPE) tissue. This has a major issue in that without addition of uracil-N-glycosylase, some of the Ts in the original DNA are deaminated to Cs, which upon sequencing and calling variants can show up as mutations, either as T>C transitions or G>A (from the opposite strand, due to PCR in the library prep). I do not have any idea how long these samples were stored without UNG before sequencing.

TVC (Torrent Variant Caller) gives a deamination metric (essentially, sum of T>C and G>A variants over all variants called), and for our samples, the highest value seen is ~0.92. Naively postprocessing the variants show that for these samples, C>T/T>C transitions overwhelm the remaining variants among my samples.

Application

My use case is this: these are medical samples, which have been inspected by a pathologist (hence the FFPE treatment), and I want to determine which variants are predictive of outcome, hence I have two potentially contradictory goals: reduce false positives and capture the rarer variants which may hold predictive power.

Question

• Given the IonTorrent variant calling pipeline (sequencing > BAM file > TVC > VCF file with deamination statistic):
  • (Q1) Is there a way of correcting the output VCF to remove transition and deamination errors?
  • (Q2) Alternatively, is there a set of filters to use in bcftools to reduce the effect of the errors on how variants are called?

Answer 1

Sites that were deaminated are often called abasic sites.

Q1: How to obtain DNA damage metrics?

You might want to check out mapDamage - this software is used to perform these sorts of analyses and produce plots of types of DNA damage along sequencing reads.

Q2: How to correct for abasic sites and other DNA damage errors?

You can use your choice of kmer-based or alignment based error correction software. The errors will be random and if the sequencing depth is sufficient the erroneous kmers at the ends of reads will be corrected. However, depending on what you are trying to do downstream with your reads you may be better served by not correcting the reads, then dealing with transitions/damage later.

Answer 2

Having inspected one of my samples (with estimated deamination metric, the proportion of called SNPs consisting of C>T or G>A transitions, of 0.837), this looks like a straightforward filtering procedure. Looking at the minor allele fraction (AF) and the quality score (-10 log10(p-value)) per type of SNP: it can be seen that the C>T and G>A transitions likely caused by deamination (and PCRing up the deaminated DNA to get the libraries) all cluster in the bottom right of the plot - low AF and low QUAL, so setting thresholds to the left (higher QUAL) and above (higher AF) the groups will remove these from the called variants.