Detect mutation context in a read of a sam file

Question

After sequencing (Illumina) of some DNA, I generated a sam file through alignment of a fastq file (using Bowtie2). Instead of using variant calling programs, I want to know the specific variances, compared to the reference genome, in each specific read.

For example, I'm only interested in C to A mutations. I would like to include reads that contain MORE than 5 C to A mutations. Afterwards I would like to analyze the specific sequence context of that mutation and for doing that I need to know the 3 nucleotides before and after the C to A mutation, which gives a context of NNNCNNN.

I looked at samtools mpileup or bcftools, but I can't find the exact thing that I want to analyze. Any ideas?

Thank you!

Answer 1

That kind of custom fiddling with reads and variants is very cumbersome, non-standard and also error-prone. Do a standard variant calling pipeline and then filter for the mutations that you want. Then extract the variant position (so the coordinates) and get the variant context from the reference genome. Using individual reads for this is not advisable as only a properly called variant is reliable. The power of NGS is the validation of an observation by redundancy (=many reads covering a variant). With individual reads your variant can simply be noise, there is hardly a way to tell. Alignment will confidently map the read to the "context" as you call it, so in my opinion there is not really an advantage of calling the "context" directly from individual reads. You could better extract every confident C>A mutation, then get the three bases up/downstream and then make some kind of position frequency matrix to create a motif.

Example, just for illustration, imagine these were for variants (A in the middle) with three bases up/downstream:

CTAATTG
CTCATAG
CTTATAG
CCTAGTG

The motif would then be this, you can easily see whether there is some conservation in the "context" of your variant.ad

The variant calling you could do with an end-to-end pipeline such as nf-core sarek, this will eventually give you a VCF file(s) with the variant coordinates. From there simply filter for the variants you want, get the coordinate, add three bases up/downstream to start/end and then use something like bedtools getfasta to get the sequences from the reference genome.

Answer 2

If you want to look at a few, specific reads manually, you can do this on VarSome. We have recently released a new feature that allows you to enter single NGS or Sanger reads and get the list of variants detected in that read. For example, here:

In this example, there was one variant identified when comparing to the human hg38 genome assembly:

• CYP24A1(NM_000782.5):c.1186C>T which is at chr20:54158136 G⇒A (note that the CYP24A1 gene is on the - strand which is why the C to T variant in HGVS notation is a G to A in the genome). You can then click on the variant link(s) to investigate the variant's effect and annotation.

Obviously, this isn't practical for many reads, but it can be handy if you just want to look at a few manually. More importantly, NGS is very error prone (which is why we sequence at high depth), and you should absolutely not assume that just because a variant is present in a single read it was also actually present in the source DNA you sequenced.

---

Disclaimer: I work for the company that builds VarSome.