Tag Archives: downloads

Accessing the Genome Browser Programmatically Part 2 – Using the Public MySQL Server and gbdb System

If you missed part 1 about obtaining sequence data, you can catch up here. Note: We now have an API which can also perform many of these functions.

The UCSC Genome Browser is a large repository of data from multiple sources, and if you want to query that annotation data, the easiest way to get started is via the Table Browser. Choose the assembly and track of interest and click the “describe table schema” button, which will show the MySQL database name, the primary table name, the fields of the table and their descriptions. If the track is stored not in MySQL but as a binary file (like bigBed or bigWig) in /gbdb, it will show a file name, e.g. "Big Bed File: /gbdb/dm6/ncbiRefSeq/ncbiRefSeqOther.bb". If this is the case, skip directly to the Accessing the gbdb directory system section below. Otherwise, the track data is either a single MySQL table or a set of related tables, which you can either download as gzipped text files from the “Annotation Database” section on our downloads page (here’s the GRCh37/hg19 listing) and work on them locally, or use the public MySQL server and issue MySQL queries remotely. Generally speaking, the format for most of our tables is similar to the formats described here, e.g., in bed (“chrom chromStart chromEnd”) format, and we do not store any sequence or contigs in our databases, which means you’ll need to use the instructions in Part 1 of this blog series in order to get any raw sequence data.

Accessing the public MySQL server
The best way to showcase the public MySQL server is to show some examples — here are a few to get you started:
1. If you want to download some transcripts from the new NCBI RefSeq Genes track, you can use the following command:

$ mysql -h genome-mysql.soe.ucsc.edu -ugenome -A -e "select * from ncbiRefSeq limit 2" hg38
+-----+-------------+-------+--------+---------+-------+----------+--------+-----------+--------------------------------------------------------------------+--------------------------------------------------------------------+-------+---------+--------------+------------+-----------------------------------+
| bin | name        | chrom | strand | txStart | txEnd | cdsStart | cdsEnd | exonCount | exonStarts                                                         | exonEnds                                                           | score | name2   | cdsStartStat | cdsEndStat | exonFrames                        |
+-----+-------------+-------+--------+---------+-------+----------+--------+-----------+--------------------------------------------------------------------+--------------------------------------------------------------------+-------+---------+--------------+------------+-----------------------------------+
| 585 | NR_046018.2 | chr1  | +      |   11873 | 14409 |    14409 |  14409 |         3 | 11873,12612,13220,                                                 | 12227,12721,14409,                                                 |     0 | DDX11L1 | none         | none       | -1,-1,-1,                         |
| 585 | NR_024540.1 | chr1  | -      |   14361 | 29370 |    29370 |  29370 |        11 | 14361,14969,15795,16606,16857,17232,17605,17914,18267,24737,29320, | 14829,15038,15947,16765,17055,17368,17742,18061,18366,24891,29370, |     0 | WASH7P  | none         | none       | -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1, |
+-----+-------------+-------+--------+---------+-------+----------+--------+-----------+--------------------------------------------------------------------+--------------------------------------------------------------------+-------+---------+--------------+------------+-----------------------------------+

2. If you are interested in a particular enhancer region, for instance “chr1:166,167,154-166,167,602”, and want to find the nearest genes within a 10kb range, then the following query will do the job:

$ chrom="chr1"
$ chromStart="166167154"
$ chromEnd="166167602"
$ mysql -h genome-mysql.soe.ucsc.edu -ugenome -A -e "select \
   e.chrom, e.txStart, e.txEnd, e.strand, e.name, j.name as geneSymbol from ncbiRefSeqCurated e,\
   ncbiRefSeqLink j where e.name = j.id AND e.chrom='${chrom}' AND \
      ((e.txStart >= ${chromStart} - 10000 AND e.txStart <= ${chromEnd} + 10000) OR \ (e.txEnd >= ${chromStart} - 10000 AND e.txEnd <= ${chromEnd} + 10000)) \
order by e.txEnd desc " hg38
+-------+-----------+-----------+--------+----------------+------------+
| chrom | txStart   | txEnd     | strand | name           | geneSymbol |
+-------+-----------+-----------+--------+----------------+------------+
| chr1  | 166055917 | 166166755 | -      | NR_135199.1    | FAM78B     |
| chr1  | 166055917 | 166166755 | -      | NM_001320302.1 | FAM78B     |
| chr1  | 166069298 | 166166755 | -      | NM_001017961.4 | FAM78B     |
+-------+-----------+-----------+--------+----------------+------------+

