Question: Workflow Or Tutorial For Snp Calling?
gravatar for Matthieu
9.7 years ago by
Nova Scotia, Canada
Matthieu460 wrote:

I am looking for a good workflows, readings or tutorial for SNP calling. I read some other posts on this topic, but I would like a more detailed explanation. Population genomics and sequence data are new to me (I have a general CS and biology background). It might just be me, but these tools are not as straightforward or as documented as I'd like. Any links or explanations would be good!

So far, my situation is as follows:

  • I have Illumina sequence reads for a highly polymorphic species
  • I aligned these reads using BWA against a reference genome with default parameters, but I am not sure if I should change parameters (if so, which ones?) due to the highly polymorphic data
  • I am unsure of the next step, I will probably be using SAMtools or GATK... I tried making an mpile up but got really confused after that.
  • I should also be accessing SNP quality..what tools are used for that? I already see some sequencing errors when browsing the data.

As you can tell, I am totally new with this. It is pretty exciting so I want to learn and be able to do some of these things! Thanks in advance.

edit: I also get so confused with some of the output, more detailed documentation on that would be nice as well!

ADD COMMENTlink modified 4.8 years ago by alexisdereeper30 • written 9.7 years ago by Matthieu460

You referenced several other posts but not this one ; you may find it helpful.

ADD REPLYlink modified 15 months ago by _r_am32k • written 9.7 years ago by David Quigley11k

ah, I missed that. It is very helpful. Thank you, David!

ADD REPLYlink written 9.7 years ago by Matthieu460
gravatar for moorem
5.3 years ago by
United Kingdom
moorem230 wrote:

GATK is certainly the best variant caller and it's the difference is important because with indels in particular, they can be hugely different. 

I've spoken to a few people who have had issues with GATK so I've actually written a blog post with a workflow that works for GATK. The major issues seemed to be using GATK's best practices, particularly their add or replace read groups module which sometimes causes it to fall over later. By adding the read group (you can make it whatever you like) when you run bwa, then this workflow actually works. 

