Dockerized process for variant calling on RNA-Seq data.
The steps for the workflow are:
- Call a small set of variants using Freebayes (v9.9.2-27-g5d5b8ac)
- Annotate variants using SnpEff (4.3r, GRCh38.86)
docker run --rm -v /data/:/data \
hbeale/mini_var_call \
/data/sampleID.sortedByCoord.md.bam \
/data/ref/GCA_000001405.15_GRCh38_no_alt_analysis_set.fa
The docker can be obtained from docker.io with the command.
docker push hbeale/mini_var_call
The source code is available at https://github.com/UCSC-Treehouse/rna_variant_call.
Because there are few of these variants, it is common to have no variants reported for a given sample.
Because these are cannonical variants that may be supported by few reads, Freebayes is set to maximum sensitivity. Review stats and bam alignment to determine whether you think they are real.
This type of error is expected and is not of concern:
Could not find any mapped reads in target region chr2:29222346..29222347
Could not find any mapped reads in target region chr2:29222406..29222407
The docker takes as arguments the locations of a coordinate-sorted bam and the GCA_000001405.15_GRCh38_no_alt_analysis_set.fa. The fasta file must have a correponding index file. The bam file may have a corresponding index file (bai). It will greatly increase the speed (15 min -> 3 min for the file I tested time on), but indexing takes long enough that it doesn't doesn't improve the overall time required to index the bam file for this purpose alone.
The output file will be in the mounted directory containing test.bam.
sampleID.mini.ann.vcf
The following commands can be used to review variants with a quality score above zero. We expect this to contain false positive variants, so review the full data if the preliminary results are of interest.
vcf=sampleID.sortedByCoord.md.mini.ann.vcf
cat $vcf | grep -v ^# | \
awk '$6 !~ /E/ && $6 > 0 { print }' | \
sed 's/^.*;EFF=\([^)]*\)).*/\1/' | \
cut -f6,4 -d"|" --output-delimiter " " | awk '{ print $2 " " $1}'
Generates this type of output:
JAK2 R683S
NRAS G12D
The code repository contains data that can be run to confirm that your installation is working correctly.
Commands to call variants:
cd /data/tmp
git clone https://github.com/UCSC-Treehouse/rna_variant_call.git
docker run --rm -v /data/:/data \
hbeale/mini_var_call \
/data/tmp/rna_variant_call/test/test.bam \
/data/ref/GCA_000001405.15_GRCh38_no_alt_analysis_set.fa
Commands to review variants with a quality score above zero:
vcf=/data/tmp/rna_variant_call/example/Jak2Example.mini.ann.vcf
cat $vcf | grep -v ^# | \
awk '$6 !~ /E/ && $6 > 0 { print }' | \
sed 's/^.*;EFF=\([^)]*\)).*/\1/' | \
cut -f6,4 -d"|" --output-delimiter " " | awk '{ print $2 " " $1}'
Expected output
JAK2 R683S
00:00:00 SnpEff version SnpEff 4.3r (build 2017-09-06 16:41), by Pablo Cingolani
00:00:00 Command: 'ann'
00:00:00 Reading configuration file 'snpEff.config'. Genome: 'GRCh38.86'
00:00:00 Reading config file: /data/snpEff.config
00:00:00 Reading config file: /root/snpEff/snpEff.config
00:00:01 done
00:00:01 Reading database for genome version 'GRCh38.86' from file '/root/snpEff/./data/GRCh38.86/snpEffectPredictor.bin' (this might take a while)
00:00:35 done
00:00:35 Loading interactions from : /root/snpEff/./data/GRCh38.86/interactions.bin
00:00:52 Interactions: 1793688 added, 0 skipped.
00:00:52 Building interval forest
00:01:00 done.
00:01:00 Genome stats :
#-----------------------------------------------
# Genome name : 'Homo_sapiens'
# Genome version : 'GRCh38.86'
# Genome ID : 'GRCh38.86[0]'
# Has protein coding info : true
# Has Tr. Support Level info : true
# Genes : 58051
# Protein coding genes : 20423
#-----------------------------------------------
# Transcripts : 198002
# Avg. transcripts per gene : 3.41
# TSL transcripts : 166906
#-----------------------------------------------
# Checked transcripts :
# AA sequences : 0 ( 0.00% )
# DNA sequences : 163038 ( 82.34% )
#-----------------------------------------------
# Protein coding transcripts : 94384
# Length errors : 13357 ( 14.15% )
# STOP codons in CDS errors : 51 ( 0.05% )
# START codon errors : 11250 ( 11.92% )
# STOP codon warnings : 7169 ( 7.60% )
# UTR sequences : 91549 ( 46.24% )
# Total Errors : 23099 ( 24.47% )
#-----------------------------------------------
# Cds : 704604
# Exons : 1182163
# Exons with sequence : 1182163
# Exons without sequence : 0
# Avg. exons per transcript : 5.97
# WARNING! : Mitochondrion chromosome 'MT' does not have a mitochondrion codon table (codon table = 'Standard'). You should update the config file.
#-----------------------------------------------
# Number of chromosomes : 524
# Chromosomes : Format 'chromo_name size codon_table'
# 'HSCHR1_2_CTG3' 248975002 Standard
# 'HSCHR1_1_CTG31' 248973653 Standard
# 'HSCHR1_2_CTG31' 248971826 Standard
...
# 'KI270742.1' 186739 Standard
# 'GL000205.2' 185591 Standard
# 'GL000195.1' 182896 Standard
# 'KI270736.1' 181920 Standard
# 'KI270733.1' 179772 Standard
... and so on
00:01:08 done.
00:01:08 Logging
00:01:09 Checking for updates...
00:01:10 Done.