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 the Mapping tutorial and/or the Variant calling tutorial.

Getting situated

You can start this tutorial two ways:

1. If you have an intro_to_mapping directory with output from the previous tutorials with aligned reads or variant calls, then you should use those files for this tutorial.

2. If you do not have any results, you can copy some "canned" ones that we provide. Copy the entire contents of this directory back to your local machine:

$BI/intro_to_mapping/igv

Then skip down to #Launching IGV.

Prepare a GFF feature file for the reference sequence

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 with this command:

bp_genbank2gff3.pl NC_012967.1.gbk

Take a look at the original Genbank file and the new GFF3 file and try to get a handle on what is going on in this conversion.

(NB: I had to "fix" bp_genbank2gff3.pl to work on our GenBank file, by removing a line in the Perl script that caused an error. If someone knows of an easier way to get a GFF file from a file downloaded from Genbank, please share!)

grep '^gi' blast_jag1.o586038 > jag1_blast.out
module load perl
module load bioperl
bl2gff.pl jag1_blast.out > jag1_blast.out.gff

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.

We're going to copy all of these into a new directory called IGV to make it easier to just transfer the ones that we want. Change into your intro_to_mapping directory and run these commands if you just came from mapping:

mkdir IGV
cp NC_012967.1.gbk.gff IGV
cp bowtie/bowtie.sam IGV/bowtie.sam
cp bowtie2/bowtie2.sam IGV/bowtie2.sam
cp bwa/SRR030257.sam IGV/bwa.sam

mkdir IGV
cp NC_012967.1.gbk.gff IGV
cp samtools_bowtie/NC_012967.1.fasta IGV
cp samtools_bowtie/NC_012967.1.fasta.fai IGV
cp samtools_bowtie2/SRR030257.sorted.bam IGV/bowtie2.sorted.bam
cp samtools_bowtie2/SRR030257.sorted.bam.bai IGV/bowtie2.sorted.bam.bai
cp samtools_bowtie/SRR030257.sorted.bam IGV/bowtie.sorted.bam
cp samtools_bowtie/SRR030257.sorted.bam.bai IGV/bowtie.sorted.bam.bai
cp samtools_bwa/SRR030257.sorted.bam IGV/bwa.sorted.bam
cp samtools_bwa/SRR030257.sorted.bam.bai IGV/bwa.sorted.bam.bai
cp comparison/* IGV

Now, copy this entire IGV directory back to your local Desktop machine.

pwd

scp -r <username>@lonestar.tacc.utexas.edu:</full/path/to/IGV/> .

Launching IGV

There are two ways; Launching IGV in your web browser or by downloading the binaries locally and running IGV from your machine.

In Browser

Navigate a web browser to this page:http://www.broadinstitute.org/software/igv/download

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

wget http://www.broadinstitute.org/igv/projects/downloads/IGV_2.1.22.zip
unzip IGV_2.1.22.zip
cd IGV_2.1.22
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.

Load mapped reads into IGV

From the main window of IGV, click on File -> Load 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:

• Navigate more quickly. The 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 and control-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.

  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.

You can do this from within IGV:

1. Choose File -> Run igvtools....
2. Choose "index" from the commands drop-down menu.
3. Select your *.vcf file (Ex: bowtie.vcf) for "Input File"
4. 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

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

Find the dbSNP accession number for the SNP apparent in the two 1000 genomes project trios in the intron between exons 8 and 9 of the GABBR1 gene.

Steps:

1. Download and install the Integrative Genome Viewer from the Broad Institute.
2. Select "Human hg18" as the reference genome
3. Get some data: File -> “Load from Server…” -> 1000 genomes -> CEU and YRI trios
4. Find the gene and the right exons
5. Zoom in until you find the SNP
6. Load and look at the SNP track: File -> Load from server -> Annotations -> Variants -> dbSNP

This is whole genome coverage data; later we'll look at exome data.