This was written for microbes so ploidy is set to 1 for HaplotypeCaller, but this can be set to whatever you like (do check the documentation). Originally posted here(



You’ll need to install picardtools, GATK, bwa and optionally, vcffilter for this workflow. Picardtools and GATK are simply .jar files so that’s no problem while you probably already have bwa installed, otherwise installation is well documented!

-The workflow

This workflow begins with short read (fastq) files and a fasta reference. First a sequence dictionary is created, short reads are aligned to the reference and read group information provided, resulting sequence alignment map (sam) file sorted and converted to binary alignment map (bam) format, duplicates marked, bam file sorted, indel targets identified, indels realigned and variants called. Simple!

For simplicity an example set of commands are provided here, where the reference is reference.fasta and the short reads are R1.fastq.gz and R2.fastq.gz. You will need to enter the paths and versions of the software being used at each step and your actual file names.

#Create sequence dictionary
java -jar~/bin/picard-tools-1.8.5/CreateSequenceDictionary.jar REFERENCE=reference.fasta OUTPUT=reference.dict

#Align reads and assign read group
bwa mem -R “@RG\tID:FLOWCELL1.LANE1\tPL:ILLUMINA\tLB:test\tSM:PA01” reference.fasta R1.fastq.gz R2.fastq.gz > aln.sam

#Sort sam file
java -jar ~/bin/picard-tools-1.8.5/SortSam.jar I=aln.sam O=sorted.bam SORT_ORDER=coordinate

#Mark duplicates
java -jar ~/bin/picard-tools-version/MarkDuplicates.jar I=sorted.bam O=dedup.bam METRICS_FILE=metrics.txt

#Sort bam file
java -jar ~/bin/picard-tools-version/BuildBamIndex.jar INPUT=dedup.bam

#Create realignment targets
java -jar ~/bin/GATK3.3/GenomeAnalysisTK.jar -T RealignerTargetCreator -R reference.fasta -I dedup.bam -o targetintervals.list

#Indel realignment
java -jar ~/bin/GATK3.3/GenomeAnalysisTK.jar -T IndelRealigner -R PA01.fasta -I dedup.bam -targetIntervals targetintervals.list -o realigned.bam

#Call variants (HaplotypeCaller)
java -jar ~/bin/GATK3.3/GenomeAnalysisTK.jar -T HaplotypeCaller -R reference.fasta -I realigned.bam -ploidy 1 -stand_call_conf 30 -stand_emit_conf 10 -o raw.vcf

The resulting vcf file will contain your variant calls!

Then you can optionally filter the variants:

#Filter variants
~/bin/vcflib/bin/vcffilter -f ‘DP > 9’ -f ‘QUAL > 10’ raw.vcf > filtered.vcf

Or first split the raw.vcf file into SNPs and indels:

#Extract SNPs
java -jar ~/bin/GATK3.3/GenomeAnalysisTK.jar -T SelectVariants -R reference.fasta -V raw.vcf -selectType SNP -o snps.vcf

#Extract Indels
java -jar ~/bin/GATK/GenomeAnalysisTK.jar -T SelectVariants -R reference.fasta -V raw.vcf -selectType INDEL -o indels.vcf

I also have a neat perl wrapper to automate this workflow over many short read files and would be happy to make this available if people are interested! Please do comment with any questions or issues and I’ll do my best to resolve them!

ADD COMMENTlink modified 5.3 years ago • written 5.3 years ago by moorem230

wow thanks for the workflow! I guess this will be my starting point in this domain. have you already automated this in a perl script as you have described ? and one more question, are there certain parameters that need to be considered when working with reads from human genome?

ADD REPLYlink written 4.7 years ago by tariktj0

Just to give you an idea, this can be automated in any language be it Perl/Python/bash as per your comfort, just to be sure about the Base Recalibration step, I do not see that here . Take a look at this link . The Base Recalibration is a standard best practice measure suggested. Another thing on the filtration is that one can discard low quality variants or even for that fact make a plot of various quality scores and then discard the first quantile or even lesser. You can also use the mentioned one in the GATK website but it depends upon the quality of your data and to what extent you had a coverage that can best suit those parameter. I would always suggest filtering of variants should be dependent to each one to its own data based on the distributions you have.

ADD REPLYlink written 4.7 years ago by ivivek_ngs5.0k


Can I also use this on RNA-seq data or does this only work on genome studies? is it okay to use GATK indel realignment on a RNA-seq data? Thanks!

ADD REPLYlink written 4.4 years ago by jfo0

Depends on what context you want to use this work flow on RNA-Seq, if you want to find variants from RNA-Seq data then GATK workflow can be used but it is better to use the recommended one which is here. If you can make your own workflow it is fine. Else you can refer to the automated script here I made. But you can always do the filtration based on your distributions as well. Feel free to ask.

ADD REPLYlink written 4.4 years ago by ivivek_ngs5.0k
gravatar for Pablo
9.7 years ago by
Pablo1.9k wrote:
  • BWA using defaults it's probably OK.

  • If you have a SAI file from the previous step, you need to convert it to SAM or BAM. Do something like (assuming your reference genome is hg37.fasta)

    bwa samse hg37.fasta s.sai s.fastq > s.sam

  • Then, Create BAM file (we assume you installed SamTools)

    samtools view -S -b s.sam > s.bam

  • Sort BAM file (will create s_sort.bam)

    samtools sort s.bam s_sort

  • Call variants: I.e. Create VCF file (BcfTools is part of samtools distribution)

    samtools mpileup -uf hg37.fasta s_sort.bam | bcftools view -vcg - > s.vcf

There is a lot more (like local realignment, etc.). But if this is your first time doing it, you should start with the basics.

ADD COMMENTlink written 9.7 years ago by Pablo1.9k
gravatar for Jim Vallandingham
9.7 years ago by
Kansas City, MO
Jim Vallandingham340 wrote:

We are working on a SNP pipeline now. You might find my work in progress pipeline useful.

The pipeline currently starts with an alignment from BWA. It uses GATK for SNP calling.

Briefly, the flow involves:

  • Realigning the BAM file using GATK's RealignerTargetCreator and IndelRealigner
  • Optional recalibration step (we currently don't use)
  • SNP calling using GATK's UnifiedGenotyper
  • Indel calling using the UnifiedGenotyper
  • basic filtering of the resulting VCF files
    • Right now we use some basic metrics to attempt to filer out low quality SNPs. I'm sure this step could be improved
  • Annotate called and filtered SNPs
    • Currently we use a custom script to add gene / transcript / exon / intron and other information from Ensembl.
    • Recently (yesterday) I found snpEff from another BioStar discussion. We will be using it in the future for this kind of annotation.

We are still fleshing out the details on filtering and such, but it might be a good starting point to for executing GATK in a working order

ADD COMMENTlink modified 8.8 years ago • written 9.7 years ago by Jim Vallandingham340

I've been asked to do something similar and would use the same programs as you, so I will have a look and see if I can get yours working. snpEff is very good for annotation as I expect you have found. Although I would look to expand and include a structural variation software such as pindel or Delly (embl uses this) because GATK does not detect large inserts well. Thanks very much for posting this, hopefully looking through your code will accelerate the process of putting my pipeline or modifying yours to what we want.

ADD REPLYlink written 7.4 years ago by rob234king600
gravatar for Khader Shameer
9.4 years ago by
Manhattan, NY
Khader Shameer18k wrote:

I strongly recommend this recent article from authors of GATK.

It covers various aspect associated with SNP calling in detail. At the same time do refer the software manual/wiki for up-to-date options incorporated in the toolkit.

alt text

ADD COMMENTlink written 9.4 years ago by Khader Shameer18k
gravatar for rob234king
7.5 years ago by
UK/Harpenden/Rothamsted Research
rob234king600 wrote:

I have put together a tutorial website with four core tutorials on it, RNA-Seq, ChIP-Seq, Genome assembly, and SNP calling that may be of use to you.

This website was created to share bioinformatics tutorials.

ADD COMMENTlink modified 5.8 years ago • written 7.5 years ago by rob234king600

Login - Register is a major put-off, maybe make it optional?

ADD REPLYlink written 7.5 years ago by zx87549.9k

Thanks for the input, I know registering is a bit of a pain but the website is based upon groups that contain tutorials and I use the login to only show those tutorials from groups you are a member of to make it manageable in the future if it got bigger. You can use false details if your concerned, email etc it's just used for this purpose. If you register it will automatically log you in and if go to "groups" on side bar and join the "Cranfield University" group, then click on "group tutorials" in left side bar you can access the tutorials I have put up. I think they need a little re-formatting but have some useful demonstrations of tools like RSATs etc. Thanks again for the reply appreciate it.

ADD REPLYlink written 7.5 years ago by rob234king600

You got me curious, but the login system is not working…

ADD REPLYlink written 7.4 years ago by Andre Elias80

Thanks I'll check it out it seems the server needs restarting, there is a memory leak on one or more of the deployments on this server. I ask for it to be reset and investigate.

ADD REPLYlink written 7.4 years ago by rob234king600

A number of student are submitting to this server over the next month which is resulting in permgen errors, some kind of memory leak somewhere or due to multiple submissions. Issue corrected for the moment.

CUBELP2 is a static website and CUBELP a sharing platform, typing their names and "bioinformatics" which should locate them in google.

ADD REPLYlink modified 7.4 years ago • written 7.4 years ago by rob234king600

Please fix the link. It doesn't work anymore!

ADD REPLYlink written 5.8 years ago by SmallChess540

ADD REPLYlink written 5.8 years ago by rob234king600
gravatar for Anil
9.4 years ago by
Anil10 wrote:

Use GALAXY software. It is free and user friendly and u will get most of the step in oe programme only.

ADD COMMENTlink written 9.4 years ago by Anil10
gravatar for Ming Tang
7.5 years ago by
Ming Tang2.6k
Houston/MD Anderson Cancer Center
Ming Tang2.6k wrote:

have a look at this course material from UT-Austin

ADD COMMENTlink written 7.5 years ago by Ming Tang2.6k
gravatar for alexisdereeper
4.8 years ago by
alexisdereeper30 wrote:

You might want to filter and explore your VCF file for different kinds of analyses (diversity, haplotypes, population structure). You can then make use of the online SNP pipeline:

ADD COMMENTlink written 4.8 years ago by alexisdereeper30

I tried to use a vcf 4.2 file but does not give any kind of output, the images are all blank. The file size was also around 21 MB so that kinda find since I ran the example of the polymorphism vcf file which seemed to be ok. Does it work with new vcf4.2 format? Would like to try on my variant file which I generated for one sample from RNA-Seq and should be having a lot of false positives as there was no matched normal , still would like to see the various statistics, so was trying the tool but it does not produce anything.

ADD REPLYlink written 4.8 years ago by ivivek_ngs5.0k

It should work with vcf4.2 format. Maybe there was a special caracter in the sample names that blocked the VCFtools program...

ADD REPLYlink written 4.7 years ago by alexisdereeper30
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1074 users visited in the last hour