| Table of Contents |
|---|
...
If you have completed 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.
Tip I strongly encourage everyone to complete these tutorials rather than relying on the canned data. I routinely work with raw fastq files, map them, and call variants on them. Outside of this course I use IGV maybe 1 time a year, and I don't think I'm an uncommon case. 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 language bash title This block is collapsed because I want to strongly strongly discourage anyone from jumping straight to this option. The mapping and variant calling tutorials are way way more important than this visualization tutorial. collapse true $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>@stampede2.tacc.utexas.edu:/corral-repl/utexas/BioITeam/gva_course/mapping/IGV . # to copy to a local computer skipping the step of copying to a stampede2 directory and secure copying from there.
If you are sure you this is how you want to proceed and another offer of assistance with any problems you may be having with the mapping or SNV tutorial, won't tempt you, then skip down to #Launching IGV.
...
IGV likes its reference genome files in GFF (Gene Feature Format) rather than the fasta or gbk formats we've been working with. While you may assume this is a job for our old friend bp_seqconvert.pl, that script actually doesn't deal with GFF files. So, we're going to use another tool for sequence format conversion called Readseq. Install the readseq package from bioconda to a new your GVA2021 conda environment to get started.
| Code Block | ||||||
|---|---|---|---|---|---|---|
| ||||||
conda create -n GVA-readseqactivate GVA2021 conda install -c bioconda readseq readseq --version |
...
readseq is a java based program which means it is envoked in a very different manner than anything we have worked with thus far. Luckily for us the conda package actually includes a wrapper allowing us to envoke the command simply by typing the readseq name like all the other programs we have worked with.
| Tip | ||
|---|---|---|
| ||
This is one of only 2 java based programs that this course covers. As the readseq wrapper conda provided makes this so much easier to envoke, we will use it. It is recommended to look back at a previous last years tutorial to see how this was handled without the wrapper incase you encounter a java based program in your own work that doesn't have such a helpful wrapper, and need to know where to start. |
...
| Code Block | ||
|---|---|---|
| ||
cds mkdir GVA_IGV_Tutorial cd GVA_IGV_Tutorial readseq $SCRATCH/GVA_bowtie2_mapping/NC_012967.1.gbk -f GFF -o NC_012967.1.gbk.gff |
A final oddity of the readseq program is that rather than displaying any kind of status message, or being silent when executed, the program actually displays the version of the readseq program itself. This is something that initially made me assume the conversion had failed. 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 using commands like head and tail.
| Tip | ||
|---|---|---|
| ||
You may notice that the output file appended a ".gff" ending to the ".gbk" ending rather than replacing it. This is can be done to demonstrate the order of operations performed on the file (in this case taking a gbk file and converting it to a gff file) a longer list of operations, such as sequential filtering a vcf file for frequency above 90%, with mapq scores above 20, on chromosome 7, between 10,000,000 and 190,000,000 bp might result in file contents looking like the following: sample.vcf some programs do not like unknown or effectively nonsense file endings, in which case you may need to append ".vcf" to the names above . |
Copy files to your local computer
...
Copy files to your local computer
Again, since IGV is an interactive graphical viewer program that we'll be running on our local computer, we need to get the files we want to visualize onto your desktop machine.
...
| Anchor | ||||
|---|---|---|---|---|
|
| Warning |
|---|
For the remainder of the tutorial, work on your local machine. NOT TACC! |
Load genome into IGV
| Note | ||
|---|---|---|
| ||
Most of the differences will be trivial differences between windows and mac versions, or difficulties in getting the identical level of zoom in/out on the genome, but there may be some that are important and impact learning. If you are unclear if the difference is important get my attention on zoom and we'll work through it. |
| Warning |
|---|
For the remainder of the tutorial, work on your local machine. NOT TACC! |
...
. |
From the main window of IGV:, click
...
on Genomes > Create .genome File...
- select the fasta version of your genome (NC_012967.1.fasta)
...
and you should be presented with the following window that you can fill out with information specific to your system.
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 must 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
...
- Choose Tools > Run igvtools....
- Choose "index" from the commands drop-down menu.
- Select the
SRR030257.vcffile for "Input File" - Click the "run" button in the middle of the window.
It will look like nothing has happened aside from the appearance of "Done" in the messages box, but you can now close the igvvtools "Run" window and choose choose File > Load from File. If you navigate to your IGV directory, you will now see a brand new SRR030257.vcf.idx file. You can now load the SRR030257.vcf file file, and it will show up as a new track near the top of your window. Note that it is the VCF file, not the newly made index file that you are selecting. like with the other files you loaded above, IGV requires you to have an index in the same directory as the file you are loading, but it is the file, not its index that you load.
...
.
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 as indexing BAM files in particular can be a computationally hefty task.
You are now free to investigate different areas and their alignments in the genome.
...
- Close IGV (if you have it open from the first tutorial with your mapping, SNV, and SV data) and reopen it.
- 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)
- Load the bam files you downloaded: File > Load from File… and select HCC1143.normal.21.19M-20M.bam
- Turn on dbSNP annotations File > Load from Server… >Annotation > Variation and Repeats > dbSNP 1.4.7
- Right click on the track name on the left and select sort alignments by start location
There are 2 mutations visible in the chr21:19,479,237-19,479,814 region answer the following questions:
Expand title Are both SNPS supported by reads mapping to both the forward and reverse DNA strand (hint: make sure reads are colored by strand)? Yes, both forward and reverse reads (red and blue if colored by strand) contain the SNPs compared to the reference
Expand title Which is more likely to be related to disease? why? The one on the left does not correspond to a dbSNP entry and is therefore more likely to be related to disease state.
There are 2 SNPs visible in the chr21:19,666,833-19,667,007 region. Answer the following questions:
Expand title Two mutations very close together is often a case of poor alignment scores. Is that the case here (remember this is human data)? 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
Expand title Is either likely to be related to disease? Neither is likely to be related to disease (or at least not to rare disease) as both mutations have previously been identified as naturally occurring by dbSNP
Expand title What is going on in the chr21:19,324,469-19,331,468 region? 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
Expand title What is going on in the chr21:19,102,154-19,103,108 region? 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.
Info title How can we identify this region as an AluY element There are several methods that could accomplish this. 1 as pointed out in class would be to pull the sequence from this region and blast that region, but that sounds like a lot of work. A second and much better way would be to turn on some additional tracks in IGV. In this case Alu elements are identified as 'SINE' elements in the repeat masker data base.
To turn this on, File > Load from Server… >Annotation > Variation and Repeats > Repeat Masker > SINEInformation about what Repeat masker is doing and where it comes from can be found here: http://genome.ucsc.edu/cgi-bin/hgTrackUi?db=hg19&g=rmsk
What other interesting things can you find?
An additional tutorial from another group working with the same human data can be found here if interested.
...
- Tablet - a lightweight NGS data browser
- Visualize mapped data at UCSC genome browser
In your work if you come across another (particularly if you find it very useful, or easier to use than IGV) I'd love to hear from you.