Difference between revisions of "SeqShop: Sequence Mapping and Assembly Practical, June 2014"
Line 27: | Line 27: | ||
</div> | </div> | ||
− | == Examining GotCloud Align Input Files == | + | == Examining [[GotCloud]] Align Input Files == |
=== Examining Raw Sequence Reads : FASTQs === | === Examining Raw Sequence Reads : FASTQs === | ||
FASTQ : standard file format provided to you by those who did the sequencing. | FASTQ : standard file format provided to you by those who did the sequencing. | ||
Line 95: | Line 95: | ||
=== Reference Files === | === Reference Files === | ||
− | Reference files can be downloaded with GotCloud or from other sources | + | Reference files can be downloaded with [[GotCloud]] or from other sources |
* See [[GotCloud: Genetic Reference and Resource Files]] for more information on downloading/generating reference files | * See [[GotCloud: Genetic Reference and Resource Files]] for more information on downloading/generating reference files | ||
Revision as of 13:39, 15 June 2014
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
- 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.
- Open putty
- Start->Enter "putty" in the search and select "PuTTY" from the program list
- Configure PuTTY in the PuTTY Configuration window
- Host Name:
seqshop-server.sph.umich.edu
- 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
- 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 /home/mktrost/seqshop/setup.txt
- You won't see any output after running
source
- It silently sets up your environment
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 ${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?
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 the
space bar
to jump down a page
Use 'q'
to exit out of less
q
- Do you remember the parts of a FASTQ?
Look at the paired read:
less -S ${IN}/fastq/HG00551.SRR190851_2.fastq
Remember, use 'q'
to exit out of less
q
- Do you notice something in common?
Reference Files
Reference files can be downloaded with GotCloud or from other sources
- See GotCloud: Genetic Reference and Resource Files for more information on downloading/generating reference files
For alignment, you need:
- Reference genome FASTA file
- Contains the reference base for each position of each chromosome
- Additional information on the FASTA format: http://en.wikipedia.org/wiki/FASTA_format
- 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 ${REF}
Let's read the reference FASTA file (all reference bases for the chromosome):
less ${REF}/human.g1k.v37.chr22.fa
Remember, use 'q'
to exit out of less
q
- Where is the reference sequence?
GotCloud FASTQ Index File
The FASTQ index 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
- For recalibration we need to know which reads were in the same run.
FASTQ Index 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 index file I prepared for this tutorial:
less -S ${IN}/align.index
Remember, use 'q'
to exit out of less
q
- Which samples had multiple runs?
- Hard to read the index? Need a hint?
- Use cut to extract just the MERGE_NAME & RGID fields cut -f 1,4 ${IN}/align.index
How do you point GotCloud to your index file?
- Command-line
--index_file
option
- or
- 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}/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 ${IN}/gotcloud.conf
Use the space bar
to advance if the whole file isn't displayed.
- If your input and references are at different paths than specified, what would you change?
Run GotCloud Align
Now that we have all of our input files, we need just a simple command to run them
${GC}/gotcloud align --conf ${IN}/gotcloud.conf --numcs 2
--numcs
means to run 2 samples at a time.- How many you can run concurrently depends on your system
This should take 1-3 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 ${OUT}
Quality Control Files
Let's take a look at our quality control output directory:
ls ${OUT}/QCFiles
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
- Check the 'FREEMIX' column for genotype-free estimate of contamination
less -S ${OUT}/QCFiles/HG00551.genoCheck.selfSM
Remember, use 'q'
to exit out of less
q
- Is there evidence of sample contamination?
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?
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?
BAM Files
Binary Sequence Alignment/Map (SAM) Format
- Maps reads to Chromosome/Position
- For a detailed explanation of the SAM/BAM format, see:
- SAM/BAM Spec: http://samtools.github.io/hts-specs/SAMv1.pdf
- Additional information I put together as I started working with SAM/BAM: SAM
Let's look at the BAMs (aligned reads that are ready for variant calling):
ls ${OUT}/bams
Let's examine at the first 5 lines of the BAM file:
${GC}/bin/samtools view -h ${OUT}/bams/HG00551.recal.bam|head -n 5
- What are the chromosome and position of the first record in the BAM file?
Logging Off
To logout of seqshop-server, type:
exit
And close the windows.
When done, log out of the Windows machine.