In this tutorial, we're going to view the aligned reads and variants that we called in the past two lessons in the Integrative Genomics Viewer from the Broad Center. You'll need the output from Introduction to mapping (bowtie, BWA) and Introduction to variant calling (SAMtools).the Mapping tutorial and/or the Variant calling tutorial.
| Table of Contents |
|---|
Getting situated
You can start this tutorial two ways
...
If you do not have the results, you can copy some "canned" ones that we provide using these commands. Remove any existing intro_to_mapping directory in your $SCRATCH space before running these commands, if you have one.
| Code Block |
|---|
cds cp -r /corral-repl/utexas/BioITeam/$BI/ngs_course/intro_to_mapping/data intro_to_mapping cd intro_to_mapping cp -r /corral-repl/utexas/BioITeam$BI/ngs_course/intro_to_mapping/samtools_* . cp -r /corral-repl/utexas/BioITeam/$BI/ngs_course/intro_to_mapping/comparison . |
...
IGV likes its reference genome files in GFF (Gene Feature Format). Unfortunately, our old friend bp_seqconvert.pl doesn't do GFF. Fortunately, it's cousin bp_genbank2gff3.pl does.
| Code Block |
|---|
module load bioperl cp /corral-repl/utexas/BioITeam/ngs_course/scripts/bp_genbank2gff3.pl . ./bp_genbank2gff3.pl NC_012967.1.gbk |
| Expand | ||||
|---|---|---|---|---|
| ||||
NC_012967.1.gbk.gff |
...
| Expand | |||||
|---|---|---|---|---|---|
| |||||
Another useful trick with either IGV or UCSC: displaying your own BLAST results: BioPerl allows for super-easy conversion from blast output to a gff file; IGV and the UCSC browser both understand GFF files. The short script Let's use the blast result we had from the earlier test for the JAG1 gene to show you how. You'll need to provide the input file - it's the ".oNNNNNN" output file from your blast job.
The resulting jag1_blast.out.gff can be moved to your local machine and opened in IGV. Load the human reference first though! |
...
| Code Block |
|---|
wget http://www.broadinstitute.org/igv/projects/downloads/IGV_2.1.1422.zip unzip IGV_2.1.1422.zip cd IGV_2.1.1422 java -Xmx2g -jar igv.jar |
Load genome into IGV
From the main window of IGV, click on File -> Import Genome and you should be presented with the following window.
Enter the ID and Name of the Genome you are working with and select the path to your *.fasta file (the index, *.fai file needs to be in the same directory), then select the path to your *.gff file for the Gene File.
...
- Why are some reads different colors? Hint: Try changing the display options to show read pairs and editing some of the distance constraints.
- What is a typical mapping quality (MQ) for a read? Convert this to the probability that it is mismapped.
Expand title Remember the formula for a Phred quality score? The estimated probability that a read is mapped incorrectly is 10^(-MQ/10).
- Can you find a variant where the sequenced sample differs from the reference? This is going to be like looking for a needle in a haystack. Fortunately, we are going to learn how to use variant callers tomorrow and then we'll be able to zoom right to areas where there are discrepancies between reads and the reference genome that might indicate there were mutations in the sequenced E. coli.
Load variant calls into IGV
We're really interested in places in the genome where we think there are mutations. You can load the VCF file to check out those spots, but first you need to (guess what?) index it.
...
- Check out the rbsA gene region? What's going on here?
- What is going on in the pykF gene region? You might see red read pairs. What does that mean?
- Go to coordinate 2,698,092. Compare the bowtie and BWA alignments. Can you explain what's going on here?
- Go to coordinate 475,263. Compare the bowtie and BWA alignments. What's the story here?
Advanced exercise: human data scavenger hunt
Data from the CEU trio from the 1000 Genomes Project can be found directly from the Broad's server for IGV.
...