Be careful of "sort -k1,1 -k2,2n -u"

When you attempted to sort and extract the unique genomic regions using "sort -k1,1 -k2,2n -u", you might make a mistake by missing the region with the same chr and start, but different end position.

The right way should be  "sort -k1,1 -k2,2n -k3,3n -u" or  "sort -k1,1 -k2,2n | sort -u"

How to extract the gap region in human genome?

Just notice that I should avoid the gap region, esp. when we generate a random background as your null distribution using tools such as bedtools shuffle.

Short answer: go below UCSC Table Browser link and choose to save as a bed file
 http://genome.ucsc.edu/cgi-bin/hgTables?clade=mammal&org=Human&db=hg19&hgta_group=allTables&hgta_track=hg19&hgta_table=gap&hgta_regionType=genome&hgta_outputType=primaryTable

As below table shown, 8.28% of hg19 assembly are simply gap.

Gap (gap) Summary Statistics

	item count 457 item bases 239,845,127 (8.28%) item total 239,845,127 (8.28%) smallest item 47 average item 524,825 biggest item 30,000,000

altColor in UCSC track setting

Many track types allow setting a color range that varies from color to altColor. For instance the CpG Island tracks use the altColor setting to display the weaker islands, while the stronger ones are rendered in color. If altColor is not specified, the system will use a color halfway between that specified in the color tag and white instead.

Be aware that wiggles with negative values are drawn in altColor not color as positive values are.

Hierarchical structure of UCSC Genome Browser track hub

UCSC Genome Browser tracks can be organized into groups by using the container multiWig, compositeTrack on, and superTrack on lines. Supertracks can contain composite tracks and container multiWigs, but not vice versa. With supertracks, composite tracks, and container multiWigs, children will inherit the settings from their parents, but can override their parent settings within their own stanzas.

Here is their hierarchical relationship:

superTrack on

|==== child tracks (bam, bigBed, bigWig or vcfTabix)

|==== container multiWig

     |---- child track (bigWig)

|==== compositeTrack on

     |==== view

         |---- child track with subGroups setting (bam, bigBed, bigWig or vcfTabix, but not mix)

         |---- child track with subGroups setting (bam, bigBed, bigWig or vcfTabix, but not mix) 

superTrack also allow to mix other supported types of hub tracks: bam, bigBed, bigWig or vcfTabix.

a bigWigSummary bug

Write down a bigwigsummary bug I found today. It's found when I attempted to get the max value in a region using -type=max and dataPoints=1:

$ bigWigSummary test.bw -type=max chr12 54070173 54072173 1

13.3672

$ bigWigSummary test.bw -type=max chr12 54070173 54072173 10

0.944904 1.02475 0.568405 0.741671 1.43119 1.08896 0.705965 0.542034 0.380971 0.591934

As you see, if I use dataPoints=1, the max value is 13.3672 and when I use dataPoints=10 the max value is 1.43119. So there must be something wrong, since the max value should not change no matter how many data points we check. Visualizing the bigwig in UCSC Genome Browser shows that 1.43119 is correct for this case. Interestingly, 13.3672 is indeed a peak summit, but not for this region, rather a region upstream. I don’t know why bigWigSummary take the summit from region outside. This only happened when dataPoints=1. 

I put the data below. You can download to test:

https://www.dropbox.com/sh/q963gxj9e1kj8ly/AAAG-LqCT_WQi6dgMm0lhoDIa?dl=0

People reported similar error for bigwigSummary, for example

bigWigSummary outputs different values as bigWigAverageOverBed: UCSC team explained "The reason for this is that the summary levels have some rounding error and some border conditions when extracting data over relatively small regions." They suggested to use bigWigAverageOverBed if you want the highest level of accuracy. But bigWigAverageOverBed won't output the max and also it's nothing with high or low accuracy, but rather a bug. 

Tao Liu also reported another bug for bigwig when it's converted from wig in compressed manner (by default), and suggested to fix it by using -unc when converting wig to bigwig. I've tried to use -unc when converting from bedGraph to bigwig, the bug is still there. 