3. If you need to get gene names and their lengths for RNA-seq read normalization, you can use the following query:

$ mysql -h genome-mysql.soe.ucsc.edu -u genome -A -e “ \
  select l.name, kr.value, psl.qEnd - psl.qStart as length \
  from   refGene r, hgFixed.refLink l, knownToRefSeq kr, knownCanonical kc, refSeqAli psl \
  where  r.name = l.mrnaAcc and r.name = kr.value and kr.name = kc.transcript \
         and r.name = psl.qName group by kr.value limit 3” hg38
+-------+-----------+--------+
| name  | value     | length |
+-------+-----------+--------+
| A2M   | NM_000014 |   4920 |
| NAT2  | NM_000015 |   1317 |
| ACADM | NM_000016 |   2622 |
+-------+-----------+--------+

In addition to our download site and public MySQL server hosted here in California, we have also recently added support for a download site (http://hgdownload-euro.soe.ucsc.edu) and public MySQL server (genome-euro-mysql.soe.ucsc.edu) hosted in Europe, which will speed up downloads for many of our users.

Please follow the Conditions for Use when querying the public MySQL servers.

Many of the command line utilities available on our utilities downloads server are also able to interact with our databases or download files, like mafFetch (as long as your ~/.hg.conf file is present as discussed below):

$ mafFetch xenTro9 multiz11way region.bed stdout
##maf version=1
##maf version=1 scoring=blastz
a score=0.000000
s xenTro9.chr9     15946024 497 +  80437102 ACTAT...
e galGal5.chr14     1678315   0 -  15595052 I
e xenLae2.chr9_10L 13130032 2034 - 117834370 I

a score=2992.000000
s xenTro9.chr9     15946521 145 +  80437102 TCATC...
s xenLae2.chr9_10L 13132066 148 - 117834370 TTATC...

Note: Only the first 5 bases on each line and only the first 10 lines are shown for brevity.

Here we are directly querying the mutliz11way table for the Xenopus tropicalis xenTro9 assembly, no need to download the entire alignment file to the local disk and query manually. Commands of this nature usually require a special private .hg.conf file in the user’s home directory (note the leading dot). This configuration file contains a couple key=value lines that most of our programs can parse and then use to access the public MySQL server. This page contains a sample .hg.conf file that can be used by most of the command line utilities to direct them to access either our US MySQL server or our European MySQL server. That sample .hg.conf is certainly enough to get started, but for more information about the various Genome Browser configuration options, please see the comments in the ex.hg.conf and minimal.hg.conf files.


Accessing the gbdb directory system
The third method of grabbing our data is via the /gbdb/ directory system. This location, browsable here, holds most of the bigBed, bigWig, and other large data files that we do not keep directly in MySQL databases/tables. There are many utilities available for manipulating these files, and most of them are able to work on remote files, for example:

$ bigBedToBed -chrom=chr1 -start=5563837 -end=5564370 http://hgdownload.soe.ucsc.edu/gbdb/hg38/crispr/crispr.bb stdout 
chr1    5563870    5563893        55    +    5563870    5563890    0,200,0    255,255,0    128,128,0    CAAGTGGAATCAGGATGCCT    GGG    55    72% (57)    52% (46)    10    60    MIT Spec. Score: 55, Doench 2016: 72%, Moreno-Mateos: 52%    3345002138
chr1    5563878    5563901        59    +    5563878    5563898    0,200,0    0,200,0    128,128,0    ATCAGGATGCCTGGGATATG    TGG    59    63% (54)    61% (50)    6    63    MIT Spec. Score: 59, Doench 2016: 63%, Moreno-Mateos: 61%    22777603204

Also note that we have all of this data available via rsync as well, so the following command will work to download the crispr.bb file referenced above:

$ rsync -vh hgdownload.soe.ucsc.edu::gbdb/hg38/crispr/crispr.bb
-rw-rw-r--  1466266135 2017/03/30 14:31:48 crispr.bb

sent 33 bytes  received 70 bytes  206.00 bytes/sec
total size is 1.47G  speedup is 14235593.54

If you are interested in say, Human GRCh37/hg19 gbdb data, then all you have to do is change the “hg38” at the end of the template http://hgdwonload.soe.ucsc.edu/gbdb/hg38 url to “hg19”, resulting in http://hgdwonload.soe.ucsc.edu/gbdb/hg19. This holds for all databases at UCSC, like mm10 or bosTau8.

Summary
Just as in part 1, if you are going to continually request parts of the same files or table over and over again, it is best to download the file from our downloads server and operate on it locally. All of our track data, including MySQL tables and bigBed/Wig/BAM files are hosted on our downloads server at http://hgdownload.soe.ucsc.edu. Generally speaking bigBeds/bigWigs/BAMs and other binary files are located in the hgdownload.soe.ucsc.edu/gbdb/ location discussed earlier, while MySQL table data in gzipped plain text format can be found at http://hgdownload.soe.ucsc.edu/goldenPath/$db (where $db is a database name like hg19 or hg38) or via queries against the public MySQL server directly.

Stay tuned for part 3 of this programmatic access series — controlling the Genome Browser image!

If after reading this blog post you have any public questions, please email genome@soe.ucsc.edu. All messages sent to that address are archived on a publicly accessible forum. If your question includes sensitive data, you may send it instead to genome-www@soe.ucsc.edu.

Accessing the Genome Browser Programmatically Part 1 – How to get sequence from the UCSC Genome Browser

Note: We now have an API which can also perform many of these functions.

As the number of bioinformaticians have grown since the inception of the UCSC Genome Browser in 2000, there has been an increased need for programmatic access to the data and tools hosted at UCSC. Although there is no true API developed by UCSC (yet), there are a number of ways to interface with the UCSC Genome Browser, some more efficient than others. The intention of this blog post series is to explain some of the preferred ways to access the commonly requested Genome Browser data and tools and to add a bit of explanation of the architecture of the UCSC Genome Browser in general. The three most common requests are 1) how to download a single stretch of sequence in FASTA format, 2) how to download multiple ranges of sequence, and 3) how to get basic statistics on the nucleotides in a sequence. If you want the in-depth examples and explanation, skip down, but if you’re crunched for time, all you really need to know is the following three Q&As:

Q: How do I extract some sequence?
A: The best choice is to use the twoBitToFa command, available for your system here (Windows 10 users can use the linux.x86_64/ binaries in the Windows Subsystem for Linux). Here’s an example:

$ twoBitToFa http://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.2bit:chr1:100100-100200 stdout
>chr1:100100-100200
gcctagtacagactctccctgcagatgaaattatatgggatgctaaatta
taatgagaacaatgtttggtgagccaaaactacaacaagggaagctaatt

Q: What if I have a list of coordinates?
A: Again use twoBitToFa, this time with the -bed option (also check out the post on coordinate systems):

$ cat input.bed
chr1 4150100 4150200 seq1
chr1 4150300 4150400 seq2
$ twoBitToFa http://hgdownload.soe.ucsc.edu/goldenPath/mm10/bigZips/mm10.2bit -bed=input.bed stdout
>seq1
gcatcccagtcctgatactggaaaattcatttagtgacaagcgagggcca
cttgggattctctcacccccatatttaggagaccttattagggtcacctt
>seq2
tatccccttccctccccaccagatactacaattcacatcatactctgtcc
cccagtctacccataaaatctattctatttacctctccaaacgaagatct

Q: How do I count A, C, G, T?
A: twoBitToFa followed by faCount (available from the same location as twoBitToFa):

