Open main menu

Genome Analysis Wiki β

SeqShop: Sequence Mapping and Assembly Practical, December 2014

Contents

Introduction

Main Workshop wiki page: SeqShop: December 2014

See the introductory slides for an intro to this tutorial.

Goals of This Session

  • What we want to learn
    • Basic sequence data file formats (FASTQ, BAM)
    •  How to generate aligned sequences that are ready for variant calling from raw sequence reads
    •  How to evaluate the quality of sequence data
    •  How to visualize sequence data to examine the reads aligned to particular genomic positions

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 will setup some environment variables to point you to

    • GotCloud program
    • Tutorial input files
    • Setup an output directory
    source /net/seqshop-server/home/mktrost/seqshop/setup.txt
    
    • You won't see any output after running source
      • It silently sets up your environment

    Look at setup.txt

    cat /net/seqshop-server/home/mktrost/seqshop/setup.txt
    
    • setup.txt screenshot

    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

    Download and untar file containing the example data used in the practicals:

    mkdir -p ~/seqshop
    cd ~/seqshop
    wget http://csg.sph.umich.edu/mktrost/seqshopExampleDec2014.tar.gz
    tar xvf seqshopExampleDec2014.tar.gz
    

    You will see the names of all the files included in the example data scrolling on the screen as they are unpacked from the tar file.

    Download & Build GotCloud

    If you do not already have GotCloud:

    • download, decompress, and build the version of gotcloud that was tested with this tutorial:
    wget https://github.com/statgen/gotcloud/archive/gotcloud.1.15.tar.gz
    tar xvf gotcloud.1.15.tar.gz
    mv gotcloud-gotcloud.1.15 gotcloud
    cd gotcloud/src
    make
    cd ../..
    

    Remember the path to gotcloud/ that is what you will need to set your GC variable to.

    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.

    1. Point to where you installed GotCloud
    2. Point to where you installed the seqshop files
    3. Point to where you want the output to go
    Using bash (replace the paths below with the appropriate paths):
    export GC=~/seqshop/gotcloud
    export SS=~/seqshop/example
    export OUT=~/seqshop/output
    Using tcsh (replace the paths below with the appropriate paths):
    setenv GC ~/seqshop/gotcloud
    setenv SS ~/seqshop/example
    setenv OUT ~/seqshop/output

    Examining GotCloud Align Input Files

    Examining Raw Sequence Reads : FASTQs

    FASTQ : standard file format provided to you by those who did the sequencing.

    For more information on the FASTQ format, see: http://en.wikipedia.org/wiki/FASTQ_format

    For this tutorial, we will use FASTQs for 4 1000 Genome samples

    • Subset of FASTQs - should map to chromosome 22 36000000-37000000
    ls ${SS}/fastq/
    

    There are 24 fastq files: combination of single-end & paired-end.

    Can you tell which files are single-end and which are paired-end?
    • View answer:
      • Paired-end files have a _1.fastq or _2.fastq extension
      • This convention isn't mandatory, but something similar is common
      • Fastqsm.png


    Look at a couple of FASTQs:

    less -S ${SS}/fastq/HG00551.SRR190851_1.fastq
    

    less is a Linux command that allows you to look at a file.

    • -S option prevents line wrap
    • Use the arrow (up/down/left/right) keys to scroll through the file
    • Use the space bar to jump down a page

    Use 'q' to exit out of less

    q
    
    Do you remember the parts of a FASTQ?
    • No, remind me:
    • Fastq.png


    Look at the paired read:

    less -S ${SS}/fastq/HG00551.SRR190851_2.fastq 
    

    Remember, use 'q' to exit out of less

    q
    
    Do you notice something in common?
    • View answer:
      • Paired-end reads have matching read names with a different extensions
      • This convention isn't mandatory, but something similar is common
      • Fastq3.png

    Reference Files

    Reference files can be downloaded with GotCloud or from other sources

    For alignment, you need:

    1. Reference genome FASTA file
    2. VCF (variant call format) files with chromosomes/positions
      • dbsnp - used to skip known variants when recalibrating
      • hapmap - used for sample contamination/sample swap validation

    Take a look at the chromosome 22 reference files included for this tutorial:

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

    Let's read the reference FASTA file (all reference bases for the chromosome):

    less ${SS}/ref22/human.g1k.v37.chr22.fa
    

    Remember, use 'q' to exit out of less

    q
    
    Where is the reference sequence?
    • Answer:
      • The ends of a chromosome are 'N' - unknown bases
      • Let's look at 5 lines of the file starting at line 300,000
       tail -n+300000 ${SS}/ref22/human.g1k.v37.chr22.fa |head -n 5
      

      Fasta.png

    If you want to access the FASTA file by position, you can use samtools faidx command

    ${GC}/bin/samtools faidx ${SS}/ref22/human.g1k.v37.chr22.fa 22:36000000 | less
    

    or

    ${GC}/bin/samtools faidx ${SS}/ref22/human.g1k.v37.chr22.fa 22:36000000-36000100
    

    GotCloud FASTQ List File

    The FASTQ list file is created by you to tell GotCloud about each of your FASTQ files:

    • Where to find it
    • Sample name
      • Each sample can have multiple FASTQs
      • Each FASTQ is for a single sample
    • Run identifier (optional)
      • For recalibration we need to know which reads were in the same run.

    FASTQ List Format:

    • Tab delimited
    • Starts with a header line
    • One line per single-end read
    • One line per paired-end read (only 1 line per pair).

    Let's look a look at the FASTQ list file I prepared for this tutorial:

    less -S ${SS}/fastq.list 
    

    Remember, use 'q' to exit out of less

    q
    
    Which samples have multiple paired end reads?
    • Need a reminder of the format?
    • FqindexNew.png

      • Note: in the screenshots, the fields are shifted into clear columns to make it easier to read
        • When you view the file, the fields will not line up in neat columns and it can be hard to read
    • Hard to read the index? Need a hint?
      • Use cut to extract just the SAMPLE & FASTQ2 fields
      cut -f 1,3 ${SS}/fastq.list
      
    • Answer:
      • HG00553 & HG00640
      • They have multiple FASTQ2 files listed
      • FqListFASTQ2.png


    How do you point GotCloud to your FASTQ list file?

    • Command-line --list option
    or
    • Configuration file FASTQ_LIST setting.

    The command-line setting takes precedence over the configuration file setting.

    GotCloud Configuration File

    This file is created by you to configure GotCloud for your data.

    • Default values are provided in ${GC}/bin/gotcloudDefaults.conf
      • Most values should be left as the defaults
    • Specify values in your configuration file as:
      • KEY = value
    • Use $(KEY) to refer to another key's value
    • If a KEY is specified twice, the later value is used
    • Does not have access to environment variables
    • '#' indicates a comment

    Let's look at the configuration file I created for this test:

    more ${SS}/gotcloud.conf
    

    Use the space bar to advance if the whole file isn't displayed.

    If your references are in a different path than what is specified, what would you change?
    • Answer:
      • You would change REF_DIR to the new path
      • GcConfNew.png

    Run GotCloud Align

    AlignDiagram.png

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

    • When running at home if you don't have 4 CPUs, reduce the --numjobs setting (it will take longer to run).
    ${GC}/gotcloud align --conf ${SS}/gotcloud.conf --numjobs 4 --base_prefix ${SS} --outdir ${OUT}
    
    • ${GC}/gotcloud runs GotCloud
    • align tells GotCloud you want to run the alignment pipeline.
    • --conf tells GotCloud the name of the configuration file to use.
      • The configuration for this test was downloaded with the seqshop input files.
    • --numjobs means to run 4 samples at a time.
      • How many you can run concurrently depends on your system.
    • --base_prefix tells GotCloud the prefix to append to relative paths.
      • The Configuration file cannot read environment variables, so we need to tell GotCloud the path to the input files, ${SS}
      • Alternatively, gotcloud.conf could be updated to specify the full paths
    • --outdir tells GotCloud where to write the output.
      • This could be specified in gotcloud.conf, but to allow you to use the ${OUT} to change the output location, it is specified on the command-line

    GcalignStartNew.png

    This should take about 1 minute to run.

    It should end with a line like: Processing finished in 54 secs with no errors reported

    • The WARNING messages are just to let you know that the default Read Group field settings are being used.


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


    Inside GotCloud align, not only sequence alignment but also pre-processing of sequence data, including deduplication and base quality recalibration are performed along with quality assessment, as illustrated below.

    Gotcloud align detail.png

    Examining GotCloud Align Output

    Let's look at the output directory:

    ls ${OUT}
    

    GcalignOutMNew.png

    Quality Control Files

    Let's take a look at our quality control output directory:

    ls ${OUT}/QCFiles 
    

    GcalignOutQCm.png

    Sample Contamination/Swap

    Check for sample contamination:

    • *.selfSM : Main output file containing the contamination estimate.
      • Check the 'FREEMIX' column for genotype-free estimate of contamination
          • 0-1 scale, the lower, the better
          • If [FREEMIX] >= 0.03 and [FREELK1]-[FREELK0] is large, possible contamination
      • See VerifyBamID: A guideline to interpret output files for more information
    less -S ${OUT}/QCFiles/HG00551.genoCheck.selfSM
    

    Remember, use 'q' to exit out of less

    q
    
    Is there evidence of sample contamination?
    • Answer:
      • No, FREEMIX = 0.00000 (<0.03)

      Contam1New.png

    QC Metrics

    See: QPLOT: Diagnose sequencing quality for more info on how to use QPLOT results.

    Let's look at some quality control metrics:

    cat ${OUT}/QCFiles/HG00551.qplot.stats
    
    What is the mapping rate & average coverage for HG00551?
    • Answer
      • 98.93% Mapped
      • 7.44 MeanDepth

      QplotsNew.png

    Generate a pdf of quality metrics:

    Rscript ${OUT}/QCFiles/HG00551.qplot.R
    

    Examine the PDF:

    evince  ${OUT}/QCFiles/HG00551.qplot.pdf&
    

    It is ok if you see a warning message when opening evince. It should still open. If not, let me know. To close evince, just close the pdf window.

    Does the Empirical vs reported Phred score look as good as we would like?
    • Answer
      • No, it is well above the line
      • This is due to the small region used for recalibration
      • QplotpdfNew.png

      • Look at the PDF I produced when I ran the whole genome:
      evince ${SS}/ext/HG00551.wg.qplot.pdf&
      

      Qplotpdfwg.png

    BAM Files

    Binary Sequence Alignment/Map (SAM) Format

    Let's look at the BAMs (aligned reads that are ready for variant calling):

    ls ${OUT}/bams
    

    GcalignOutBAMm.png

    Let's examine at the first 7 lines of the BAM file using samtools view:

    ${GC}/bin/samtools view -h ${OUT}/bams/HG00551.recal.bam|head -n 7
    
    What are the chromosome and position of the first record in the BAM file?
    • Need a reminder of the format?
    • Bam.png

    • Answer
      • Chr 22, Pos: 16918656

      BamRecNew.png

    Accessing BAMs by Position

    BAM's are so big, what if we want to see a position part way through the file?

    Add a region to the view command we used above. Let's find all reads that overlap positions 36907000-36907005:

    ${GC}/bin/samtools view -h ${OUT}/bams/HG00551.recal.bam 22:36907000-36907005
    
    • Just a few reads.

    Let's visualize what reads in that area look like using samtools tview:

    ${GC}/bin/samtools tview ${OUT}/bams/HG00551.recal.bam ${SS}/ref22/human.g1k.v37.chr22.fa
    
    • Type ‘g’
      • Type 22:36907000
    • Type ‘n’ to color by nucleotide
    • Use the arrow keys to move around and look at the area.

    Understanding the syntax:

    • '.' : match to the reference on the forward strand
    • ',' : match to the reference on the reverse strand
    • ACGTN : mismatch to reference on the forward strand
    • acgtn : mismatch to reference on the reverse strand
    Do you see anything interesting?
    • Screenshot
      • We will have to remember this region when we run snpcall to see what it says.

      TviewNew.png

    Other tview commands:

    • Type '?' for a help screen
    • Type 'q' to quit tview

    Feel free to play around more and browse the BAM files.

    Other tools for BAMs

    We have developed a lot of tools that operate on BAM files.

    See Software: BamUtil Tools for a list

    • Many operations:
      • diff : diff 2 BAM files
      • stats: per positions statistics
      • bam2Fastq : convert a BAM back to a FASTQ (how I created the fastqs for this tutorial)
      • Lots of others
    • Feel free to try some out
    • If you have any questions, let me know, I wrote most of them and am happy to help.

    Logging Off

    This section is specifically for the SeqShop Workshop computers.

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

    To logout of seqshop1/2/3/4, type:

    exit
    

    To logout of seqshop-server, type:

    exit
    

    And close the windows.

    When done, log out of the Windows machine.

    Return to Workshop Wiki Page

    Return to main workshop wiki page: SeqShop: December 2014