Still looking for workout, and also report to UCSC:

https://groups.google.com/a/soe.ucsc.edu/forum/#!topic/genome/pWjcov-xQyQ

Update: One workout I found is, to use intersectBed and groupBy in bedtools on bedGraph file (rather than bigwig). Here is pseudocode:

intersectBed -a regions.bed -b signal.bedGraph -wo -sorted | groupBy -g 1,2,3 -c 7 -o max > regions.maxSignal.bed

Update2:  use the -minMax option from the latest version of bigWigAverageOverBed (from v304). See reply from UCSC group:

Thank you for contacting us. One of our engineers was able to reproduce this and says the bigWigSummary program uses the bigWig summary levels, which are lower resolution versions of the data, so when you ask for a single value in a particular range, the range that is used may include bases that are before and after the range.
If you want only the values exactly within the range, you can use (the latest version from v304) bigWigAverageOverBed like so:
echo "chr12 54070173 54072173 one" | bigWigAverageOverBed myTest.bw stdin stdout -minMax | cut -f 8

Trimmed mean and median in AWK

This can be easily done in R, but sometime you want to get it in scripting language like awk or perl in order to process the big data line by line.

Here is the code snip:

# for median
# Note: thanks to anonymous reply below, put c, j in arguments to define that they are local variables.

function median(v, c, j) 

{ 

    c=asort(v,j); 

    if (c % 2) return j[(c+1)/2]; 

    else return (j[c/2+1]+j[c/2])/2; 

}

# for trimmed mean (where p is the percentage of data to be trimmed)
function trimmedMean(v, p, c, j) 

{ 

    c=asort(v,j); 

    a=int(c*p);

    for(i=a+1;i<=(c-a);i++) s+=j[i];

    return s/(c-2*a); 

}

To use it, for example if we want to generate a merged bigwig track for a list of samples, we can take median value of all samples at each genomic position, here is it:

unionBedGraphs -i `ls *normalized.bedGraph` | awk 'function median(v) {c=asort(v,j); if (c % 2) return j[(c+1)/2]; else return (j[c/2+1]+j[c/2])/2.0; } {OFS="\t"; n=1; for(i=4;i<=NF;i++) S[n++]=$i; print $1,$2,$3, median(S)}' > AllSamples.median.normalized.bedGraph
bedGraphToBigWig AllSamples.median.normalized.bedGraph ChromInfo.txt AllSamples.median.normalized.bw

about UCSC Genome Browser track

Two points:

1. The custom track data may be compressed by any of the following programs: gzip (.gz), compress (.Z), or bzip2 (.bz2). But not for bigwig and bam.

2. In a track hub Db configuration file, up to 9 subgroup types can be defined for a composite, such as:

subGroup1 <gTag1> <gTitle1> <mTag1a=mTitle1a> [mTag1b=mTitle1b…]
subGroup2 <gTag2> <gTitle2> <mTag2a=mTitle2a> [mTag2b= mTitle2b…]
...
subGroup9 <gTag9> <gTitle9> <mTag9a=mTitle9a> [mTag9b= mTitle9b…]

But these is no such limitation (I guess so, not test yet) for the tag/title pairs in each subGroup. For example, ENCODE data trackDb put all TFs in one subGroup:
http://ftp.ebi.ac.uk/pub/databases/ensembl/encode/integration_data_jan2011/hg19/trackDb.txt

One question: How to share tracks but secure the data files in the track hub?

The current directory hierarchy for a hub is like:

