Usage

dspchip command line options

dspchip comes with a number of options to perform signal analysis of a ChIP-seq experiments. These values apply to dspchip 0.8.2

• -i,--sigA
• -c,--sigB
• --formatA
• --formatB
• -n,--name
• --pl
• --winfun
• --fir
• --threshold
• --profile
• -W,--wavelet
• -C,--corr
• --scale
• -e,--ewin
• -w,--wstep
• -s,--profstep
• -r,--ratio
• --ichrom
• --echrom
• --allchr
• --csize
• -p,--peaks
• --modelfeat
• --pvalue
• --dist
• -K
• -q,--quality
• --downA
• --downB
• --nodup
• --nosig
• --noprofile
• --save
• --stats

-i,--sigA

The first signal (A) to analyze, typically the ChIP file

-c,--sigB

The second signal (B) to analyze, typically a control

--formatA

The file format for signal A. Accepted format are bam, bw (for bigwig), bed, bar or wig. Default is bam.

--formatB

The file format for signal B. If not specified it is set equal to format A

-n,--name

The name of the experiment. This will be the prefix for all the output files. Default is default

--pl

The analysis pipeline. Pipeline steps is a string made of the following letters: F W G A E N C L S R X V T Z. Each letters mean a step:

• F: Perform FIR window smoothing. A FIR window should be specified (see --fir option)
• W: Perform discrete wavelet transform for smoothing. (see --wavelet option)
• G: Apply a gaussian blur smoothing
• A: Perform signal autocorrelation
• E: Equalize the signal so that cumulative distribution of value is (0,1)
• N: Normalize the signal values so that it is distributed with N(0, 1)
• L: Calculates the log2ratio between signal A and signal B
• S: Subtracts signal B from signal A
• X: Cross-correlate two signals (deprecated)
• V: Convolve two signals (deprecated)
• R: Calculates the magnitude of the variation (i.e. (s1 - s2) / s2)
• T: Calculates the threshold correct data accordingly (Z step is highly suggested after this)
• Z: Removes negative values

Note that L, S, X, V and R reduce the number of signals to one. Default is NFTZ

--winfun

The windowing function for profile output (that is bedgraph). This can be mean, med or max (for average, median or max). Default is mean.

--fir

The FIR window name. This can be flat, hanning, hamming, bartlett, blackman or gauss. Default is flat.

--threshold

Choose between thresholding approach. Thresholding happens at peak finding time. When thresholding is engaged, peaks won't be searched using canny-like edge detector. Accepted values are otsu and mad. This value can be specified with a keyword global or local to perform thresholding on the whole chromosome or local windows; in order to use this, the keyword should be separated by a comma (i.e. local,otsu). Default value global,otsu.

--profile

The file format for profile output. Can be bdg (for gzipped bedgraph) or bw (for bigwig). Default is bw

-W,--wavelet

The name of the wavelet to use. This has to be defined in pywt package (check here). Default is db5

-C,--corr

Whether to perform correlation analysis. This will plot a chart of correlation between two signals

--scale

Looks for original data range (min, max) and rescales processed data to a similar interval. Useful to avoid meaningless values such as, say, 0.00001231.

-e,--ewin

The expected size of the feature under investigation. Default 1000

-w,--wstep

The step size for windowing functions. Default 100

-s,--profstep

The step size for output profile. Default 100

-r,--ratio

The ratio between the minimum value and the maximum value in a peak to consider it a single peak or two separate peaks. This essentially checks the depth of the "valley" separating peaks. Default 0.4

--ichrom

A comma separated list of chromosomes to be included in the analysis

--echrom

A comma separated list of chromosomes to be excluded from the analysis

--allchr

Whether to consider all the chromosomes in file. This included _random chromosomes which are usually excluded _

--csize

When bed, bar or wig file is specified, a chromosome size table is needed to allocate space for each chromosome. bigwig and bam files already contain a chromosome size index. The file has to be in the format

chrN N

where chrN is any chromosome and N is its size. If you have 2bit files from UCSC genome browser you can get the size just issuing

$ twoBitInfo genomefile.2bit chrlen.csv

where genomefile.2bit is unique to each genome build (e.g. hg18.2bit, mm9.2bit)

-p,--peaks

Whether to call peaks. This is not set by default.

--modelfeat

Specify the feature to model for p-value estimation. Can be area, height or length. The peak finding will be looking at that feature to build value distribution. Default is area

--pvalue

The p-value filter for peaks. Default 1.0

--dist

Tries to infer the proper distribution fitting on data. Distributions evaluated are: Normal, LogNormal, Gamma, PowerLaw. Most likely will be a Gamma distribution Default hist.

-K

In order to define a max window to define separate peaks, this parameters will set the max distance as K ** expWindow. Default is 0.5. **

-q,--quality

When input data is BAM, filter alignments with quality less than the specified threshold. Default is 0.

--downA

The ratio for downsample data from signal A. This is useful when data are not balanced. Downsample is performed using an uniform distribution. Works only with interval files (bam, bed). Default is 1.0.

--downB

The ratio for downsample data from signal B.

--nodup

A naive deduplication filter for sorted interval files (bam, bed).

--nosig

Do not perform signal analysis.

--noprofile

Do not output a profile output. Useful when looking for peaks only.

--save

Save numpy arrays for each step. For debugging purposes.

--stats

Save data statistics. For debugging purposes. This is likely to be removed in the future.

[Anonymous]

Originally posted by: Vladimir...@gmail.com

My code to get the chromosome sizes:

import urllib2

UCSC_chrom_URL = r'http://hgdownload.cse.ucsc.edu/goldenPath/%s/database/chromInfo.txt.gz'

def get_chrom_length ( dbname ):
    """dbname is the database name in UCSC, like hg18, mm9...
    """
    # first try to find a local file called dbname
    try:
        fhd = open(dbname,"r")
        chrom_len = {}
        for l in fhd:
            fs = l.split()
            try:
                chrom_len[fs[0]] = int(fs[1])
            except:
                pass
        fhd.close()
    except:
        # get chromosome length from UCSC db download page
        f = urllib2.urlopen(UCSC_chrom_URL % (dbname))
        tmpfname = tempfile.mkstemp(prefix="qcmanychip")[1]
        tmpf = open(tmpfname,'w')
        tmpf.write(f.read())                # write file content to temp file
        tmpf.close()
        f.close
        # read it
        fhd = gzip.open(tmpfname,'r')
        chrom_len = {}
        for l in fhd:
            fs = l.split()
            try:
                chrom_len[fs[0]] = int(fs[1])
            except:
                pass
        fhd.close()
        os.unlink(tmpfname)
    return chrom_len

-Tao