Overview

The Integrative Genomics Viewer (IGV) from the Broad Center allows you to view several types of data files involved in any NGS analysis that employs a reference genome, including how reads from a dataset are mapped, gene annotations, and predicted genetic variants.

Learning Objectives

In this tutorial, we're going to learn how to do the following in IGV:

• Create a custom genome database (usually used for microbial genomes) or load a pre-existing genome assembly (usually used for the genomes of model organisms and higher Eukaryotes).
• Load output from mapping reads to a reference genome.
• Load output from calling genetic variants.
• Navigate the view of the genome and interpret the display of this data.

Theory

Because NGS datasets are very large, it is often impossible or inefficient to read them entirely into a computer's memory when searching for a specific piece of data. In order to more quickly retrieve the data we are interested in analyzing or viewing, most programs have a way of treating these data files as databases. Database indexes enable one to rapidly pull specific subsets of the data from them.

The Integrative Genomics Viewer is a program for reading several types of indexed database information, including mapped reads and variant calls, and displaying them on a reference genome. It is invaluable as a tool for viewing and interpreting the "raw data" of many NGS data analysis pipelines.

Workflow 1: Viewing E. coli data in IGV

Data files

You can start this tutorial two ways:

1. If you have a mapping directory with output from the Mapping tutorial and the SNV calling tutorial, then you should use those files for part 1 of this tutorial. You can proceed with either one alone or with both.

2. Expand
   title If you have not done the other tutorials and want a "canned" data set provided for you, click here for example files.
   Code Block $BI/gva_course/mapping/IGV # location of example files cp -r /corral-repl/utexas/BioITeam/gva_course/mapping/IGV . # example command to copy to current directory scp -r username@ls5.tacc.utexas.edu:/corral-repl/utexas/BioITeam/gva_course/mapping/IGV . # to copy to a local computer skipping the step of copying to a lonestar directory and secure copying from there. 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 deal with GFF. So, we're going to show you another tool for sequence format conversion called Readseq. We've already installed it into the $BI/bin directory so you don't have to, but here we provide the steps that can be used to install it in a local directory.

...

To use it you need to first download the file readseq.jar linked from here. To get this onto TACC easily, use:

...

Overview

The Integrative Genomics Viewer (IGV) from the Broad Center allows you to view several types of data files involved in any NGS analysis that employs a reference genome, including how reads from a dataset are mapped, gene annotations, and predicted genetic variants.

Learning Objectives

In this tutorial, we're going to learn how to do the following in IGV:

• Create a custom genome database (usually used for microbial genomes) or load a pre-existing genome assembly (usually used for the genomes of model organisms and higher Eukaryotes).
• Load output from mapping reads to a reference genome.
• Load output from calling genetic variants.
• Navigate the view of the genome and interpret the display of this data.

Theory

Because NGS datasets are very large, it is often impossible or inefficient to read them entirely into a computer's memory when searching for a specific piece of data. In order to more quickly retrieve the data we are interested in analyzing or viewing, most programs have a way of treating these data files as databases. Database indexes enable one to rapidly pull specific subsets of the data from them.

The Integrative Genomics Viewer is a program for reading several types of indexed database information, including mapped reads and variant calls, and displaying them on a reference genome. It is invaluable as a tool for viewing and interpreting the "raw data" of many NGS data analysis pipelines.

Workflow 1: Viewing E. coli data in IGV

Data files

You can start this tutorial two ways:

1. If you have a mapping directory with output from the Mapping tutorial and the SNV calling tutorial, then you should use those files for part 1 of this tutorial. You can proceed with either one alone or with both.

2. Expand
   title If you have not done the other tutorials and want a "canned" data set provided for you, click here for example files.
   Code Block $BI/gva_course/mapping/IGV # location of example files cp -r /corral-repl/utexas/BioITeam/gva_course/mapping/IGV . # example command to copy to current directory scp -r username@ls5.tacc.utexas.edu:/corral-repl/utexas/BioITeam/gva_course/mapping/IGV . # to copy to a local computer skipping the step of copying to a lonestar directory and secure copying from there. 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 deal with GFF. So, we're going to show you another tool for sequence format conversion called Readseq. We've already installed it into the $BI/bin directory so you don't have to, but here we provide the steps that can be used to install it in a local directory.