myHub/ - directory containing track hub files

     hub.txt -  a short description of hub properties
     genomes.txt - list of genome assemblies included in the hub data
     hg19/ - directory of data for the hg19 (GRCh37) human assembly
          trackDb.txt - display properties for tracks in this directory
          dnase.html - description text for a DNase track 
          dnaseLiver.bigWig - wiggle plot of DNase in liver
          dnaseLiver.bigBed - regions of active DNase
          dnaseLung.bigWig - wiggle plot of DNase in lung
          dnaseLung.bigWig - regions of active DNase
          ...
          rnaSeq.html - description text for an RNAseq track
          rnaSeqLiver.bigWig - wiggle plot of RNAseq data in liver
          rnaSeqLiver.bigBed - intron/exon lists for liver
          rnaSeqLung.bigWig - wiggle plot of RNAseq data in lung
          rnaSeqLung.bigBed - intron/exon lists for lung
     hg18/ - directory of data for the hg18 (Build 36) human assembly
          trackDb.txt - display properties for tracks in this directory
          dnase.html - description text for a DNase track 
          dnaseLiver.bigWig - wiggle plot of DNase data in liver
          dnaseLiver.bigBed - regions of active DNase
          dnaseLung.bigWig - wiggle plot of DNase data in lung
          dnaseLung.bigWig - regions of active DNase
          ...
          rnaSeq.html - description text for an RNAseq track
          rnaSeqLiver.bigWig - wiggle plot of RNAseq data in liver
          rnaSeqLiver.bigBed - intron/exon lists for liver
          rnaSeqLung.bigWig - wiggle plot of RNAseq data in lung
          rnaSeqLung.bigBed - intron/exon lists for lung

The UCSC webpage also indicates that "unlisted hubs are in no way secure." But this is definitely a unsolved problem. Maybe the only solution is to set up your own local mirror?

How to use the "textHistogram -aveCol" and textHist2 in the UCSC executables

The Jim Kent's executable tools are general quite friendly to use. But for some options of some tools, it's a bit too terse. For example, the "-aveCol" option in textHistogram says

-aveCol=N - A second column to average over. The averages
will be output in place of counts of primary column.

Also, for textHist2, it says

textHist2 - Make two dimensional histogram table out

of a list of 2-D points, one per line.

usage:

   textHist2 input

options:

   -xBins=N - number of bins in x dimension

   -yBins=N - number of bins in y dimension

   -xBinSize=N - size of bins in x dimension

   -yBinSize=N - size of bins in x dimension

   -xMin=N - minimum x number to record

   -yMin=N - minimum y number to record

   -ps=output.ps - make PostScript output

   -psSize=N - Size in points (1/72th of inch)

   -labelStep=N - How many bins to skip between labels

   -margin=N - Margin in points for PostScript output

   -log    - Logarithmic output (only works with ps now)

   -postScale=N (default 1.000000) - What to scale by after normalization

What do this mean? For example, I have userFile with 2 columns like,
1 8
2 9
3 6
5 3

Here is reply from Brooke Rhead in UCSC Genome Bioinformatics Group:
The -aveCol option prints the average values of the items in each bin instead of the number of items in each bin. (The default bin size is 1.) So, if you run textHistogram on your example input without the -aveCol option and the default bin size, you get the number of items in each bin, which is 1 item for all bins except #4:

$ textHistogram userFile
1 ************************************************************ 1
2 ************************************************************ 1
3 ************************************************************ 1
4 0
5 ************************************************************ 1

If you instead use -aveCol=2, you get the average of the *value* in column 2. Since bin size is still equal to 1, this amounts to printing the value in column 2 in each bin:

$ textHistogram -aveCol=2 userFile
1 ***************************************************** 8.000000
2 ************************************************************ 9.00000
3 **************************************** 6.000000
4 0.000000
5 ******************** 3.000000

If you specify a smaller bin size, the option might seem more useful:

$ textHistogram -binSize=2 -aveCol=2 userFile
0 ************************************************************ 8.00000
2 ******************************************************** 7.500000
4 *********************** 3.000000

The bins are taken from the first column: 0-1, 2-3, and 4-5, and the values are from the second column. So, for instance, the bin labeled "2" above contains the average of the values from this part of your file:
2 9
3 6
. . . bin "2" contains the average of 9 and 6, or 7.5.


The textHist2 program makes more sense if you specify the number of bins and bin sizes to use. For example, if we specify the following options on your same input file:

1 8
2 9
3 6
5 3

we get:

$ textHist2 -xBins=6 -yBins=10 userFile
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 0 0 1
0 0 0 0 0 0
0 0 0 0 0 0
0 0 0 1 0 0
0 0 0 0 0 0
0 1 0 0 0 0
0 0 1 0 0 0

