ABRF2010 Tutorial

From genomewiki
Jump to navigationJump to search

Tutorial: Viewing sequencing data in the UCSC Genome Browser

Presented March 20, 2010 as part of the workshop Tools to Facilitate Management, Analysis and Visualization of 2nd Generation Sequencing Data at ABRF 2010.

This page serves as a step-by-step guide to the demo of UCSC Genome Browser tools presented at the workshop. Earlier in the workshop, Charlie Nicolet shared some sequencing and analysis results on chromosome 21 for ChIP-Seq of RNA polymerase II to HeLa cells. Thanks to Charlie for allowing us to post copies of his data and use them as a basis for this tutorial.

Genome Graphs with .sgr data

Some of the result files from the ChIPSeq toolset that Charlie Nicolet showed us earlier are in the .sgr format, i.e. tab-separated text that specifies chromosome, position and a numeric signal value. This format can be uploaded to UCSC's Genome Graphs tool for genome-at-a-glance view of the data.

Steps

  1. In a new web browser window or tab, go to http://genome.ucsc.edu/ .
  2. On the left side of the page, there is a blue sidebar. Click on "Genome Graphs", the 7th link from the top of the sidebar.
  3. In the assembly menu, change the selection to "Mar. 2006 (NCBI36/hg18)".
  4. Click the upload button.
  5. In the input labeled "name of data set", type or paste "redbin".
  6. Paste this URL: http://genome-test.cse.ucsc.edu/ABRF2010/chr21_extended.txt_redbin.sgr in the box labeled "Paste URLs or data", and click submit.
  7. You should see confirmation screen; click OK to return to the Genome Graphs display.
  8. Now you should see a layout of human chromosomes. Genome Graphs doesn't show the data until you select it. In the graph menu on the left, change "--nothing--" to "redbin 1". (If you don't see blue data appearing over chromosome 21, go to http://genome.ucsc.edu/cgi-bin/hgGenome?hgS_doOtherUser=submit&hgS_otherUserName=AngieHinrichs&hgS_otherUserSessionName=ABRF2010_GG to catch up.)
  9. Now the data appear -- in this case, for chr21 only. Click on a peak.
  10. Click on a peak to jump to the Genome Browser, displaying a 1,000,000-base region of the genome corresponding to the location you clicked in Genome Graphs.
  11. Note the track at the top labeled "redbin 1" on the left and "User Supplied Track 1" above -- that is the data shown in Genome Graphs. Click on the top label "User Supplied Track 1" to make the peaks taller and more visible. (If you don't see blue peaks labeled "redbin 1" on the left, go to http://genome.ucsc.edu/cgi-bin/hgTracks?hgS_doOtherUser=submit&hgS_otherUserName=AngieHinrichs&hgS_otherUserSessionName=ABRF2010_redbin_1 to catch up.)

In general, peaks correspond to upstream regions of genes, which makes sense for RNA polymerase II ChIP-seq.

Transforming .sgr into bedGraph for Genome Browser viewing

The Genome Browser can draw a somewhat nicer display of .sgr data, if we do a simple transform of the .sgr data into a similar format called bedGraph. bedGraph has both start and end genomic coordinates, so scores can be drawn over regions as opposed to single bases. Here is a command that translates .sgr to bedGraph, expanding each data point to cover a 30-base region: 

   perl -wpe '@w=split("\t"); $_ = "$w[0]\t" . ($w[1]-1) . "\t" . ($w[1]+29) . "\t$w[2]";' \
     chr21_extended.txt_redbin.sgr > redbin.bed

If you would like to show only datapoints above a certain threshold, say 5, you can make the command slightly more complicated: 

   perl -wpe '@w=split("\t"); $_ = "$w[0]\t" . ($w[1]-1) . "\t" . ($w[1]+29) . "\t$w[2]"; \
              if ($w[2] < 15) {$_ = "";}' \
     chr21_extended.txt_redbin.sgr > redbinGt15.bed

This line can be added to the beginning of redbinGt15.bed, in order to add labels to the track display and configure default appearance (see http://genome.ucsc.edu/goldenPath/help/bedgraph.html for more information about settings). A line break is used for readability here, but must be removed if you are going to edit this track line into a file: 

   track name=redbinGt15 description="redbin items with score >= 15, expanded to 30 bases" \
     type=bedGraph viewLimits=0:100 autoScale=off maxHeightPixels=128:20:11

The file has been saved here: http://genome-test.cse.ucsc.edu/ABRF2010/redbinGt15.bed .

Steps

  1. In the Genome Browser, click the "manage custom tracks" button below the image.
  2. Click on "add custom tracks"
  3. Paste these two lines into the box labeled "Paste URLs or data" and click submit:
http://genome-test.cse.ucsc.edu/ABRF2010/redbinGt15.bed
browser position chr21:42,886,797-43,886,796
  1. You should see a confirmation page like this (image)
  2. To go to the genome browser, you can click on "Genome Browser" in the top blue bar, or click on "chr21" in the Pos column of the table showing your custom tracks, or on the "go to genome browser" button. Note: in the absence of a "browser position" line, only the top blue bar link will return you to the position you were looking at before.
  3. Now the new custom track shows grayscale representations of the score. You can click on the track title ("redbin items with score > 15") to expand the display to bar graph, which should look a lot like the "redbin 1" track created in Genome Graphs.

(If you don't see two tracks labeled "redbin 1" and "redbinGt15" over the UCSC Genes track, go to http://genome.ucsc.edu/cgi-bin/hgTracks?hgS_doOtherUser=submit&hgS_otherUserName=AngieHinrichs&hgS_otherUserSessionName=ABRF2010_redbinGt15 to catch up.)

Comparing with UCSC-hosted data

Now we will take a look at some ChIP-seq data from the ENCODE project. The "Yale TFBS" track in the Regulation section combines ChIP-seq data from several labs, and it happens to include data for RNA polymerase II ("Pol2") in HeLa-S3 cells, like Charlie Nicolet's chr21 data. We can visually compare the two datasets to see how well they correspond.

Steps

  1. In the Genome Browser, scroll down below the image to the Regulation section. Click the "+" button on the left side of the blue "Regulation" header to expand.
  2. Click on the link titled "Yale TFBS". This will take you to the track controls.
  3. In the large table under "Select subtracks by cell line and factor",
  4. Hit the "-" button in upper left corner to deselect all subtracks. Note that HeLa-S3 cell line is in 3rd column of checkboxes.
  5. Scroll down to factor Pol2. select the 3rd check box. Only HeLa/Pol2 is selected for display.
  6. Scroll back up to the top of the page. Change visibility to "full" and click submit.
  7. In the Genome Browser, HeLa/Pol2 peaks and signal appear below "Spliced ESTs".

The intervening tracks make it a little difficult to compare the ENCODE data with Charlie's, so we will adjust the display to see only the tracks we want to compare.

  1. Click the "hide all" button under the image.
  2. In the Custom Tracks section below the image, change the visibilities of redbin 1 and redbinGt15 to "full".
  3. In the Genes and Gene Prediciton Tracks section, change the visibility of UCSC Genes to "dense".
  4. In the Regulation section, change the visibility of Yale TFBS to "full".
  5. Click the "refresh" button on any blue section header.

(If you don't see a clear correspondence between the track on the bottom and the tracks on the top, go to http://genome.ucsc.edu/cgi-bin/hgTracks?hgS_doOtherUser=submit&hgS_otherUserName=AngieHinrichs&hgS_otherUserSessionName=ABRF2010_TFBS to catch up.)

Viewing alignments in BAM format

Charlie Nicolet has shared the ELAND export file that accompanies chr21_extended.txt, chr21_export.txt. This can be transformed into BAM for viewing in the UCSC Genome Browser. The samtools package includes a script that converts ELAND export files into SAM (the text version of BAM). Export files often have fasta names instead of sequence names, e.g. "chr21.fa" instead of "chr21", so the ".fa" must be removed: 

   export2sam.pl --read1=chr21_export.txt \
   | perl -wpe 's/(chr.*)\.fa/$1/' \
     > chr21.sam

Edit chr21.sam to add these header lines at the beginning (note: the words must be separated by tab characters, not spaces): 

@HD	VN:1.0	GO:none	SO:coordinate
@SQ	SN:chr21	LN:46944323

Use samtools to convert SAM to BAM, sort the BAM, and index the sorted BAM (see http://genome.ucsc.edu/goldenPath/help/bam.html for more info): 

   samtools view -bSo chr21.unsorted.bam chr21.sam
   samtools sort chr21.unsorted.bam chr21
   samtools index chr21.bam

Now the sorted chr21.bam and its index file chr21.bam.bai must be put on a web server (ftp, http, or https). For example, if you create files my.bam and my.bam.bai, and you have a personal web directory http://my.edu/mylab/~me/, you can copy or link the files to http://my.edu/mylab/~me/my.bam and http://mylab/~me/my.bam.bai. Then you would use only the first URL (http://my.edu/mylab/~me/my.bam) in a track definition line like the one below. For this example, we will use http://genome-test.cse.ucsc.edu/ABRF2010/chr21.bam .

Steps

  1. In the Genome Browser, click the "manage custom tracks" buttom below the image.
  2. Click the "add custom tracks" button.
  3. Paste this text into the box:
track name=chr21_export type=bam bigDataUrl=http://genome-test.cse.ucsc.edu/ABRF2010/chr21.bam visibility=pack
4. Click submit. Click "Genome Browser" in the blue top bar.

Now you should see an almost-solid dark blue bar with a few dark red stripes, labeled "chr21_export" above and to the left. (If you don't, got to http://genome.ucsc.edu/cgi-bin/hgTracks?hgS_doOtherUser=submit&hgS_otherUserName=AngieHinrichs&hgS_otherUserSessionName=ABRF2010_chr21_export_1MB to catch up.) There are so many alignments in this 1MB window that the items can't be drawn individually -- the image would become too tall. Let's view a smaller genomic region to get a better look at the alignments.

  1. Click and drag on the base position track to select a window around the couple largest peaks. (If it doesn't work for you, paste chr21:43,169,278-43,191,176 into the position box and click "jump".)
  2. Now we are viewing about 22,000 bases covering two peaks, and there are still so many alignments that the track is forced into "squish" display mode: items are shown at half-height and read names are suppressed. Let's zoom in again, to the right side of the left peak (or jump to position chr21:43,172,778-43,173,354).
  3. Even when viewing 577 bases, there is still quite a pileup of read alignments to view. The read names take up a lot of room, so let's hide them. Click on the light blue bar to the left of the chr21_export track, to go to track controls.
  4. There is a checkbox to the right of "Display read names" -- click to clear the checkbox, and then click submit.
  5. Mismatches appear as light red tickmarks. They are a bit easier to see if the items are drawn in grayscale, so again click on the light blue bar on the left to return to track controls.
  6. Under "additional coloring modes", select "Use gray for [alignment quality]", and click submit.
  7. Note the trend for mismatches to appear more often in alignments with lower quality, and also the relative frequency of lower-quality reads at the edge of the peak. (If you don't see that, go to http://genome.ucsc.edu/cgi-bin/hgTracks?hgS_doOtherUser=submit&hgS_otherUserName=AngieHinrichs&hgS_otherUserSessionName=ABRF2010_chr21_export_grayQual to catch up.)
  8. Zoom in once again to focus on the region with the most mismatching items (chr21:43,172,834-43,172,931). Now we see the exact mismatches between reference and read bases.
  9. Let's take a look at base qualities. Once again, click the light blue bar on the left to return to track controls.
  10. Change the selection to "Use gray for [base qualities]", and click submit. Some of the mismatches, but not all, almost disappear due to the light shades used for low quality bases. (If you don't see that, go to http://genome.ucsc.edu/cgi-bin/hgTracks?hgS_doOtherUser=submit&hgS_otherUserName=AngieHinrichs&hgS_otherUserSessionName=ABRF2010_chr21_export_baseQual).

In some regions, stacks of red tickmarks (mismatches) indicate HeLa variants with respect to the reference genome assembly. Here are several examples.

  1. Paste chr21:37,366,724-37,366,992 into the position box and click "jump".
  2. To compare with variants in dbSNP, scroll down to the Variation section and change the visibility of SNPs (130) to pack. Click "refresh" on any blue section header.
  3. The SNP track appears at the bottom of the image. SNPs classified "coding non-synonymous" by dbSNP are colored red. (If you don't see that, click http://genome.ucsc.edu/cgi-bin/hgTracks?hgS_doOtherUser=submit&hgS_otherUserName=AngieHinrichs&hgS_otherUserSessionName=ABRF2010_chr21_export_SNPs to catch up.)
  4. To see a variant not in dbSNP, paste chr21:37,367,831-37,367,845.
  5. To see a heterozygous variant adjacent to a dbSNP 3-base indel, paste chr21:37,367,478-37,367,492.
  6. This 110-base region contains three know variations from the reference assembly: chr21:37,366,860-37,366,969.

Where to learn more

Retrieved from "https://genomewiki.ucsc.edu/index.php?title=ABRF2010_Tutorial&oldid=9747"