wget https://sourceforge.net/projects/readseq/files/readseq/2.1.19/readseq.jar 

Readseq is written in java which makes it a little more complicated to use, but the general command to run the software is one of these (note that you do need to include the entire path, not just the "readseq.jar" name):

...

To do the conversion that we want, use this command:

...

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 at TACC. Indexing BAM files can be a computationally hefty task. 

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: SRR030257.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

...

Interested in determining the probability that a read is not where it should be? What is a typical mapping quality (MQ) for a read?

Can you find a variant where the sequenced sample differs from the reference? This would be like looking for a needle in a haystack if not for the use of variant callers and the control-f and control-b options 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
  title Coordinate 161,041. What gene is this in and what is the effect on the protein sequence?
  Gene is pcnB, mutation is a snp

• Expand
  title Coordinate 3,248,957. What gene is this in and what is the effect on the protein sequence?
  Gene is infB, mutation is a snp

• Expand
  title Coordinate 3,894,997. What type of mutation is this?
  Deletion of the rbsD gene

...

Navigate to coordinate 3,289,962. Compare the results for different alignment programs and settings. Can you explain what's going on here?

...

There is a 16 base deletion in the 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?

...

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.

...

Workflow 2: Viewing Human Genome Data in IGV

Now that you've familiarized yourself with IGV using a "simple" bacteria, let's look at something a "little" more complex: the human genome.

Advanced exercise: human data scavenger hunt 

Data from the 1000 Genomes Project can be found directly from the Broad's server for IGV. There are now MANY genomes available this way.

Find one or more dbSNP accession numbers for SNPs apparent in one of the two 1000 genomes project trios in the GABBR1 gene.

Steps:

1. Download and install the Integrative Genome Viewer from the Broad Institute. 
2. Select "Human hg19" as the reference genome from the top left drop down (you may need to select "more" to have hg19 as an option)
3. Get some data: File > Load from Server… > 1000 genomes > Alignments > ACB > exome > HG01880
4. Navigate to the rightmost exons of the GABBR1 gene.
5. Zoom in until you find some SNPs. (Hint look just to the left of the 2nd exon).
6. What type of library is this? (Hint: zoom out)
7. If you knew this was a cancer patient, consider how strongly you would think this may be a potentially causative mutation.
8. Imagine it was actually in the exon rather than just into the intron... would that make you consider it more?
9. Load and look at the SNP track: File > Load from server > Annotations > Variants and Repeats > dbSNP 1.4.7
10. The track may load with the Refseq genes, making it useful to resize that window to view both the gene and the dbSNP information simultaneously.
11. Consider if this makes you think it more likely or less likely that this is a causative mutation.

Optional Tutorial Exercises ...

samtools faidx NC_012967.1.fasta
samtools view -b -S -o bowtie/SRR030257.bam bowtie/SRR030257.sam
samtools sort bowtie/SRR030257.bam -o bowtie/SRR030257.sorted
samtools index bowtie/SRR030257.sorted.bam

...

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 a hypothetical test for the JAG1 gene to show you how. You'll need to provide the output file from your blast job as the input to the script (in this case "blast_jag1.058638".

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

...

. Indexing BAM files can be a computationally hefty task. 

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: SRR030257.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.

• Interested in determining the probability that a read is not where it should be? What is a typical mapping quality (MQ) for a read?

  Expand
  title Click here for the formula.
  The estimated probability that a read is mapped incorrectly is 10^(-MQ/10). Where MQ is the mapping quality.

• Can you find a variant where the sequenced sample differs from the reference? This would be like looking for a needle in a haystack if not for the use of variant callers and the control-f and control-b options 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
  title Some interesting example coordinates
  Expand title Coordinate 161,041. What gene is this in and what is the effect on the protein sequence? Gene is pcnB, mutation is a snp Expand title Coordinate 3,248,957. What gene is this in and what is the effect on the protein sequence? Gene is infB, mutation is a snp Expand title Coordinate 3,894,997. What type of mutation is this? Deletion of the rbsD gene Expand title Check out the rbsA gene region? What's going on here? There was a large deletion. Can you figure out the exact coordinates of the endpoints? Navigate to coordinate 3,289,962. 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 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. See if you can find more interesting locations. There are ~40 mutations total in this sample MOST of which are false positives.