The first column in the input file specifies the x-axis in the two dimensional histogram, and the second column specifies the y-axis. The numbering of the rows and columns starts at 0 in both cases. The first line of your input:

1 8

Causes a 1 to appear in the second column from the right and the 9th row from the bottom, or column 1 and row 8, if you start counting from zero. The "2 9" line causes the 1 to appear in column 2, row 9, and so on.

grouping multiWig as subtracks

This is not a tip, but a request (to a UCSC genome feature).

I was trying to set multiWig containers as a subtrack of a composite (as below), but it does not work.


track CAGE
compositeTrack on
shortLabel CAGE
longLabel CAGE Tracks
type bed 3
visibility full
priority 20
subGroup1 view Views SGuniq=Signal_Uniqmap SGall=Signal_Allmap
subGroup2 stage Stage 10dpp=10dpp 20dpp=20dpp
subGroup3 treatment treatment wt=Wild_Type
dimensions dimensionX=treatment dimensionY=stage
sortOrder stage=+ treatment=+ view=+
dragAndDrop subTracks
configurable on

track CageUniqmapsignal
shortLabel Signal
view SGuniq
type bigWig
autoScale on
alwaysZero on
windowingFunction maximum
maxHeightPixels 50:50:11
parent CAGE

# --------- 10dpp ---------
track 10dpp_wt_CAGE_uniq
container multiWig
shortLabel CAGE_uniq
longLabel 10dpp_wt_CAGE_uniq
aggregate transparentOverlay
showSubtrackColorOnUi on
subGroups view=SGuniq stage=10dpp treatment=wt
parent CAGE

track 10dpp_wt_CAGE_uniq.plus
bigDataUrl http://10dpp_wt_CAGE_uniq+.bw
parent 10dpp_wt_CAGE_uniq
color 0,0,255 

track 10dpp_wt_CAGE_uniq.minus
bigDataUrl http://10dpp_wt_CAGE_uniq-.bw
parent 10dpp_wt_CAGE_uniq
color 255,0,0
# --------- tracks for other stages... ---------


Just found that UCSC has not implemented it yet. Here is the message:

"It is not currently possible to combine bigBeds and bigWigs into a multiWig
track. They can be combined into a multiple view composite. The bigWigs
could be combined into a multiWig, but that multiWig cannot itself be part
of a composite. The bigWigs could be displayed both in a multiWig and as
subtracks of a separate composite that also contains bigWigs.

In the future there may be a "container" like multiWig that will contain
both bigBeds and bigWigs. Separately, in the future we may support a
multiWig as one subtrack of a composite. However, neither of these is
currently supported."

Hope this can be implemented soon!!!

(If anyone has a solution to the problem, please let me know. Thanks very much)

convert bed to gtf, gtf to bed

bed --> gtf:

bedToGenePred input.bed stdout | genePredToGtf file stdin output.gtf

Here is an example from UCSC wiki (http://genomewiki.ucsc.edu/index.php/Genes_in_gtf_or_gff_format):

Some gene tracks are in a bed format in the database, perhaps with extra columns past the standard bed format. In this case, extract the standard bed columns, convert it to a genePred and then to a gtf. For example

mysql --user=genome --host=genome-mysql.cse.ucsc.edu -A -N -e "select chrom,chromStart,chromEnd,name,score,strand,thickStart,thickEnd from wgRna;" hg19 | bedToGenePred stdin stdout | genePredToGtf file stdin wgRna.gtf

gtf --> bed:

download gtf2bed.pl from Eric's site:

http://code.google.com/p/ea-utils/source/browse/trunk/clipper/gtf2bed

perl gtf2bed.pl input.gtf > output.bed

If simply converting gtf to bed6 format (no exons/intron info), the following bash line can work:

fgrep -w transcript gencode.v17.annotation.gtf | sed 's/[";]//g;' | awk '{OFS="\t"; print $1, $4-4,$5,$12,0,$7}'

or a bed9 format with additional info on gene/type/name

fgrep -w transcript gencode.v17.annotation.gtf | sed 's/[";]//g;' | awk '{OFS="\t"; print $1, $4-4,$5,$12,0,$7,$18,$14,$10}' 

Also, I edited Eric's code a bit to output a more meaningful gene ID, e,g.

($transid) = $f[8]=~ /transcript_id "([^"]+)"/;
($geneid) = $f[8]=~ /gene_id "([^"]+)"/;
($gene_type) = $f[8]=~ /gene_type "([^"]+)"/;
($gene_name) = $f[8]=~ /gene_name "([^"]+)"/;
($trans_type) = $f[8]=~ /transcript_type "([^"]+)"/;
$id="${gene_name}__${geneid}__${transid}__${gene_type}.${trans_type}";