$ twoBitToFa http://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.2bit:chr1:100100-100200 stdout | faCount stdin
#seq    len     A       C       G       T       N       cpg
chr1:100100-100200      100     37      17      21      25      0       0
total   100     37      17      21      25      0       0

Run twoBitToFa or faCount with no arguments to get a usage message and view all of their options:

$ faCount
faCount - count base statistics and CpGs in FA files.
...


The most efficient way to get sequence from UCSC Genome Browser

The most common data request we receive is a request for FASTA sequence or sequences, making it a fitting subject for part 1 of this blog series about programmatic access to the Genome Browser. If you are browsing a region in the genome browser and you want to get a FASTA sequence for just the region you are browsing, using the keyboard shortcut ‘vd’ (v then d for view DNA) is probably the easiest way. But what about when you want to get sequences for a list of regions? What about if you need your web application to download the sequence? You could download sequence interactively with the Table Browser, although the solution is somewhat cumbersome: first you must make a custom track of the region(s) you would like sequence for, and then use the “output format: sequence” option with your custom track selected as the primary track. Fortunately, there is a much easier approach – downloading the 2bit file for your organism of interest and then using the twoBitToFa command on it like so:

$ wget http://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.2bit
$ twoBitToFa hg38.2bit:chr1:100100-100200 stdout
>chr1:100100-100200
gcctagtacagactctccctgcagatgaaattatatgggatgctaaatta
taatgagaacaatgtttggtgagccaaaactacaacaagggaagctaatt

The twoBitToFa command is available from the list of public utilities, in the directory appropriate to your operating system. twoBitToFa even accepts a URL to our downloads server as the 2bit argument, so if you wanted to grab some mm10 sequence, or even a list of sequences, you can just query the downloads server directly like so:

$ cat input.bed
chr1 4150100 4150200 seq1
chr1 4150300 4150400 seq2
$ twoBitToFa http://hgdownload.soe.ucsc.edu/goldenPath/mm10/bigZips/mm10.2bit -bed=input.bed stdout
>seq1
gcatcccagtcctgatactggaaaattcatttagtgacaagcgagggcca
cttgggattctctcacccccatatttaggagaccttattagggtcacctt
>seq2
tatccccttccctccccaccagatactacaattcacatcatactctgtcc
cccagtctacccataaaatctattctatttacctctccaaacgaagatct

Note that “stdout” in the above commands is a special option (along with the corresponding “stdin”) that tells the majority of UCSC commands to read/write from/to /dev/stdin and /dev/stdout instead of the required filenames, and is exemplified by the following common usage of generating some quick statistics on a region like chr1:100100-100200:

$ twoBitToFa http://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.2bit:chr1:100100-100200 stdout | faCount stdin
#seq    len     A       C       G       T       N       cpg
chr1:100100-100200      100     37      17      21      25      0       0
total   100     37      17      21      25      0       0

The twoBitToFa and URL to hgdownload 2bit combo is important because our downloads server is significantly more robust than our DAS CGI, can support more requests, and won’t slow the main site down for other users. We’ve also noticed that our DAS server often receives many requests for the same sequence, so for those of you providing software where the same query will be made multiple times, consider whether it would be more efficient to download an entire 2bit file to your local disk, rather than send the same query thousands of times to our servers.

Summary
twoBitToFa and faCount are two useful utilities, among the many other hundreds of tools available, that are useful for extracting sequence data. While not as preferable to working with locally downloaded files, twoBitToFa can also work with URLs to 2bit files, such as those on the UCSC Genome Browser download site. Stay tuned for part 2 of this programmatic access series — Using the Genome Browser public MySQL server and gbdb.

If after reading this blog post you have any public questions, please email genome@soe.ucsc.edu. All messages sent to that address are archived on a publicly accessible forum. If your question includes sensitive data, you may send it instead to genome-www@soe.ucsc.edu.