View source for Big Bed/Wig QA

== Overview of bigBed/bigWig tracks in the Browser ==

For more information on these two file types see the following links:
* [http://genome.ucsc.edu/goldenPath/help/bigBed.html bigBed]
* [http://genome.ucsc.edu/goldenPath/help/bigWig.html bigWig]

There are two different types of bigBed/bigWig tracks in the Genome Browser; those that have BBI Link tables and those that don't. Those with BBI Link tables will often have a single-line table that points to the bigBed/bigWig file. For example, the "tgpPhase1AccessibilityPilotCriteria" table in hg19:
<pre>
mysql> select * from tgpPhase1AccessibilityPilotCriteria;
+-----------------------------------------------------------+
| fileName                                                  |
+-----------------------------------------------------------+
| /gbdb/hg19/1000Genomes/paired.end.mapping.1000G..pilot.bb |
+-----------------------------------------------------------+
</pre>

Newer bigBed/bigWig tracks often don't even have an associated BBI link table. For example, the "tgpPhase3AccessibilityPilotCriteria" track for hg19 doesn't have an associated table:
<pre>
mysql> show tables like "%tgpPhase3AccessibilityPilotCriteria%";
Empty set (0.04 sec)
</pre>
But rather has a "bigDataUrl" line in trackDb/human/hg19/trackDb.ra that points to the /gbdb/ locaiton of the file:
<pre>
    track tgpPhase3AccessibilityPilotCriteria  
    subTrack tgpPhase3Accessibility  
    shortLabel 1000G Accs Pilot  
    longLabel 1000 Genomes Project Phase 3 Paired-end Accessible Regions - Pilot Criteria  
    bigDataUrl /gbdb/hg19/1000Genomes/phase3/20141020.pilot_mask.whole_genome.bb
</pre>

More background on phasing out BBI Link tables can be found in Redmine [http://redmine.soe.ucsc.edu/issues/15059 #15059].

== QA-ing bigBed/bigWig tracks ==

Many of our standard QA scripts (countPerchrom.csh, qaGbTracks, etc.) either don't work or don't really do much when they encounter one of these bigBed/bigWig tracks.

However, you can use some Genome Browser command-line utilities to extract some information from these files and do some basic QA of these files.

===validateFiles===
First, you can use <tt>validateFiles</tt> to perform some basic validation of the bigBed/bigWig files. Here's the basic usage for <tt>validateFiles</tt>:
<pre>
validateFiles -chromInfo=FILE -options -type=FILE_TYPE file1 [file2 [...]]
</pre>
You can replace <tt>-chromInfo=FILE</tt> with <tt>-chromDb=$db</tt> to tell validate files to use the chromInfo table rather than the chrom.sizes. 
It's particularly nice to be able to just specify the db rather than having to specify the path to the chrom.sizes file (often in /hive/data/genomes/$db/chrom.sizes).

You will also need to specify the file type using <tt>-type=TYPE</tt>. For bigBed/bigWig here are the relevant sections from the usage message:
<pre>
       bigWig       : Big Wig
                      See http://genome.ucsc.edu/goldenPath/help/bigWig.html

       bigBedN[+P]  : bigBED N  or bigBED N+ or bigBED N+P, similar to BED
                      See http://genome.ucsc.edu/goldenPath/help/bigBed.html
</pre>

The "[+P]" for bigBed files will be the extra fields past the standard N BED files in your track. For example, the file /gbdb/hg19/bbi/patBulk.bb is a bigBed12+13. This means the first 12 fields are the standard BED fields and the last 13 are specialized fields that are specified in an "AutoSQL" or ".as" file. (For more on AutoSQL, see the [http://genomewiki.ucsc.edu/index.php/AutoSql AutoSQL Genomewiki Page].)

In addition to the required options described above, here are some options that are useful when QA-ing bigBed/bigWig files:
<pre>
   -as=fields.as                If you have extra "bedPlus" fields, it's great to put a definition
                                of each field in a row in AutoSql format here. Applies to bed-related types.

   -tab                         If set, expect fields to be tab separated, normally
                                expects white space separator. Applies to bed-related types.

   -doReport                    Output report in filename.report
</pre>

Note that the report files produced by <tt>-doReport</tt> are placed in the same directory as the bigBed/bigWig file. For example, if you ran <tt>validateFiles</tt> with the <tt>-doReport</tt> option on /gbdb/hg19/bbi/patNonBulk.bb, then the report files would end up in /gbdb/hg19/bbi/.

You can run <tt>validateFiles</tt> without any arguments to see the full usage message.

Here's an example of a command you might use to validate a bigBed file:
<pre>
validateFiles -chromDb=hg19 -type=bigBed12+13 -as=/hive/data/genomes/hg19/bed/patents/patSummary.as -tab /gbdb/hg19/bbi/patNonBulk.bb
</pre>

The program will report any errors it finds to stderr. If no problems are found, you should see a message like:
<pre>
Error count 0
</pre>

===bigBedInfo and bigWigInfo===

Both <tt>bigBedInfo</tt> and <tt>bigWigInfo</tt> perform similar functions for their respective file types, which is to provide an overview of the data in the bigBed or bigWig file. 
While <tt>bigBedInfo</tt> is primarily dicussed in this section, <tt>bigWigInfo</tt> and it's options and output should be fairly similar.

To start, we'll take a look at the usage statement for <tt>bigBedInfo</tt> and then discuss some of the more useful options for QA:
<pre>
bigBedInfo - Show information about a bigBed file.
usage:
   bigBedInfo file.bb
options:
   -udcDir=/dir/to/cache - place to put cache for remote bigBed/bigWigs
   -chroms - list all chromosomes and their sizes
   -zooms - list all zoom levels and their sizes
   -as - get autoSql spec
   -extraIndex - list all the extra indexes
</pre>
<tt>bigWigInfo</tt> has a similar usage message.

The <tt>-chroms</tt> option is useful in that it will print a list of all chromosomes names found in the bigBed/bigWig file. 
You can compare this list to the full list of chromosomes for that database to see which chromosomes might be missing data. 
If you're not sure whether or not a chromosome should have data on it, ask the engineer responsible for the track/table/file.

Here's an example of a command you could use to get some basic summary info on a bigBed file:
<pre>
$ bigBedInfo -chroms /gbdb/hg19/bbi/patNonBulk.bb
</pre>

And here's the output of that command:
<pre>
version: 4
fieldCount: 25
hasHeaderExtension: yes
isCompressed: yes
isSwapped: 0
extraIndexCount: 0
itemCount: 684,876
primaryDataSize: 134,443,250
zoomLevels: 0
chromCount: 24
        chr1 0 249250621
        chr10 1 135534747
        chr11 2 135006516
        chr12 3 133851895
        chr13 4 115169878
        chr14 5 107349540
        chr15 6 102531392
        chr16 7 90354753
        chr17 8 81195210
        chr18 9 78077248
        chr19 10 59128983
        chr2 11 243199373
        chr20 12 63025520
        chr21 13 48129895
        chr22 14 51304566
        chr3 15 198022430
        chr4 16 191154276
        chr5 17 180915260
        chr6 18 171115067
        chr7 19 159138663
        chr8 20 146364022
        chr9 21 141213431
        chrX 22 155270560
        chrY 23 59373566
basesCovered: 1,152,938,050
meanDepth (of bases covered): 20.639954
minDepth: 1.000000
maxDepth: 9068.000000
std of depth: 186.314036
</pre>

The <tt>itemCount</tt> line indicates the number of items in the bigBed file, similar to the number of lines in a table. It may be useful to compare this count with that of the previous version of the file or table. 

The <tt>basesCovered</tt>, as the name suggests, is the number of bases in the genome covered by items in the track. 
You can divide this number by the number of bases in the genome (found on the hgGateway "Sequences" page) to get an estimate of coverage, similar to featureBits. 
So in this case, our coverage would be <tt>(1,152,938,050/3,137,161,264)*100 = 36.75%</tt>.

The <tt>meanDepth</tt>, <tt>minDepth</tt>, <tt>maxDepth</tt>, and <tt>std of depth</tt> lines convey information about the depth of coverage for items in the bigBed file. 
For, example, the <tt>meanDepth</tt> above would indicate that, on average, every base covered by elements in the bigBed file is covered by 20 items.

===Replicating featureBits, countPerChrom, and checkTableCoords measurements===
Since most of the QA tools don't do much upon encountering bigBed files, you can convert your files to regular bed files and replicate the measurements yourself like so:
====featureBits====
If your track is a bigBed based track, then <tt>featureBits</tt> will not work. First you will need to turn your bigBed into a bed file, then, if necessary, account for tabs v. spaces, and lastly, if working with a bed12(+), split your bed into exons before running through <tt>featureBits</tt>:
<pre>
$ bigBedToBed bed.bb out.bed
$ # if necessary:
$ # sed -i 's/ /_/g' out.bed
$ # split into exons (must be bed12 or less, change cut command if you have a bed4+ or 8+, etc)
$ cut --complement -f13- out.bed | bedToExons stdin out.exons.bed
$ featureBits [db] -countGaps out.exons.bed [gap]
</pre>
If your track is a composite or super track of many bigBeds, you can run the above steps in a loop like so:
<pre>
$ for file in ../*.bed; do printf "%s\n" "$file"; sed 's/ /_/g' $file | \
cut --complement -f13- | bedToExons stdin $file.exons.bed; \
featureBits hg38 -countGaps $file.exons.bed; \
featureBits hg38 -countGaps $file.exons.bed gap; done
../unipChain.bed
33610384 bases of 3209286105 (1.047%) in intersection
0 bases of 3209286105 (0.000%) in intersection
</pre>
====countPerChrom====
<tt>countPerChrom.csh</tt> runs on MySQL tables, and uses the separate utility <tt>graph.csh</tt> to actually make the nice histograms. You can run <tt>graph.csh</tt> manually on a two column file of chromosome counts computed from a grep like so:
<pre>
for chr in {1..22} X Y M; do printf "chr%s %s\n" "$chr" "$(grep -wc "^chr${chr}" unipChain.bed)"; done > hg38.unipChain.chromCounts
graph.csh hg38.unipChain.chromCounts 35
 chr1 xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
 chr2 xxxxxxxxxxxxxxxxxxxxx              
 chr3 xxxxxxxxxxxxxxxxxx                 
 chr4 xxxxxxxxxxxxx       
...
</pre>
If your track is a composite or super track of many bigBeds, you can run the above steps in a loop like so:
<pre>
$ for file in *.bed; do for chr in {1..22} X Y M; do printf "chr%s %s\n" "$chr" "$(grep -wc "^chr${chr}" $file)"; done > hg38.$file.chromCounts; done
$ for file in *.chromCounts; do graph.csh $file 35 > $file.graph; echo "each x=35" >> $file.graph; done
$ cat hg38.unipChain.bed.chromCounts.graph 

 chr1 xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
 chr2 xxxxxxxxxxxxxxxxxxxxx              
 chr3 xxxxxxxxxxxxxxxxxx                 
 chr4 xxxxxxxxxxxxx       
...
</pre>
====positionalTblCheck and checkTableCoords====
<tt>positionalTblCheck</tt> checks that tables are sorted via chromosome and chromStart. <tt>bedToBigBed</tt> does this check before running, so if a bigBed exists it will have to be sorted. Similarly, <tt>bedToBigBed</tt> also checks that all the coordinates in the input bed file are valid for the chrom.sizes file provided, so the <tt>checkTableCoords</tt> check can be skipped as well.

== Other Tools That May Be Useful for QA of big* tracks ==

* <tt>bigWigToBedGraph</tt> or <tt>bigWigToWig</tt>— these programs convert a bigWig file to ASCII bedGraph/WIG format.
* <tt>bigWigSummary</tt> — this program extracts summary information from a bigWig file.
* <tt>bigBedToBed</tt> — this program converts a bigBed file to ASCII BED format.
* <tt>bigBedSummary</tt> — this program extracts summary information from a bigBed file. 



[[Category:Browser QA]]