How to get tRNA/rRNA/mitochondrial gene GTF file

Cufflinks/Tophat ask for a GTF file to mask the abundant transcripts (e.g. tRNA/rRNA/chrM). Here is the step to get such a file:

Go to UCSC Table browser:

For tRNA/rRNA:

    * Select "All Tables" from the group drop-down list
    * Select the "rmsk" table from the table drop-down list
    * Choose "GTF" as the output format
    * Type a filename (e.g. "rRNA.tRNA.gtf") in "output file" so your browser downloads the
result
    * Click "create" next to filter
    * Next to "repClass," type rRNA
    * Next to free-form query, select "OR" and type repClass = "tRNA"
    * Click submit on that page, then get output on the main page

For chrM genes:

    * Select "All Tables" from the group drop-down list
    * Select the "knownGene" table from the table drop-down list
    * Choose "GTF" as the output format
    * Type a filename (e.g. "chrM.gtf") in "output file" so your browser downloads the result
    * Click "create" next to filter
    * Next to "chrom," type chrM
    * Click submit on that page, then get output on the main page

cat the two files from above steps. That's it. 
P.S.  I just noticed that for galgal3 (chicken), UCSC does not have rRNA/tRNA in the repeat database. 

How to tell which library type to use (fr-firststrand or fr-secondstrand)?

First of all, as a bioinformatian, you should ask the data producer (e.g. the one who prepared the RNAseq library) which protocol they used to generate the data.

Tophat manual page has listed the general strand-specific protocol:

Library Type	Examples	Description
fr-unstranded	Standard Illumina	Reads from the left-most end of the fragment (in transcript coordinates) map to the transcript strand, and the right-most end maps to the opposite strand.
fr-firststrand	dUTP, NSR, NNSR	Same as above except we enforce the rule that the right-most end of the fragment (in transcript coordinates) is the first sequenced (or only sequenced for single-end reads). Equivalently, it is assumed that only the strand generated during first strand synthesis is sequenced.
fr-secondstrand	Ligation, Standard SOLiD	Same as above except we enforce the rule that the left-most end of the fragment (in transcript coordinates) is the first sequenced (or only sequenced for single-end reads). Equivalently, it is assumed that only the strand generated during second strand synthesis is sequenced.

In case you don't know the library-type, you can still figure it out by yourself. Tophat FAQ page provided a solution for that (http://tophat.cbcb.umd.edu/faq.html#library_type). But more simply (comparing to running 1M reads first), you can choose few reads and BLAT to genome and infer the library-type from the mapping result.

Generally, reads from the left-most end of RNA fragment (always from 5´ to 3´) are always mapped to transcript-strand, and (for pair-end sequencing) reads from the right-most end are always mapped to the opposite strand. See the arrows direction in the below schema. This is because the sequencer always read from 5´ to 3´.

Summary of library type protocols (for Tophat/Bowtie)

But regarding to which strand the RNA fragment is synthesized from, this involves different strand-specific protocols. Thanks to the illustration figure (see below) from Zhao Zhang, we could see that for example dUTP method is to only sequence the strand from the first strand synthesis (the original RNA strand is  degradated due to the dUTP incorporated), so the /2 read is from the original RNA strand.

Strand-specific library protocols (Credit: Zhao Zhang)

