SeqShop: Sequence Mapping and Assembly Practical, June 2014

From Genome Analysis Wiki
Jump to navigationJump to search

Introduction

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


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
  2. 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

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

      • PuttyS.png

  4. 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
  5. 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 /home/mktrost/seqshop/setup.txt

View setup.txt

Setup.png


Alternatively, if you would like to change the output directory, copy the file, make the modifications and source your own file: 

cp /home/mktrost/seqshop/setup.txt ~/setup.txt
nedit ~/setup.txt
source ~/setup.txt

(You can use your favorite editor instead of nedit. I typically use emacs, but nedit is more like Windows.)

Examining Raw Sequence Reads

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 ${IN}/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 ${IN}/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 zless if the file is compressed.

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 ${IN}/fastq/HG00551.SRR190851_2.fastq
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

    GotCloud Alignment Pipeline

    AlignDiagram.png

    Why GotCloud?

    • Easy to learn & run
      • All-in-one sequence analysis pipeline
      • You don’t have to know the details of individual component
    • Robust parallelization
      •  Automatic partition of multi-sample jobs
      •  Reliable and fault-tolerant parallelization via GNU make
        • Restart from where it stopped upon unexpected crash
    • Cloud & Cluster-friendly
      • Supports multiple clusters such as MOSIX, Slurm, & SGE
      • Amazon instances allow running large-scale jobs without having your own cluster

    Sequence Processing Recommendations

    • Be consistent within a project
      • Process all samples with same pipeline
        • Batch effect may arise if different pipelines are used due to pipeline differences
      •  Use the same configuration within a project

    Examining GotCloud Align Input Files

    Sequence Data Files : FASTQs

    We already looked at those in: Examining Raw Sequence Reads

    Reference Files

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

    • Partial reference for this example
    ls ${GC}/reference/chr7

    RefDir.png


    VCF files

    • List of chromosome/position
    • Used for:
      • dbsnp - recalibration skips known variants
      • hapmap - used for sample contamination/sample swap validation
      • variant filtering


    Let's read the first 5 lines of the genome reference FASTA file (all reference bases for a chromosome): 

    head -n 5 ${GC}/reference/chr7/human.g1k.v37.chr7.fa

    Fasta3.png

    The start of the chromosome is all N's, so let's look at a later section (reading 5 lines starting at line 2000): 

    tail -n+2000 ${GC}/reference/chr7/human.g1k.v37.chr7.fa |head -n 5

    Fasta3 (copy).png

    Additional information on the FASTA format: http://en.wikipedia.org/wiki/FASTA_format

    See GotCloud: Genetic Reference and Resource Files for more information on downloading/generating reference files.

    GotCloud FASTQ Index File

    You need to tell GotCloud about each FASTQ file

    • Full path
    • Sample name
      • Each sample can have multiple FASTQs
      • Each FASTQ is for a single sample

    The FASTQ index file is created by you to direct GotCloud to your FASTQ files, providing additional information for them.

    • tab delimited
    • columns may be in any order
    • starts with a header line
    • one line per single-end read
    • one line per paired-end read (only 1 line per pair).


    Required Columns

    Column Name Description Recommended Value
    MERGE_NAME
    • Base name for the resulting BAM file for the sample
    • Used to group multiple fastqs or fastq pairs into a single BAM
    		 Sample Name
    FASTQ1
    • Name of the fastq or the first in the pair if paired-end. (Only 1 line per pair)
    		 path/fastq1
    FASTQ2
    • Name of the 2nd fastq in paired-end reads.
    • Column is not required if all fastqs are single-end
    • '.' if the column is used, but this line is single-ended.
    		 path/fastq2


    The following columns are optional and used to populate the Read Group Information in the BAM file.

    • RGID field is required if using any of these fields, the others are optional.

    What is a Read Group?

    • Groups reads together
    • Used for recalibration
      • Each sequencing run should get a different ReadGroup
    • Typically a new name for each fastq pair/group

    If you do not want the field for:

    • any fastq, leave the column out of the header line
    • a single line, use a '.'


    Optional Columns

    Column Name Description Recommended Value
    RGID Read Group ID Run ID
    SAMPLE Sample Name Sample Name
    LIBRARY Library
    • separate FASTQs for a sample that were prepped separately
    		 if you don't know or it is all the same, use Sample Name
    CENTER Center Name Name of the sequencing center producing the FASTQ
    PLATFORM Platform CAPILLARY, LS454, ILLUMINA,

    SOLID, HELICOS, IONTORRENT, or PACBIO

    Your sequencing core may provide to you a file with information to fill in these columns.

    For our example, we have sequence.index which contains the information from 1000 Genomes for the FASTQs we are processing. 

    less -S ${GC}/inputs/fastq/sequence.index

    In this file, we want the SAMPLE_NAME, FASTQ_FILE, RUN_ID, LIBRARY_NAME, CENTER_NAME, INSTRUMENT_PLATFORM (columns 10, 1, 15, 6, 13).

    • You can use perl/awk/linux to extract these fields & format as necessary.
    • I prepared a perl script that you can use:
    perl ${GC}/scripts/genIndex.pl > ${SETUP}/align.index

    Let's look at the index file: 

    less -S ${SETUP}/align.index

    Align index.png

    The command-line --fastq option or the configuration file FASTQ_PREFIX setting can be used to specify a prefix to the FASTQ1/FASTQ2 file paths.

    This file is specified either via the command-line --index_file parameter or via the configuration file INDEX_FILE 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}/gotcloud/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
    • Keys to override:
    Key Name Description
    Index File Settings - pointing GotCloud to your data
    INDEX_FILE Path to the FASTQ index file that you created
    • Alternatively, this can be specified on the command-line as --index_file
    FASTQ_PREFIX Prefix to be added to the FASTQ files in INDEX_FILE
    • Not required
    BAM_INDEX Path to the BAM index file
    • to be created by alignment
    • to be used for snp calling
    Reference File Settings - telling GotCloud where to find your reference files
    REF_DIR Path to your reference files
    • You don't have to use this, you can specify the full path for each file
    REF Path/filename of the FASTA reference file
    • If different than default: $(REF_DIR)/human.g1k.v37.fa
    DBSNP_VCF Path/filename of the DBSNP file
    • If different than default: $(REF_DIR)/dbsnp_135.b37.vcf.gz
    HM3_VCF Path/filename of the HapMap3 file
    • If different than default: $(REF_DIR)/hapmap_3.3.b37.sites.vcf.gz
    OMNI_VCF Path/filename of the OMNI file
    • If different than default: $(REF_DIR)/1000G_omni2.5.b37.sites.PASS.vcf.gz
    INDEL_PREFIX Path/filename base of the indels file
    • If different than default: $(REF_DIR)/1kg.pilot_release.merged.indels.sites.hg19

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

    more ${GC}/inputs/gotcloud.conf

    It already points to your align file.

    Run GotCloud Align

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

    ${GC}/gotcloud/gotcloud align --conf ${GC}/inputs/gotcloud.conf --numcs 2
    • --numcs means to run 2 samples at a time.
      • Depends on your system

    GcalignStart.png

    This should take < 4 minutes to run.

    It should end with a line like: Processing finished in 133 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.

    Examining GotCloud Align Output

    Let's look at the output directory: 

    ls ${OUTPUT}

    GcalignOutM.png

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

    ls ${OUTPUT}/bams

    GcalignOutBAMm.png

    BAM Files:

    • Binary Sequence Alignment/Map (SAM) Format
    • Maps reads to Chromosome/Position
    • For a detailed explanation of the SAM/BAM format, see:
    • Consists of:
      • Header
        • Starts with '@'
        • Records - one for each sequence read

    Let's examine a BAM file: 

    samtools view -h ${OUTPUT}/bams/

    BAM.png

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

    ls ${OUTPUT}/QCFiles

    GcalignOutQCm.png

    Check for sample contamination:

    • *.selfSM : Main output file containing the contamination estimate.
      • If you are only interested in checking sample contamination:
        • Check the 'FREEMIX' column for genotype-free estimate of contamination
        • Check the 'CHIPMIX' column for contamination estimates with external genotypes (if provided)
    • *.selfRG : Same output to .*selfSM, but separated by readGroup (which might be helpful for library-level examination)
    • *.depthSM : depth distribution of reads covering the marker position of the input VCF, across all readGroups.
    • *.depthRG : depth distribution of reads covering the marker position of the input VCF, per readGroups.
    less -S ${OUTPUT}/QCFiles/HG00551.genoCheck.selfSM

    Contam1.png

    Next, let's look at some quality control metrics: 

    cat ${OUTPUT}/QCFiles/HG00551.qplot.stats
    • 99.16% mapping rate
    • 94.01% high quality bases
    • 7x coverage
    • 31.3% A, 31.3% T
    • 18.7% C, 18.7% G


    Generate the pdf's of our quality metrics: 

    Rscript ${OUTPUT}/QCFiles/HG00551.qplot.R
Rscript ${OUTPUT}/QCFiles/HG00553.qplot.R
Rscript ${OUTPUT}/QCFiles/HG00640.qplot.R
Rscript ${OUTPUT}/QCFiles/HG00641.qplot.R

    Examine the PDF: 

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

    The first plot: Empirical vs reported Phred score does not look as good as we would like.

    • This is due to the small region used for recalibration

    Look at the PDF I produced when I ran the whole genome: 

    evince ${GC}/example/HG00551.wg.qplot.pdf&

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

    Retrieved from "http://genome.sph.umich.edu/w/index.php?title=SeqShop:_Sequence_Mapping_and_Assembly_Practical,_June_2014&oldid=10607"