Show presence of known mutation in RNA-seq data

Question

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:

1. Run STAR 2-pass

2. Extract alignments of interest, e.g.:

   samtools view input.bam "Chr10:18000-45500" > output.sam

3. 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?

Answer 1

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
<mpileup> 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.

Answer 2

If you only have one gene and you only need to do this once then the simplest possible workflow is to generate the aliment using STAR (optimally with the two pass method) and open the two resultant .bam files in IGV, the coverage column at the stop above your alignment track should clearly show the counts of reference and non-reference supporting reads. The IGV documentation shows a screen shot the pop-up coverage indicator and allele specific read count.

Answer 3

SCmut is a method to detect cell-level mutation from single-cell RNA-sequencing, published by Vu et al. (2019)

https://github.com/nghiavtr/SCmut