Taking a real example, first getting some reads (in fasta format) from the paired-end sequencing fastq file using command like:

$ zcat ~/nearline/rnaseq/BU/Jul2012/Sample_3576_H_01.R1.fastq.gz | sed 's/@//g;s/ /_/g' | awk '{if(NR%4==1)print ">"$0;if(NR%4==2) print $0;}' | head

$ zcat ~/nearline/rnaseq/BU/Jul2012/Sample_3576_H_01.R2.fastq.gz | sed 's/@//g;s/ /_/g' | awk '{if(NR%4==1)print ">"$0;if(NR%4==2) print $0;}' | head

Blatting them in UCSC Genome Browser

Below is screenshot for top hits of one pair of reads. They mapped to exons of OS9 genes (the left one is /1 and right one is /2, with opposite direction). We see that /1 mapped to transcript direction, /2 mapped to opposite direction, which means it can only be fr-secondstrand or fr-unstrand (cannot be fr-firststrand).

Continuing to look at other reads in the file, we can find examples like these:

where /2 mapped to transcript strand and /1 mapped to the opposite strand. Combining with the observation from above, we can conclude that this is a fr-unstrand library.

get UCSC images for a list of regions in batch

Here is my working R code for the task. It can be simplified as 3 lines.

# example of controling individual track

#theURL="http://genome.ucsc.edu/cgi-bin/hgTracks?db=mm9&wgRna=hide&cpgIslandExt=pack&ensGene=hide&mrna=hide&intronEst=hide&mgcGenes=hide&hgt.psOutput=on&cons44way=hide&snp130=hide&snpArray=hide&wgEncodeReg=hide&pix=1000&refGene=pack&knownGene=hide&rmsk=hide"

# example of controling via session

theURL="http://genome-preview.ucsc.edu/cgi-bin/hgTracks?hgS_doOtherUser=submit&hgS_otherUserName=sterding&hgS_otherUserSessionName=mm9_piRNA&hgt.psOutput=on&pix=1000"

# read regions

toPlot=read.table("../data/piRNA.clusters.coordinates.cpg.bed", header=T)

## paralle version

# library(multicore)

# mclapply(1:nrow(toPlot), function(i) screenshotUCSC(theURL, "", as.character(toPlot$chr[i]), toPlot$start[i]-2000, toPlot$end[i]+1999, paste("region_", i, "_", toPlot$name[i],".pdf", sep="")), mc.cores=10)

# anti-robot version 

# UCSC Policy: Program-driven use of this software is limited to a maximum of one hit every 15 seconds and no more than 5,000 hits per day.

for(i in 1:nrow(toPlot)){

    screenshotUCSC(theURL, "", as.character(toPlot$chr[i]), toPlot$start[i]-2000, toPlot$end[i]+1999, paste("region_", i, "_", toPlot$name[i],".pdf", sep=""))

    Sys.sleep(5) 

}

# merge script

mergePDF("piRNAs_ucsc_screeshot.pdf", list.files(pattern="region_.*.pdf"))

try(system("rm region_.*.pdf"))

####### lib ############

# Here is an R script wrote by Aaron Statham which saves UCSC to pdfs -

# you can choose which genome and tracks to display by altering the 'url' parameter. 'trackfile' is the url of a file describing the custom tracks (beds/bigwigs) to display

mergePDF <- function(output="merged.pdf", sourcefiles=c("source1.pdf","source2.pdf","source3.pdf"))

{

    # create the command string and call the command using system()

    # merging command from http://hints.macworld.com/article.php?story=2003083122212228

    command=paste("gs -q -dNOPAUSE -dBATCH -sDEVICE=pdfwrite",paste("-sOutputFile",output, sep="="), paste(sourcefiles, collapse=" "),sep=" ")

    try(system(command))

}

#  Reference: http://www.biostars.org/post/show/6132/batch-viewing-of-ucsc-browser-graphic/

