In this tutorial, we're going to view the aligned reads and variants that we called in the past two lessons in the Integrated Genomics Viewer from the Broad Center. You'll need the output from Introduction to mapping (bowtie, BWA) and Introduction to variant calling (SAMtools).
You can start this tutorial two ways
1. If you have an introduction_to_mapping directory from doing the previous tutorials, then you can use it for this tutorial.
cd $SCRATCH/intro_to_mapping |
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.
cds cp -r /corral-repl/utexas/BioITeam/ngs_course/intro_to_mapping/data intro_to_mapping cd intro_to_mapping cp -r /corral-repl/utexas/BioITeam/ngs_course/intro_to_mapping/samtools_* . cp -r /corral-repl/utexas/BioITeam/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.
module load bioperl cp /corral-repl/utexas/BioITeam/ngs_course/scripts/bp_genbank2gff3.pl . ./bp_genbank2gff3.pl NC_012967.1.gbk |
|
NC_012967.1.gbk.gff |
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!)
|
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! |
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:
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.
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_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.
|
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:
|
There are two ways; Launching IGV in your web browser or by downloading the binaries locally and running IGV from your machine.
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.
wget http://www.broadinstitute.org/igv/projects/downloads/IGV_2.1.14.zip unzip IGV_2.1.14.zip cd IGV_2.1.14 java -Xmx2g -jar igv.jar |
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.
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.
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:
bowtie.vcf) for "Input File"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 index BAM and FASTA files the same way inside of IGV if you haven't already created indexes for them.
There are a lot of things you can do in IGV. Here are a few:
control-f and control-b to jump forward and backward within that list of features. Try this on the variant calls track.page-up page-down, home, end.See the IGV Manual for more tips and how to load other kinds of data.
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:
This is whole genome coverage data; later we'll look at exome data.
|
rs29220 |