strand bias artifacts

[katestan]

We have generated targeted NanoSeq data using the human core exome twist panel (~1000X raw coverage, 200X duplex coverage). We are using the same sample as reference with VAF cutoff of 0.1 for somatic variant calling (var_v=0.1). The pipeline completed without errors using these parameters:

// dsa parameters
params.dsa_d = 2
params.dsa_q = 30
params.dsa_M = 0
// variantcaller parameters
params.var_a = 50
params.var_b = 0
params.var_c = 0.02
params.var_d = 2
params.var_f = 0.9
params.var_i = 1.0
params.var_m = 8
params.var_n = 3
params.var_p = 0
params.var_q = 45
params.var_r = 144
params.var_v = 0.1
params.var_x = 8
params.var_z = 25
// indel parameters
params.indel_rb = 2
params.indel_t3 = 136
params.indel_t5 = 8
params.indel_z = 25
params.indel_v = params.var_v
params.indel_a = params.var_a
params.indel_c = params.var_c

I am noticing some strange results in the final $sample.vcf.gz output by the NanoSeq pipeline.
Certain variants appear in all 8 of my samples. When I inspect these variants in IGV using the $sample.opmk.bam as input I notice that there is strand bias, i.e. all of the alternate alleles appear on the reverse strand. There are also multiple alternative alleles visible on IGV at the locus (all appearing on the reverse strand) that do not appear in my $sample.vcf.gz because they probably don't meet 2x2 duplex criteria.

For one of my $sample.opkm.bam, the IGV breakdown at this locus looks like:
Total count: 1247
A : 39 (3%, 0+, 39- )
T : 25 (2%, 0+, 25- )
C : 268 (21%, 0+, 268- )
G : 915 (73%, 610+, 305- )
N : 0

And the $sample.vcf.gz from the NanoSeq pipeline looks like:
chr5 154834880 . G C . PASS TRI=ACT>G;TIMES_CALLED=6;TYPE=snv;DUPLEX_VAF=0.0535714;BAM_VAF=0.00502513;BAM_VAF_BQ10=0.0840108;DEPTH_NORM_FWD=591;DEPTH_NORM_REV=60;DEPTH_FWD=8;DEPTH_REV=6.5;DUPLEX_COV=112;BAM_MUT=3;BAM_COV=597;BAM_MUT_BQ10=62;BAM_COV_BQ10=738;RB=chr5:154834649:154834922:ACC:TGA,chr5:154834707:154834907:GTA:TCT,chr5:154834718:154834922:CTA:GGT,chr5:154834720:154834908:TTT:TAG,chr5:154834729:154834964:GCA:GCT,chr5:154834731:154834930:ATC:GTG;QPOS=42,27,42,28,84,50;DPLX_ASXS=117,118,116,116,118,116;DPLX_CLIP=0,0,0,0,0,0;DPLX_NM=1,1,1,1,1,1;BULK_ASXS=118,118,118,118,118,118;BULK_NM=0,0,0,0,0,0

Questions:

(1) Can you still get a 2x2 duplex call if all of the supporting reads are reverse orientation? Could that happen if both PCR duplicates on positive strand are reverse AND both PCR duplicates on the negative strand are reverse?

(2) Is there a filter in the NanoSeq pipeline that can help me filter out this kind of strand bias artefact? I see that there is
DEPTH_FWD="Read bundle forward reads depth"
DEPTH_REV="Read bundle reverse reads depth"
What do these mean exactly and are they used for filtering?

[fa8sanger]

That “bias” is not a bias, it’s just how read bundles work. If you look at EDF 1d it may be clearer (Abascal et al 2021). I know its counterintuitive but the two original DNA strands will map in the same orientation, the only difference will be that one strand will be that the forward from one strand will be the first read (F1R2 pair) while the forward from the other strand will be the second read (F2R1 pair). So, supporting each mutation you will only see F1 + F2 (or R1 + R2). There can only be a mixture F and R when the two reads are overlapping.

In your cases the mutation has been called 6 times. When a mutation is called multiple times you can expect to see a mixture of Fwd and Rev. However, hybridisation capture (exome panel) introduces many biases in mapping orientation, and usually, in regions flanking a bait you only see one orientation.

In the final VCF, the fields
DEPTH_FWD="Read bundle forward reads depth"
DEPTH_REV="Read bundle reverse reads depth”
should be DEPTH_F1R2 and DEPTH_F2R1, it’s an old mistake

[katestan]

Thank you for your reply.

OK I understand that the two original DNA strands map in the same orientation per read bundle.

To be clear: the DEPTH_F1R2 and DEPTH_F2R1 is the average number of F1R2 and F2R1 per bundle before collapsing the reads? Otherwise it doesn't make sense, because they are not integers.

I am still unsure of how to handle these variants, as they are clearly artefacts but still appear in my VCF. Or are you suggesting that they are real given that it's normal to see only one orientation for the ALT in regions flanking a bait with capture hybridisation?

[fa8sanger]

The VCF you are looking at is the output of snv_merge_and_vaf_calc.R. This script summarises the information in muts.vcf and indel.vcf. If a mutation is seen multiple times (in different read bundles), in this VCF we provide either the average (DEPTH_FWD/REV) or just the list of individual values (e.g. QPOS). If you want to look at the individual values for DEPTH_FWD/REV check muts.vcf.

Why do you think this is an artefact? There can be artefacts for different reasons: 1) low AS-XS (not the case); 2) read position bias (in QPOS, not this case; typically associated with nearby germline indels); 3) SNP leakage or contamination from another human or even from non human DNA. For the latter, have you compared the burdens before and after masking?