...
Take a look at the contents of the original Genbank file and the new GFF file and try to get a handle on what is going on in this conversion.
If you have only done the mapping tutorial and NOT the variant calling tutorial
You will need to index your reference FASTA and convert your SAM output files into sorted and indexed BAM files. The "why?" behind these steps is described more fully in the Variant calling tutorial. If you are in your intro_to_mapping
directory, these commands will perform the necessary steps.
Warning | ||
---|---|---|
| ||
Repeat the last three commands for each SAM output file that you want to visualize in IGV. |
Copy files to your desktop
IGV is an interactive graphical viewer program. You can't run it on TACC, so we need to get the relevant files back to your desktop machine.
They include:
- Indexed reference FASTA files
- GFF reference sequence feature files
- Sorted and indexed mapped read BAM files
- VCF result files
- ... and possibly many other types of files.
The easiest way to to this is probably to copy everything you want to transfer into a new directory called IGV
. Since many of the tutorial output files had the same names (but resided in different directories) be careful to give them unique destination names when you copy them into the new directory together.
For starters, you could change into your intro_to_mapping
directory and run commands like these if you just came from the Mapping tutorial:
Code Block | |||||
---|---|---|---|---|---|
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! |
If you have only done the mapping tutorial and NOT the variant calling tutorial
You will need to index your reference FASTA and convert your SAM output files into sorted and indexed BAM files. The "why?" behind these steps is described more fully in the Variant calling tutorial. If you are in your mapping
directory, these commands will perform the necessary steps.
Warning | |||||||||
---|---|---|---|---|---|---|---|---|---|
| |||||||||
|
Now, copy this entire IGV directory back to your local Desktop machine.
...
In the terminal connected to Lonestar, figure out the complete path to the IGV directory.
Code Block |
---|
pwd
|
Open a new terminal window on your Desktop. Fill in the parts in brackets <> in this command:
Code Block |
---|
scp -r <username>@lonestar.tacc.utexas.edu:</full/path/to/IGV/> .
|
...
Launching IGV
Warning |
---|
For the remainder of the tutorial, work on your local machine. NOT TACC! |
There are two ways; Launching IGV in your web browser or by downloading the binaries locally and running IGV from your machine.
Locally on the classroom machines booted in Linux
This downloads the IGV executable and tells the command line to launch it (via the java command).
Code Block |
---|
wget http://www.broadinstitute.org/igv/projects/downloads/IGV_2.3.5.zip
unzip IGV_2.3.5.zip
cd IGV_2.3.5
java -Xmx2g -jar igv.jar
|
In a Web browser
Navigate a web browser to this page:http://www.broadinstitute.org/software/igv/download. You will need to register your email address to use this option!
Go ahead and click on the "Launch with 2 GB" option. This will download a "Java Web Start" file that you can launch by locating it on your Desktop and double-clicking.
Locally on your own Mac or Windows computer
Use this link to download IGV:
Repeat the last three commands for each SAM output file that you want to visualize in IGV. |
Copy files to your desktop
IGV is an interactive graphical viewer program. You can't run it on TACC, so we need to get the relevant files back to your desktop machine.
They include:
- Indexed reference FASTA files
- GFF reference sequence feature files
- Sorted and indexed mapped read BAM files
- VCF result files
- ... and possibly many other types of files.
The easiest way to to this is probably to copy everything you want to transfer into a new directory called IGV
. Since many of the tutorial output files had the same names (but resided in different directories) be careful to give them unique destination names when you copy them into the new directory together.
For starters, you could change into your mapping
directory and run commands like these if you just came from the mapping tutorial:
Code Block |
---|
mkdir IGV
cp NC_012967.1.fasta IGV
cp NC_012967.1.fasta.fai IGV
cp NC_012967.1.gbk.gff IGV
cp bowtie/SRR030257.sorted.bam IGV/bowtie.sorted.bam
cp bowtie/SRR030257.sorted.bam.bai IGV/bowtie.sorted.bam.bai
...
|
Now, copy this entire IGV directory back to your local Desktop machine.
Expand | ||||
---|---|---|---|---|
| ||||
In the terminal connected to Lonestar, figure out the complete path to the IGV directory.
Open a new terminal window on your Desktop. Fill in the parts in brackets <> in this command:
|
Anchor | ||||
---|---|---|---|---|
|
Launching IGV
Warning |
---|
For the remainder of the tutorial, work on your local machine. NOT TACC! |
There are two ways; Launching IGV in your web browser or by downloading the binaries locally and running IGV from your machine.
Locally on the classroom machines booted in Linux
This downloads the IGV executable and tells the command line to launch it (via the java command).
Code Block |
---|
wget http://www.broadinstitute.org/igv/projects/downloads/IGV_2.3. |
...
After unzipping, you should be able to click on igv.bat
for Windows or igv.command
on MacOSX to lauch IGV. If this is not working, you might need to try the web start.
Load genome into IGV
From the main window of IGV, click on Genomes ? Create .genome File... and you should be presented with the following window.
Enter the ID and Name of the Genome you are working with (these can be anything that makes sense to you) 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.
Load mapped reads into IGV
32.zip
unzip IGV_2.3.32.zip
cd IGV_2.3.32
java -Xmx2g -jar igv.jar
|
In a Web browser
Navigate a web browser to this page:http://www.broadinstitute.org/software/igv/download. You will need to register your email address to use this option!
Go ahead and click on the "Launch with 2 GB" option. This will download a "Java Web Start" file that you can launch by locating it on your Desktop and double-clicking.
Locally on a Mac or Windows computer
Use this link to download IGV:
http://www.broadinstitute.org/igv/projects/downloads/IGV_2.3.32.zip
After unzipping, you should be able to click on igv.bat
for Windows or igv.command
on MacOSX to lauch IGV. If this is not working, you might need to try the web start.
Load genome into IGV
From the main window of IGV, click on File ? Load from File.... Choose bowtie.sorted.bam
After importing your reference genome and loading an alignment file, your screen should look similar to the following:
And you are now free to investigate different areas and their alignments in the genome.
Navigating in IGV
There are a lot of things you can do in IGV. Here are a few:
- Zoom in using the slider in the upper right. Do this until you see mapped reads and finally individual bases appear.
- Navigate by clicking and dragging in the window. This is how you move left and right along the genome.
- Navigate more quickly. Use
page-up
page-down
,home
,end
. - Jump to the next point of interest. Click on a track name on the left side of the window (Ex: bowtie.vcf), to select it. You can then use
control-f
andcontrol-b
to jump forward and backward within that list of features. Try this on the variant calls track. - Jump right to a gene. (If you have gene features loaded.) Type its name into the search box. Try "topA".
- Load multiple BAM alignments or VCF files at once. Try this to compare a few different regions between the bowtie and BWA results.
- Change the appearance of genes. Right click on the gene track and try "expanded". Experiment with the other options.
- Change the appearance of reads. Right click on a BAM track and choose "show all bases" and "expanded". Experiment with the other options.
See the IGV Manual for more tips and how to load other kinds of data.
Exercises
...
What is a typical mapping quality (MQ) for a read? Convert this to the probability that it is mismapped.
Expand | ||
---|---|---|
| ||
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.
...
Genomes > Create .genome File... and you should be presented with the following window.
Enter the ID and Name of the Genome you are working with (these can be anything that makes sense to you) 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. Click OK and then save this *.genome file inside the same folder as your data.
Load mapped reads into IGV
From the main window of IGV, click on File > Load from File.... Choose bowtie.sorted.bam
After importing your reference genome and loading an alignment file, click on the + button in the upper right until reads appear! Then, your screen should look similar to the following:
And you are now free to investigate different areas and their alignments in the genome.
Navigating in IGV
There are a lot of things you can do in IGV. Here are a few:
- Zoom in using the slider in the upper right. Do this until you see mapped reads and finally individual bases appear.
- Navigate by clicking and dragging in the window. This is how you move left and right along the genome.
- Navigate more quickly. Use
page-up
page-down
,home
,end
. - Jump to the next point of interest. Click on a track name on the left side of the window (Ex: bowtie.vcf), to select it. You can then use
control-f
andcontrol-b
to jump forward and backward within that list of features. Try this on the variant calls track. - Jump right to a gene. (If you have gene features loaded.) Type its name into the search box. Try "topA".
- Load multiple BAM alignments or VCF files at once. Try this to compare a few different regions between the bowtie and BWA results.
- Change the appearance of genes. Right click on the gene track and try "expanded". Experiment with the other options.
- Change the appearance of reads. Right click on a BAM track and choose "show all bases" and "expanded". Experiment with the other options.
See the IGV Manual for more tips and how to load other kinds of data.
Exercises
- 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.
Expand Some interesting locations to look at for the time being... Some interesting locations to look at for the time being... - Coordinate 161,041. What gene is this in and what is the effect on the protein sequence?
- Coordinate 3,248,957. What gene is this in and what is the effect on the protein sequence?
- Coordinate 4,015,892. What is different about the reads mapped to this location?
- Coordinate 3,894,997. What type of mutation is this?
- Coordinate 1,733,647. What's going on here?
- See if you can find more interesting locations. There are ~40 mutations total in this sample.
...
You can do this from within IGV:
- Choose Tools ? > Run igvtools....
- Choose "index" from the commands drop-down menu.
- Select your *.vcf file (Ex:
bowtie.vcf
) for "Input File" - Click the "run" button.
It will look like nothing has happened, but you can now close the "Run" window and choose File ? > Load File. If you navigate to your IGV directory, you will now see a brand new bowtie.vcf.idx
file. You can now load the file bowtie.vcf
, and it will show up as a new track near the top of your window.
Tip: You can also index BAM and FASTA files the same way inside of IGV if you haven't already created indexes for them. But, it's usually easier and quicker to do this on the command line .
Exercises
...
at TACC. Indexing BAM files can be a computationally hefty task.
Exercises
Check out the rbsA gene region? What's going on here?
Expand Answer Answer There was a large deletion. Can you figure out the exact coordinates of the endpoints?
Navigate to coordinate 4753,289,288962. Compare the results for different alignment programs and settings. Can you explain what's going on here?
Expand Answer Answer There is a 16 base deletion in the gltA gltB gene reading frame.
What is going on in the pykF gene region? You might see red read pairs. What does that mean? Can you guess what type of mutation occurred here?
Expand Answer Answer The read pairs are discordantly mapped. There was an insertion of a new copy of a mobile genetic element (an IS150 element) that exists at other locations in the reference sequence.
...
If you've made it through the other exercises on your own data, take a look at some human genome re-sequencing data where the files can be loaded directly from public databases.
Advanced exercise: human data scavenger hunt
Anchor | ||||
---|---|---|---|---|
|
See this page for the human data scavenger hunt
...
- Download and install the Integrative Genome Viewer from the Broad Institute.
- Select "Human hg18" or "Human hg19" as the reference genome
- Get some data: File - > Load from Server… - > 1000 genomes - > Alignments - > CEU Trio WGS - > select those 3 samples
- Navigate to the rightmost exons of the GABBR1 gene
- Zoom in until you find some SNPs - they might be in exons or introns; there is also at least one example of a short insertion variant between exons 2 and 3
- Load and look at the SNP track: File -> Load from server -> Annotations -> Variants and Repeats -> dbSNP
This is whole genome coverage data; later we'll look at exome data.
...
rs29220, rs29222, rs28359988, rs76688565, there might be more in the locus; I got tired of looking.
Is there an alternate allele in the child which correlates with one or both of the parents? (i.e. - do genetics work?)
...
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 bl2gff.pl
does the conversion.
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.
Code Block | ||
---|---|---|
| ||
grep '^gi' blast_jag1.o586038 > jag1_blast.out
module load perl
module load bioperl
bl2gff.pl jag1_blast.out > jag1_blast.out.gff
|
...
- > Load from server > Annotations > Variants and Repeats > dbSNP
This is whole genome coverage data; later we'll look at exome data.
Expand | ||||
---|---|---|---|---|
| ||||
rs29220, rs29222, rs28359988, rs76688565, there might be more in the locus; I got tired of looking. |
Is there an alternate allele in the child which correlates with one or both of the parents? (i.e. - do genetics work?)
From here...
You can use IGV to visualize mapped reads and predicted variants from any later tutorial!
...