Workflow 2: Viewing Human Genome Data in IGV

Now that you've familiarized yourself with IGV using a "simple" bacteria, let's look at something a "little" more complex: the human genome.

Advanced exercise: human data scavenger hunt 

Throughout this class we have used the wget command several times to download files to TACC which is running bash on a linux environment. As IGV only works on your local computer rather than TACC, downloading files directly to your laptop is more practical. Unfortunately wget is a linux only command and unix does not have it so instead we shall use the curl command to download some human data to IGV and look at it. The following commands should be run in a terminal window that is not logged into tacc.

cd ~/Desktop
mkdir human_IGV
cd human_IGV
curl -O https://xfer.genome.wustl.edu/gxfer1/project/gms/testdata/bams/hcc1143/HCC1143.normal.21.19M-20M.bam
curl -O https://xfer.genome.wustl.edu/gxfer1/project/gms/testdata/bams/hcc1143/HCC1143.normal.21.19M-20M.bam.bai

These if you look at the file names, you may notice that this bam file and its index correspond to human chromosome 21 from 19 million to 20 million bases. This limited data set is to hopefully avoid IGV crashing. 

Steps:

1. Close IGV (if you have it open from the first tutorial with your mapping, SNV, and SV data) and reopen it. 
2. Select "Human hg19" as the reference genome from the top left drop down (you may need to select "more" to have hg19 as an option)
3. Load the bam files you downloaded: File > Load from File…  and select HCC1143.normal.21.19M-20M.bam
4. Turn on dbSNP annotations File > Load from Server… > Tutorials > Variants > dbSNP 1.3.1
5. Right click on the track name on the left and select sort alignments by start location
6. There are 2 mutations visible in the chr21:19,479,237-19,479,814 region answer the following questions:
   1. Are both SNPS supported by reads mapping to both the forward and reverse DNA strand (hint: make sure reads are colored by strand)?
      

   2. Which is more likely to be related to disease? why?

      Expand
      title Answers
      a. Yes, both forward and reverse reads (red and blue if colored by strand) contain the SNPs compared to the reference b. The one on the left does not correspond to a dbSNP entry and is therefore more likely to be related to disease state

      
      

7. There are 2 SNPs visible in the chr21:19,666,833-19,667,007 region. Answer the following questions:

   1. Two mutations very close together is often a case of poor alignment scores. Is that the case here (remember this is human data)?

   2. Is either likely to be related to disease? 

      Expand
      title Answers
      a. No, each read only has 1 mutation on it, these are 2 different alleles each with its own SNP relative to 'wt'. Both are reported in dbSNP b. Neither is likely to be related to disease or at least to rare disease as both mutations have previously been identified as naturally occurring by dbSNP

      
      

8. What is going on in the chr21:19,324,469-19,331,468 region?

   Expand
   title Answers
   Homozygous deletion. In the track on the left, right click and select 'view as pairs' to see linkage between R1 and R2 to see individual reads mapping to both sides of the deletion

9. What is going on in the chr21:19,102,154-19,103,108 region?

   Expand
   title Answers
   This is an example of poor alignment to a repetitive AluY element. Notice how of the read pairs that map with numerous SNPs have 1 read that maps with lots of SNPs and the other read maps with none? This is caused by mapping reads to a limited area of the whole genome, if these reads had been allowed to map to the entire genome it is very likely that both read pairs would map without SNPs somewhere else in the genome.

10. What other interesting things can you find?

Optional Tutorial Exercises ...

samtools faidx NC_012967.1.fasta
samtools view -b -S -o bowtie/SRR030257.bam bowtie/SRR030257.sam
samtools sort bowtie/SRR030257.bam -o bowtie/SRR030257.sorted
samtools index bowtie/SRR030257.sorted.bam

You can use IGV to visualize mapped reads and predicted variants from any later tutorial!

...