docs/README.chip-seq

   1 This is an updated version of the core of the ChIP-seq
   2 analysis code described in Johnson et al (2007).  It
   3 should run on any Unix-like system supporting python 2.5
   4 or better. The code is developed on Linux and MacOS X on
   5 python 2.5.
   6
   7 These scripts in the ChIPSeqMini package are now part of
   8 the ERANGE package, but are still available as a
   9 standalone package for now.
  10
  11 This code is made available as open-source, as described
  12 in the copyright file ERANGE.COPYRIGHT.
  13
  14
  15 1. REQUIREMENTS
  16 2. COMMAND LINE OPTIONS
  17 3. MAKING THE NECESSARY INPUT (RDS) FILES
  18 4. WEIGHING MULTIREADS
  19 5. RUNNING THE PEAK FINDER
  20 6. DISPLAYING DATA ONTO THE  UCSC GENOME BROWSER
  21 7. DOWNSTREAM ANALYSES
  22
  23
  24 1. REQUIREMENTS
  25
  26 1) Python 2.5 is required because some of the scripts and
  27 Cistematic (see below) need pysqlite, which is now bundled in
  28 Python.
  29
  30 2) You will also need to use Cistematic 2.3 (available at
  31 cistematic.caltech.edu) for all of the scripts that are
  32 part of the downstream analyses.
  33
  34 (optional) Use of the psyco module (psyco.sf.net) on 32-bit
  35 Linux or Mac Intel machines is highly recommended.
  36
  37 (optional) Three visualization scripts also depend on the
  38 additional package pylab (matplotlib). These scripts are:
  39 - getgosig.py
  40 - plotbardist.py
  41 - scatterfields.py
  42 You do not need to install pylab if you will be
  43 visualizing some of your analysis results differently.
  44
  45
  46 2. COMMAND LINE OPTIONS
  47
  48 You can find out more about the settings for each script
  49 by typing:
  50
  51 python $ERANGEPATH/<scriptname>
  52
  53 to see the command line options, where ERANGEPATH is the
  54 environmental variable set to the path to the directory
  55 holding the ERANGE scripts. Note that the command line
  56 options are case sensitive and that they could well
  57 fail silently.
  58
  59
  60 3. MAKING THE NECESSARY INPUT (RDS) FILES
  61
  62 You will want to first convert your read mappings to the
  63 native ERANGE read store. Please see the file
  64 README.build-rds for instructions on how to do this.
  65
  66 Build an RDS file for both the ChIP, and if available and
  67 appropriate, the control. Note that we *HIGHLY* recommend
  68 the use of a matched control sample to account for some
  69 of the general background artifacts that can be present
  70 in ChIP-seq samples (e.g. DNAse hypersensitivity,
  71 assembly collapse of some sattelite repeats, etc....).
  72
  73
  74 4. WEIGHING MULTIREADS
  75
  76 Version 3.0 of the peak finder can use multireads, i.e.
  77 reads that map equally well to more than one location
  78 in the genome, to find binding sites that are in low
  79 copy-number on-unique regions (typically less than 10).
  80
  81 ERANGE offers 3 ways to analyze these regions:
  82 (a) default weighing of 1/multiplicity
  83 (b) ignoring multireads
  84 (c) weighing of multireads based on unique reads in a
  85 given radius
  86
  87 (a) is the default in the current release of ERANGE.
  88 Simply proceed to RUNNING THE PEAK FINDER for (a) and
  89 (a). You can ignore multireads (b) by using the --nomulti
  90 flag with findall.py. For (c), use weighMultireads.py
  91 to weigh multireads based on a unique reads in the
  92 respective radius of each potential location. Once run,
  93 proceed to the section below.
  94
  95
  96 5. RUNNING THE PEAK FINDER
  97
  98 To run the peak finder without read shifting, use the
  99 following command:
 100
 101 python $ERANGEPATH/findall.py label chip.rds chip.regions.txt --control control.rds --listPeak --revbackground
 102
 103 which will run the peak finder on chip.rds / control.rds ,
 104 store the enriched region coordinates in chip.regions.txt,
 105 also store the actual local maximum in each region in the
 106 same file, and also calculate an FDR by running the
 107 finder on control.rds / chip.rds .
 108
 109 A log file (findall.log by default, change with -log)
 110 tracks the settings used to run the program as well as
 111 some of the summary statistics, which are also stored
 112 at the bottom of the regions.txt output file.
 113
 114 findall.py is tuned to conservative settings for 10-12M
 115 mappable read IPs of static, sequence-specific
 116 transcription factors in mammals with very short
 117 fragment sizes, on the order of 40-60 bp.
 118
 119 You will *NEED* to change some of the default parameters
 120 if working in smaller genomes (e.g. use smaller -spacing),
 121 if working with certain types of IPs such as histones and
 122 polymerases (test with and without --notrim and
 123 --nodirectionality), if working with rather weak IPs
 124 (e.g. --minimum and --ratio), or if working with larger
 125 fragment sizes (see the paragraph below discussing read
 126 shifting).
 127
 128 findall.py returns a per-peak p-value. By default, this
 129 is calculated using a Poisson distribution of peak RPMs
 130 (or counts, if using --raw) for each chromosome in the IP.
 131 P-value calculations can be turned off using
 132 '--pvalue none '. Alternatively, the p-value can be
 133 calculated from the background using the option
 134 '--pvalue back ', which must be combined with the option
 135 --revbackground.
 136
 137 By default, findall.py does not try to adjust the location
 138 of the reads based on half the size of the expected fragment
 139 length (the "shift"). If you believe that you need to shift
 140 your peaks, findall.py can try to pick the best shift based
 141 on the best shift for strong sites using the parameter
 142 '--shift learn '. You can also either manually specify a
 143 shift value using '--shift #bp ' or ou can calculate a
 144 "best shift" for each region using '--autoshift'. If you
 145 need to using the shift options, the recommended usage is:
 146 (i) first run findall.py with '--shift learn ', which will
 147 peak a shift if there are at least 30 regions that meet
 148 its training criteria.
 149 (ii) if (i) couldn't pick a shift, run findall.py with
 150 --autoshift and --reportshift
 151 (iii) look at the mode (most common #) for the shift
 152 (iv) rerun findall.py with --shift #bp where #bp is the mode
 153
 154 If you are storing the RDS files on an network-mounted
 155 directory, make sure to use '--cache XXXXX' to enable
 156 local caching, where is as large as appropriate as
 157 described in section 9 of README.build-rds .
 158
 159 Note that ERANGE will cache by default to /tmp, but this
 160 can be redirected to any directory pointed to by the
 161 environmental variable CISTEMATIC_TEMP.
 162
 163 To find out the current default settings and options,
 164 simply type:
 165
 166 python $ERANGEPATH/findall.py
 167
 168 for more information.
 169
 170
 171 6. DISPLAYING DATA ONTO THE  UCSC GENOME BROWSER
 172
 173 You can output bed-files of the raw reads using
 174 makebedfromrds.py and  BEDGRAPH file using
 175 makewiggle.py as described in README.build-rds .
 176
 177 You can create bed files of regions and sites (see
 178 below) using regiontobed.py and makesitetrack.py .
 179
 180
 181 7. DOWNSTREAM ANALYSES
 182
 183 Recall that Cistematic 2.3 is a required to do motif
 184 and gene-level analyses of the output of findall.py.
 185
 186 Use getallgenes.py to find the nearest gene within a
 187 radius of each binding site.
 188
 189 Use analyzego.py to do a Gene Ontology enrichment
 190 analysis of a gene list (such as from getallgenes.py).
 191 You can look at a heatmap of your GO enrichments using
 192 getgosig.py. You can also use getGOgenes.py to look at
 193 the genes with particular GO annotations.
 194
 195 To do motif-finding, use getfasta.py to get the sequences
 196 centered on the peaks of your regions of interest. For
 197 the sake of a pleasant experience, try limiting yourself
 198 to less than 100kb of combined sequence (the easiest being
 199 by picking your regions with the strongest signals).
 200
 201 Once you have a fasta file of the regions of interest, you
 202 can use findMotifs.py to find motifs using either
 203 cisGreedy (bundled with Cistematic 2.2) which is good for
 204 shorter motifs or Meme (must be installed separately -
 205 refer to the instructions on cistematic.caltech.edu for
 206 more information), which is better for longer motifs.
 207 findmotifs.py will return a set motifs in Cistematic format
 208 with a .mot extension. These motifs can then be used with
 209 getallsites.py to get the coordinates and instances of each
 210 motif in all of the regions found by the peak finder.
 211
 212 The sites can be checked against repeat-masker annotations
 213 (preloaded from UCSC with buildrmaskdb.py) using
 214 checkrmask.py. The sites for each motif can also be fed
 215 back into getallgenes.py to get genes, redo the GO analyses,
 216 etc....
 217
 218 You can use the intersect scripts (intersects.py,
 219 gointersects.py, and siteintersects.py) to compare different
 220 sets of genes/GO/site results across multiple experiments,
 221 for example.
 222
 223
 224 RELEASE HISTORY
 225
 226 version 3.2    November 2010 - updated command line options
 227 version 3.1    February 2009 - support for read shifting
 228 version 3.0    February 2009 - support for UCSC narrowPeak format in regiontobed.py
 229 version 3.0rc1 December 2008 - added parameter to control peak-trimming
 230 version 3.0b2  December 2008 - added per-peak p-value
 231 version 3.0b   November 2008 - initial release of RDS-based code
 232 with support for eland and bowtie.
 233