This is our protocol for calling putative somatic mutations from RNA-seq data.

Briefly our method can be divided into three parts,

  1. Deriving raw read-counts at all sites from the RNA-seq BAM files.
  2. Constructing site-specific error models from all the read-count data.
  3. Using the site-specific error model to identify putative somatic mutations.

More detail on how to perform each step is below.

Deriving raw read-counts at all sites from the RNA-seq BAM files

To accomplish this task we use three programs,

samtools

mpileup2readcounts

tabix

Use a recent version of samtools (ver > 1.5).

On each BAM file, run the following command, in this example the sample used is SRR106:

samtools mpileup -q 20 -Q 20 -f Homo_sapiens_assembly19.fasta SRR106.bam | mpileup2readcounts SRR106 | awk '{ if(NR==1 || \$5 + \$6 > 10) print }' | bgzip > SRR106_readcounts.txt.gz ; tabix -b 2 -e 2 -S 1 SRR106_readcounts.txt.gz

This will create a readcount file for each sample, in this case this is named as SRR106_readcounts.txt.gz. Store these files in a convenient location.

Constructing site specific error models from all the read-count data

For this step we will need to download and install bcall from here.

To generate site specific error models you can frist divide your set of samples into as many subsets as you would like or just run them as one huge subset. Dividing into smaller subsets lets you run these subsets in parallel and hence save time. The division of samples is arbitrary and not necessarily based on any biological input,

bcall prior-dump-fixed subset00 ./prior_dumps/subset00.dump SeqCap_EZ_Exome_v3_primary.bed.gz
bcall prior-dump-fixed subset01 ./prior_dumps/subset00.dump SeqCap_EZ_Exome_v3_primary.bed.gz

Here subset00 is the sample subset file and the ./prior_dumps/subset00.dump is the name of the output file where you want to store the dump. The SeqCap_EZ_Exome_v3_primary.bed.gz file restricts us to look at regions within the BED file, in this case these are regions covered by an exome capture reagent.

The first few lines of the sample subset file (for e.g subset00 in this case) looks like,

SRR1068687 ./gtex-readcounts/2016-05-27/SRR1068687_readcounts.txt.gz
SRR1068788 ./gtex-readcounts/2016-05-27/SRR1068788_readcounts.txt.gz
SRR1068808 ./2016-05-27/SRR1068808_readcounts.txt.gz
SRR1068832 ./gtex-readcounts/2016-05-27/SRR1068832_readcounts.txt.gz

The first column of this file is the sample name and the second column is the path to the readcounts file for this sample.

To merge all the subsets after running them individually, run

bcall prior-merge prior_dump_list.tsv merged.dump

Here merged.dump is the file which contains the output of all the subsets merged together.

The prior_dump_list.tsv file looks like,

subset00 ./prior_dumps/subset00.dump
subset01 ./prior_dumps/subset01.dump
subset02 ./prior_dumps/subset02.dump
subset03 ./prior_dumps/subset03.dump
subset04 ./prior_dumps/subset04.dump

The first column of this file is the name of the subset and the second column is the path to the priors for this subset.

Using the site-specific error model to identify putative somatic mutations

To identify somatic mutations using the merged prior file merged.dump generated in the last step, run

bcall call-using-merged sample_subsets/subset1 merged.dump

Here, the samples to identify somatic mutations in are specified in the subset1 file, this looks like the file we used in the earlier step,

SRR1068687 ./gtex-readcounts/2016-05-27/SRR1068687_readcounts.txt.gz
SRR1068788 ./gtex-readcounts/2016-05-27/SRR1068788_readcounts.txt.gz
SRR1068808 ./2016-05-27/SRR1068808_readcounts.txt.gz
SRR1068832 ./gtex-readcounts/2016-05-27/SRR1068832_readcounts.txt.gz

The first column is the name of the sample and the second column is the path to the readcounts for that sample. By using the readcount information and merged prior information (i.e the site-specific error model), we can identify putative somatic mutations.

These putative somatic mutations are then run through a series of steps which help remove likely artefacts.