Open main menu

Genome Analysis Wiki β

SeqShop: Variant Calling and Filtering for SNPs Practical, June 2014

Revision as of 11:07, 26 June 2014 by Mktrost (talk | contribs)

Contents

Introduction

See the introductory slides for an intro to this tutorial.


Goals of This Session

  • What we want to learn
    •  How to generate filtered variant calls for SNPs from BAMs
    • Basic variant call file format (VCF)
    •  How to examine the variants at particular genomic positions
    •  How to evaluate the quality of SNP calls

Setup in person at the SeqShop Workshop

This section is specifically for the SeqShop Workshop computers.

If you are not running during the SeqShop Workshop, please skip this section.


Login to the seqshop-server Linux Machine

This section will appear redundantly in each session. If you are already logged in or know how to log in to the server, please skip this section

  1. Login to the windows machine
  • The username/password for the Windows machine should be written on the right-hand monitor
  • Start xming so you can open external windows on our Linux machine
    • Start->Enter "Xming" in the search and select "Xming" from the program list
    • Nothing will happen, but Xming was started.
    • View Screenshot
    • Xming.png

  • Open putty
    • Start->Enter "putty" in the search and select "PuTTY" from the program list
    • View Screenshot
    • PuttyS.png

  • Configure PuTTY in the PuTTY Configuration window
    • Host Name: seqshop-server.sph.umich.edu
    • View Screenshot
    • Seqshop.png

    • Setup to allow you to open external windows:
      • In the left pannel: Connection->SSH->X11
        • Add a check mark in the box next to Enable X11 forwarding
        • View Screenshot
        • SeqshopX11.png

    • Click Open
    • If it prompts about a key, click OK
  • Enter your provided username & password as provided

  • You should now be logged into a terminal on the seqshop-server and be able to access the test files.

    • If you need another terminal, repeat from step 3.

    Login to the seqshop Machine

    So you can each run multiple jobs at once, we will have you run on 4 different machines within our seqshop setup.

    • You can only access these machines after logging onto seqshop-server

    3 users logon to:

    ssh -X seqshop1
    

    3 users logon to:

    ssh -X seqshop2
    

    2 users logon to:

    ssh -X seqshop3
    

    2 users logon to:

    ssh -X seqshop4
    

    Setup your run environment

    This is the same setup you did for the previous tutorial, but you need to redo it each time you log in.

    This will setup some environment variables to point you to

    • GotCloud program
    • Tutorial input files
    • Setup an output directory
      • It will leave your output directory from the previous tutorial in tact.
    source /home/mktrost/seqshop/setup.txt
    
    • You won't see any output after running source
      • It silently sets up your environment
      • If you want to view the detail of the setup, type
    less /home/mktrost/seqshop/setup.txt
    

    and press 'q' to finish.

    View setup.txt

    Setup.png

    Setup when running on your own outside of the SeqShop Workshop

    This section is specifically for running on your own outside of the SeqShop Workshop.

    If you are running during the SeqShop Workshop, please skip this section.

    Download the example data

    Setup your run environment

    Environment variables will be used throughout the tutorial.

    We recommend that you setup these variables so you won't have to modify every command in the tutorial.

    I'm using bash (replace the paths below with the appropriate paths):

    • Point to where you installed GotCloud
      export GC=/home/username/gotcloud
    • Point to where you installed the seqshop files
      export SS=/home/username/seqshop/
    • Point to where you want the output to go
      export OUT=/home/username/seqshop_output/

    I'm using tcsh (replace the paths below with the appropriate paths):

    • Point to where you installed GotCloud
      setenv GC /home/username/gotcloud
    • Point to where you installed the seqshop files
      setenv SS /home/username/seqshop/
    • Point to where you want the output to go
      setenv OUT /home/username/seqshop_output/


    Examining GotCloud SnpCall Input files

    Sequnce Alignment Files: BAM Files

    Per sample BAM files contain sequence reads that are mapped to positions in the genome.

    For a reminder on how to look at/read BAM files, see: SeqShop Aligment: BAM Files

    For this tutorial, we will use the 4 BAMs produced in the SeqShop: Sequence Mapping and Assembly Practical as well as with 58 BAMs that were pre-aligned to that 1MB region of chromosome 22.

    Reference Files

    Reference files can be downloaded with GotCloud or from other sources.

    For GotCloud snpcall, you need:

    1. Reference genome FASTA file
      • Contains the reference base for each position of each chromosome
        • Used to compare bases in sequence reads to the reference positions they mapped to
        • Used to identify SNPs/variations in the sequence reads
      • Additional information on the FASTA format: http://en.wikipedia.org/wiki/FASTA_format
    2. VCF (variant call format) files with chromosomes/positions
      • indel - contains known insertions & deletions to help with filtering
      • omni - used as likely true positives for SVM filtering
      • hapmap - used as likely true positives for SVM filtering and for generating summary statistics
      • dbsnp - used for generating summary statistics


    We looked at them yesterday, but you can take another look at the chromosome 22 reference files included for this tutorial:

    ls ${SS}/ref22
    
    • View Screenshot
    • RefDir.png

    GotCloud BAM Index File

    The BAM index file points GotCloud to the BAM files

    • generated by the alignment pipeline

    Look at the BAM index file the alignment pipeline generated

    cat ${OUT}/bam.index
    
    What is the path to the BAM file for sample HG00640?
    • Answer:
      • /home/YourUserName/out/bams/HG00640.recal.bam
      • Bamindex.png

    The alignment pipeline only processed 4 samples, but for snpcall, we want to run on 62 samples.

    • The other 58 samples were already aligned:
    ls ${SS}/bams
    

    Look at the BAM index for those BAMs:

    less ${SS}/bams/bam.index
    

    Remember, use 'q' to exit out of less

    q
    
    Do you notice a difference between this index and yours?
    • Answer:
      • It doesn't have a full path to the BAM file, while your index has /home/...
      • Bamindex1.png

      • That's ok, gotcloud.conf contains the path to those BAMs
      • BamindexConf.png


    We need to add these BAMs to our index

    • Append the bam.index from the pre-aligned BAMs to the one you generated from the alignment pipeline
    cat ${SS}/bams/bam.index >> ${OUT}/bam.index
    
    • ">>" will append to the file that follows it
    • Be sure to do this command just once
      • Check that your BAM index is the correct size
        wc -l ${OUT}/bam.index
        • wc -l counts the number of lines in the file
        • Should be 62

    Verify your BAM index contains the additional BAMs

     less ${OUT}/bam.index
    

    Remember, use 'q' to exit out of less

    q
    
    Do you see both sets of BAMs?
    • Annotated Screenshot:
      • If not, let me know
      • Bamindex2.png

    GotCloud Configuration File

    We will use the same configuration file as we used yesterday in GotCloud Align.

    See SeqShop: Alignment: GotCloud Configuration File for more details

    • Note we want to limit snpcall to just chr22 so the configuration already has CHRS = 22 (default was 1-22 & X).

    For more information on configuration, see: GotCloud snpcall: Configuration File

    • Contains information on how to configure for exome/targeted sequencing

    Run GotCloud SnpCall

    SnpcallDiagram.png

    Now that we have all of our input files, we need just a simple command to run:

    ${GC}/gotcloud snpcall --conf ${SS}/gotcloud.conf --numjobs 4 --region 22:36000000-37000000 --base_prefix ${SS} --outdir ${OUT}
    
    • --numjobs tells GotCloud how many jobs to run in parallel
      • Depends on your system
    • --region 22:36000000-37000000
      • The sample files are just a small region of chromosome 22, so to save time, we tell GotCloud to ignore the other regions

    Curious if it started running properly? Check out this screenshot:

    SnpcallStart.png

    This should take about 5 minutes to run.

    • After about 4 minutes of running, GotCloud snpcall will output some text to the screen. Don't worry, that is expected and is just output from some of the intermediate tools.
    • It should end with a line like: Commands finished in 289 secs with no errors reported

    If you cancelled GotCloud part way through, just rerun your GotCloud command and it will pick up where it left off.

    If you want to understand more detailed step of GotCloud SNP calling, here is a schematic picture with a little bit more details

    Gotcloudoverview.png

    Examining GotCloud SnpCall Output

    Let's look at the output directory:

    ls ${OUT}
    
    Do you see any new files or directories?
    • View Annotated Screenshot:

    GcsnpcallOut.png

    Let's look at the vcfs directory:

    ls ${OUT}/vcfs
    

    Just a chr22 directory, so look inside of there:

    ls ${OUT}/vcfs/chr22
    
    Can you identify the final filtered VCF and the associated summary file?
    • Answer & annotated directory listing:
    • Filtered VCF (SVM & hard filters): chr22.filtered.vcf.gz
    • Summary file: chr22.filtered.sites.vcf.summary

    Vcfsout.png

    Now, let's look in the split directory for the VCF with just the passing variants:

    ls ${OUT}/split/chr22
    
    Which file do you think is the one you want?
    • Answer:
    • chr22.filtered.PASS.vcf.gz

    SplitOut.png

    Filtering Summary Statistics

    cat ${OUT}/vcfs/chr22/chr22.filtered.sites.vcf.summary
    

    View Screenshot

    FilterSum.png

    To understand how to interpret the filtering summary statistics, please refer to Understanding vcf-summary output

    Filtered VCF

    Let's look at the filtered sites file.

    less -S ${OUT}/vcfs/chr22/chr22.filtered.sites.vcf
    
    • Scroll down until you find some variants.
      • Use space bar to jump a full page
      • Use down arrow to move down one line
    • Scroll right: lots of info fields, but no per sample genotype information
    What is the first filtered out variant that you find & what filter did it fail?
    • Answer:
    • It failed SVM filter

    SvmFilt.png

    Remember, use 'q' to exit out of less

    q
    


    Now, let's look at the filtered file with genotypes.

    zless -S ${OUT}/vcfs/chr22/chr22.filtered.vcf.gz
    
    • Scroll down until you find some variants.
      • Use space bar to jump a full page
      • Use down arrow to move down one line
    • Scroll right until you should see per sample genotype information

    Remember, use 'q' to exit out of less

    q
    
    • View annotated screenshot:

    SvmFiltGL.png

    Passing SNPs

    Let's look at the file of just the pass sites:

    zless -S ${OUT}/split/chr22/chr22.filtered.PASS.vcf.gz
    
    • Scroll down: they all look like they PASS

    Remember, use 'q' to exit out of less

    q
    

    Let's check if they are all PASS.

    zcat ${OUT}/split/chr22/chr22.filtered.PASS.vcf.gz |grep -v "^#"| cut -f 7| grep -v "PASS"
    

    It will return nothing since there are no non-passing variants in this file.

    Want an explanation of this command?
    • zcat ...: uncompress the zipped VCF
    • '|' : this takes the output of one command and sends it as input to the next
    • grep -v "^#" : exclude any lines that start with "#" - headers
    • cut -f 7 : extract the FILTER column (the 7th column)
    • grep -v "PASS" : exclude any rows that have a "PASS" in the FILTER column

    Compare that to the filtered file we looked at before:

    zcat ${OUT}/vcfs/chr22/chr22.filtered.vcf.gz |grep -v "^#"| cut -f 7| grep -v "PASS"
    
    Do you see any filters?
    • Answer
    • Yes
      • It should have scrolled and you should see filters like:
        • INDEL5;SVM
        • INDEL5
        • SVM

    GotCloud Genotype Refinement

    To improve the quality of the genotypes, we run a genotype refinement pipeline.

    This pipeline runs Beagle & thunder.

    Genotype Refinement Input

    The GotCloud genotype refinement pipeline takes as input ${OUT}/split/chr22/chr22.filtered.PASS.vcf.gz (the VCF file of PASS'ing SNPs from snpcall).

    The bam index and the configuration file we used for GotCloud snpcall will tell GotCloud genotype refinement everything it needs to know, so no new input files need to be prepared.

    Note: the configuration file overrides the THUNDER command to make it go faster than the default settings so the tutorial will run faster: ThunderConf.png

    Running GotCloud Genotype Refinement

    Since everything is setup, just run the following command (very similar to snpcall).

    ${GC}/gotcloud ldrefine --conf ${SS}/gotcloud.conf --numjobs 2 --region 22:36000000-37000000 --base_prefix ${SS} --outdir ${OUT}
    
    • Beagle will take about 2-3 minutes to complete
    • Thunder will automatically run and will take another 3-4 minutes

    Genotype Refinement Output

    What's new in the output directory?
    • Answer
      • beagle directory : Beagle output
      • thunder directory : Thunder output
      • umake.beagle.conf : Configuration values used for GotCloud beagle
      • umake.beagle.Makefile : GNU makefile for commands run as part of GotCloud beagle
      • umake.beagle.Makefile.log : Log of the GotCloud beagle run
      • umake.thunder.* files : Contain the configuration & steps used in GotCloud thunder

    Let's take a look at that interesting location we found in the alignment tutorial : chromosome 22, positions 36907000-36907100

    Use tabix to extract that from the VCFs:

    ${GC}/bin/tabix ${OUT}/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz 22:36907000-36907100 |less -S
    

    Remember, type 'q' to quit less.

    q
    
    Are there any variants in this region?
    • Answer:
      • Yes!
      • Positions:
        • 36907001; Ref: T, Alt: C - that's what we saw before
        • 36907098; Ref: T, Alt: C - that's what we saw before
    What is HG00551's genotype at these positions?
    1. First check which sample number HG00551 is:
    zcat ${OUT}/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz |grep "#CHROM"
    
    • That will help you figure out it's genotype.
    • Rerun the tabix command and scroll to find HG00551's genotype:
    ${GC}/bin/tabix ${OUT}/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz 22:36907000-36907100 |less -S
    
    • Answer:
      • It is the first sample
      • 0|1: Heterozygous
      • 1|1; Homozygous Alt (C)

    Remember, type 'q' to quit less.

    q
    

    Did I find interesting variants?

    The region we selected contains APOL1 gene, which is known to play an important role in kidney diseases such as nephrotic syndrome. One of the non-synonymous risk allele, rs73885139 located at position 22:36661906 increases the risk of nephrotic syndrome by >2-folds. Let's see if we found the interesting variant by looking at the VCF file by position.

    ${GC}/bin/tabix ${OUT}/vcfs/chr22/chr22.filtered.vcf.gz 22:36661906 | head -1 
    

    Did you see a variant at the position?

    ${GC}/bin/tabix ${OUT}/vcfs/chr22/chr22.filtered.vcf.gz 22:36661906 | head -1 
    22	36661906	.	A	G	18	PASS DP=409;MQ=59;NS=62;AN=124;AC=2;AF=0.013827;AB=0.4065;AZ=-0.5287;
                        FIC=-0.0092;SLRT=-0.0075;HWEAF=0.0138;HWDAF=0.0276,0.0000;LBS=36,36,0,0,1,1,0,0;
                       OBS=145,191,0,0,3,2,0,0;STR=-0.040;STZ=-0.740;CBR=0.008;CBZ=0.144;IOR=0.000;IOZ=-1.370;
                       AOI=-5.614;AOZ=-4.243;LQR=0.178;MQ0=0.000;MQ10=0.000;MQ20=0.000;MQ30=0.000;SVM=1.51214	
                  GT:DP:GQ:PL	0/0:4:28:0,12,65
    

    Let's check the sequence data to confirm that the variant really exists

    ${GC}/bin/samtools tview ${SS}/bams/HG01242.recal.bam ${SS}/ref22/human.g1k.v37.chr22.fa
    
    • Type 'g' to go to a specific position
    • Type 22:36661906 to move to the position
    • Press arrows to move between positions
    • Press 'b' if you want to color by base quality
    • Press '?' for more help

    View Screenshot

    Samtoolstviewsnp.png

    Improvements

    Let's get some information on the BEAGLE VCF:


    perl ${EXT}/bed-diff.pl --vcf1 ${REF}/1kg.omni.chr22.36Mb.vcf.gz --vcf2 ${OUT}/beagle/chr22/chr22.filtered.PASS.beagled.ALL.vcf.gz --gcRoot ${GC} --out ${OUT}/bedDiff.beagle
    

    Look at the results:

    more ${OUT}/bedDiff.beagle.summary
    
    • Results
    OVERALL:	43601	44293	0.9844
    NREF-EITHER:	19667	20359	0.9660
    NMAJ-EITHER:	14585	15277	0.9547
    
    HOMREF:	23934	100	1	0.9958
    HET:	329	11959	175	0.9596
    HOMALT:	4	83	7708	0.9888
    
    HOMMAJ:	29016	126	2	0.9956
    HET:	364	11959	140	0.9596
    HOMMIN:	3	57	2626	0.9777
    

    Now, let's see if it improved after running Thunder VCF:

    perl ${EXT}/bed-diff.pl --vcf1 ${REF}/1kg.omni.chr22.36Mb.vcf.gz --vcf2 ${OUT}/thunder/chr22/ALL/thunder/chr22.filtered.PASS.beagled.ALL.thunder.vcf.gz --gcRoot ${GC} --out ${OUT}/bedDiff.thunder
    

    Look at the results:

    more ${OUT}/bedDiff.thunder.summary
    
    • Results
    OVERALL:	43685	44293	0.9863
    NREF-EITHER:	19758	20366	0.9701
    NMAJ-EITHER:	14688	15296	0.9603
    
    HOMREF:	23927	106	2	0.9955
    HET:	286	12057	120	0.9674
    HOMALT:	6	88	7701	0.9879
    
    HOMMAJ:	28997	144	3	0.9950
    HET:	286	12057	120	0.9674
    HOMMIN:	5	50	2631	0.9795
    

    There is an improvement.

    What is GotCloud snpcall doing?

    To run GotCloud, you really just needed a single command.

    Well, that one command runs many steps. Here is a diagram of all the steps.

    SnpcallPipeline.jpg

    Aren't you glad you didn't have to configure & run each one yourself?