screenshotUCSC <- function(url, trackfile, chr, start, end, filename) {

        oldpen <- options("scipen")

        options(scipen=100)

        temp <- readLines(paste(url, "&hgt.customText=", trackfile, "&position=",chr,":",start,"-",end, sep=""))

        #cat(temp,"\n")

        pdfurl <- paste("http://genome-preview.ucsc.edu/trash",gsub(".*trash","",gsub(".pdf.*","",temp[grep(".pdf", temp, fixed=TRUE)][1])), ".pdf", sep="")

        cat(pdfurl,"\n");

        options(scipen=oldpen)

        download.file(pdfurl, filename, mode="wb", quiet=TRUE)

}

Three ways to convert bam/bed file to bigwig, separated by strand

Here are three ways to convert bam/bed to bigwig, separated by strand:


# -----------  method 1


bamToBed -i accepted_hits.bam -split > accepted_hits.bed


awk '{if($6=="+") print}' accepted_hits.bed | sort -k1,1 | bedItemOverlapCount mm9 -chromSize=ChromInfo.txt stdin | sort -k1,1 -k2,2n > accepted_hits.plus.bedGraph
awk '{if($6=="-") print}' accepted_hits.bed | sort -k1,1 | bedItemOverlapCount mm9 -chromSize=ChromInfo.txt stdin | sort -k1,1 -k2,2n | awk '{OFS="\t"; print $1,$2,$3,"-"$4}' > accepted_hits.minus.bedGraph

bedGraphToBigWig accepted_hits.plus.bedGraph ChromInfo.txt accepted_hits.plus.bw

bedGraphToBigWig accepted_hits.minus.bedGraph ChromInfo.txt accepted_hits.minus.bw

# ----------- method 2

bamToBed -i accepted_hits.bam -split > accepted_hits.bed


sort -k1,1 accepted_hits.bed | awk -v '{print $0 >> "accepted_hits.bed"$6}'
bedItemOverlapCount $index -chromSize=ChromInfo.txt accepted_hits.bed+ | sort -k1,1 -k2,2n > accepted_hits.plus.bedGraph
bedItemOverlapCount $index -chromSize=ChromInfo.txt accepted_hits.bed- | sort -k1,1 -k2,2n | awk '{OFS="\t"; print $1,$2,$3,"-"$4}' > accepted_hits.minus.bedGraph

bedGraphToBigWig accepted_hits.plus.bedGraph ChromInfo.txt accepted_hits.plus.bw

bedGraphToBigWig accepted_hits.minus.bedGraph ChromInfo.txt accepted_hits.minus.bw

# ----------- method 3

bedtools genomecov -ibam -bg -split -strand + -i accepted_hits.bam -g ChromInfo.txt > accepted_hits.plus.bedGraph
bedtools genomecov -ibam -bg -split -strand - -i accepted_hits.bam -g ChromInfo.txt > accepted_hits.minus.bedGraph

bedGraphToBigWig accepted_hits.plus.bedGraph ChromInfo.txt accepted_hits.plus.bw

bedGraphToBigWig accepted_hits.minus.bedGraph ChromInfo.txt accepted_hits.minus.bw

bigWigSummary prefers absolute path

I just got this weird error when I was running bigwigsummary on a bigwig file - it keeps saying:


Can't open ~/scratch/bw/wgEncodeOpenChromChipMcf7Pol2Aln_2Reps.norm5.rawsignal.bw to read
No such file or directory

But obviously the file is there with permission of 755 (-rwxr-xr-x).

Finally by taking Brian's suggestion, I used absolute path, instead of relative path with tilde (e.g. ~ ), the error disappeared.

Don't know if this is required by the code itself or just an usual difference between ~ and /home.

get rRNA.gtf file from UCSC Table Browser

You can get this pretty easily from the UCSC table browser (http://genome.ucsc.edu/cgi-bin/hgTables).

Select "All Tables" from the group drop-down list
Select the "rmsk" table from the table drop-down list
Choose "GTF" as the output format
Type a filename in "output file" so your browser downloads the result
Click "create" next to filter
Next to "repClass," type rRNA
Next to free-form query, select "OR" and type repClass = "tRNA"
Click submit on that page, then get output on the main page