Show presence of known mutation in RNA-seq data We have RNA-seq fastq data from control (WT) patients and a patient with a point mutation at a known location in one gene.
I'd like to retrieve the reads aligning to that gene and show the presence of the mutation.

I can think of 2 strategies:

A)

Alignment and visualization, based on GATK:


Run STAR 2-pass
Extract alignments of interest, e.g.:

samtools view input.bam "Chr10:18000-45500" &gt; output.sam
Visualize:


IGV
IGB
ASCIIGenome (see also this SE answer)



How can I count the reads showing mutation vs wild type?



B)

The second approach would simply look for the count numbers mapping (e.g. salmon) to the transcript ID(s) corresponding to mutation and WT, based on this SE answer.

Is there a better way to do this?
For counting reads I use mpileup, e.g. samtools mpileup --reference hg38.fa -r Chr10:18000-45500 input.bam, which will give base-resolution coverage for a BAM file.

I've written my own script to process mpileup output and make it easier to understand. By default it reports read coverage as a proportion of total coverage, but this can be modified by using the -counts command line argument:

$ samtools mpileup -r tig00018708_tig00000379:210665-240664:2107-2130 --reference trimMmerged_tig00018708_other.fasta local_Sampled_50M_vs_trimMmerged_tig00018708_other.bam  | /bioinf/scripts/readstomper.pl -counts

[mpileup] 1 samples in 1 input files
&lt;mpileup&gt; Set max per-file depth to 8000
Assembly,Position,Coverage,ref,cR,pR,A,C,G,T,d,i,InsMode
tig00018708_tig00000379:210665-240664,2107,10,C,10,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2108,11,G,11,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2109,11,T,11,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2110,11,T,11,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2111,12,G,12,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2112,14,G,14,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2113,13,C,13,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2114,13,G,13,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2115,13,C,2,0.15,0,0,0,11,0,0
tig00018708_tig00000379:210665-240664,2116,15,T,14,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2117,15,G,14,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2118,18,A,17,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2119,19,T,18,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2120,19,C,18,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2121,18,A,7,0.41,0,0,10,0,0,0
tig00018708_tig00000379:210665-240664,2122,18,G,18,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2123,18,T,18,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2124,18,C,18,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2125,18,C,9,0.50,0,0,0,9,0,0
tig00018708_tig00000379:210665-240664,2126,18,C,14,0.78,2,0,1,1,0,0
tig00018708_tig00000379:210665-240664,2127,18,C,18,1.00,0,0,0,0,0,0
tig00018708_tig00000379:210665-240664,2128,18,T,15,0.83,0,0,3,0,0,0
tig00018708_tig00000379:210665-240664,2129,17,A,13,0.76,0,0,4,0,0,0
tig00018708_tig00000379:210665-240664,2130,16,A,9,0.56,0,6,1,0,0,0


Note that only non-reference coverage is reported in the A/C/G/T/d/i columns. I have found that this makes it easier for subsequent data analysis, but your mileage may vary.