Browse files

Add option --SPMR to save fragment pileup per million reads to

bedGraph files.
  • Loading branch information...
1 parent f189c87 commit 21951a2a2e1048da7f09a4148b1f74cacfd36cf3 @taoliu committed May 2, 2012
Showing with 10,788 additions and 11,096 deletions.
  1. +4,579 −5,155 MACS2/IO/cBedGraph.c
  2. +10 −3 MACS2/IO/cBedGraph.pyx
  3. +95 −122 MACS2/IO/cFixWidthTrack.c
  4. +2 −2 MACS2/IO/cFixWidthTrack.pyx
  5. +2,731 −2,565 MACS2/IO/cScoreTrack.c
  6. +43 −30 MACS2/IO/cScoreTrack.pyx
  7. +4 −1 MACS2/OptValidator.py
  8. +3,271 −3,185 MACS2/cPeakDetect.c
  9. +43 −26 MACS2/cPeakDetect.pyx
  10. +10 −7 bin/macs2
View
9,734 MACS2/IO/cBedGraph.c
4,579 additions, 5,155 deletions not shown because the diff is too large. Please use a local Git client to view these changes.
View
13 MACS2/IO/cBedGraph.pyx
@@ -1,5 +1,5 @@
# cython: profile=True
-# Time-stamp: <2012-05-01 18:19:56 Tao Liu>
+# Time-stamp: <2012-05-02 16:46:19 Tao Liu>
"""Module for Feature IO classes.
@@ -856,9 +856,16 @@ class bedGraphTrackI:
return ret
- def make_scoreTrack_for_macs (self, bdgTrack2 ):
+ def make_scoreTrack_for_macs (self, bdgTrack2, float effective_depth_in_million = 1.0 ):
"""A modified overlie function for MACS v2.
+ effective_depth_in_million: sequencing depth in million after
+ duplicates being filtered. If
+ treatment is scaled down to
+ control sample size, then this
+ should be control sample size in
+ million. And vice versa.
+
Return value is a bedGraphTrackI object.
"""
cdef int pre_p, p1, p2
@@ -867,7 +874,7 @@ class bedGraphTrackI:
assert isinstance(bdgTrack2,bedGraphTrackI), "bdgTrack2 is not a bedGraphTrackI object"
- ret = scoreTrackI()
+ ret = scoreTrackI( effective_depth_in_million = effective_depth_in_million )
retadd = ret.add
chr1 = set(self.get_chr_names())
View
217 MACS2/IO/cFixWidthTrack.c
@@ -1,4 +1,4 @@
-/* Generated by Cython 0.14.1 on Tue May 1 21:57:39 2012 */
+/* Generated by Cython 0.14.1 on Wed May 2 14:11:24 2012 */
#define PY_SSIZE_T_CLEAN
#include "Python.h"
@@ -994,7 +994,6 @@ static char __pyx_k__length[] = "length";
static char __pyx_k__maxnum[] = "maxnum";
static char __pyx_k__random[] = "random";
static char __pyx_k__resize[] = "resize";
-static char __pyx_k__sample[] = "sample";
static char __pyx_k__stdout[] = "stdout";
static char __pyx_k__strand[] = "strand";
static char __pyx_k____doc__[] = "__doc__";
@@ -1035,7 +1034,6 @@ static char __pyx_k____version__[] = "__version__";
static char __pyx_k__RuntimeError[] = "RuntimeError";
static char __pyx_k__print_to_bed[] = "print_to_bed";
static char __pyx_k__get_chr_names[] = "get_chr_names";
-static char __pyx_k__random_sample[] = "random_sample";
static char __pyx_k__sample_percent[] = "sample_percent";
static PyObject *__pyx_kp_s_10;
static PyObject *__pyx_kp_s_11;
@@ -1110,13 +1108,11 @@ static PyObject *__pyx_n_s__obj;
static PyObject *__pyx_n_s__percent;
static PyObject *__pyx_n_s__print_to_bed;
static PyObject *__pyx_n_s__random;
-static PyObject *__pyx_n_s__random_sample;
static PyObject *__pyx_n_s__range;
static PyObject *__pyx_n_s__readonly;
static PyObject *__pyx_n_s__refcheck;
static PyObject *__pyx_n_s__resize;
static PyObject *__pyx_n_s__round;
-static PyObject *__pyx_n_s__sample;
static PyObject *__pyx_n_s__sample_num;
static PyObject *__pyx_n_s__sample_percent;
static PyObject *__pyx_n_s__samplesize;
@@ -6951,13 +6947,11 @@ static __Pyx_StringTabEntry __pyx_string_tab[] = {
{&__pyx_n_s__percent, __pyx_k__percent, sizeof(__pyx_k__percent), 0, 0, 1, 1},
{&__pyx_n_s__print_to_bed, __pyx_k__print_to_bed, sizeof(__pyx_k__print_to_bed), 0, 0, 1, 1},
{&__pyx_n_s__random, __pyx_k__random, sizeof(__pyx_k__random), 0, 0, 1, 1},
- {&__pyx_n_s__random_sample, __pyx_k__random_sample, sizeof(__pyx_k__random_sample), 0, 0, 1, 1},
{&__pyx_n_s__range, __pyx_k__range, sizeof(__pyx_k__range), 0, 0, 1, 1},
{&__pyx_n_s__readonly, __pyx_k__readonly, sizeof(__pyx_k__readonly), 0, 0, 1, 1},
{&__pyx_n_s__refcheck, __pyx_k__refcheck, sizeof(__pyx_k__refcheck), 0, 0, 1, 1},
{&__pyx_n_s__resize, __pyx_k__resize, sizeof(__pyx_k__resize), 0, 0, 1, 1},
{&__pyx_n_s__round, __pyx_k__round, sizeof(__pyx_k__round), 0, 0, 1, 1},
- {&__pyx_n_s__sample, __pyx_k__sample, sizeof(__pyx_k__sample), 0, 0, 1, 1},
{&__pyx_n_s__sample_num, __pyx_k__sample_num, sizeof(__pyx_k__sample_num), 0, 0, 1, 1},
{&__pyx_n_s__sample_percent, __pyx_k__sample_percent, sizeof(__pyx_k__sample_percent), 0, 0, 1, 1},
{&__pyx_n_s__samplesize, __pyx_k__samplesize, sizeof(__pyx_k__samplesize), 0, 0, 1, 1},
@@ -7232,59 +7226,38 @@ PyMODINIT_FUNC PyInit_cFixWidthTrack(void)
if (PyObject_SetAttr(__pyx_m, __pyx_n_s__logging, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 21; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
- /* "MACS2/IO/cFixWidthTrack.pyx":24
- *
- * #from array import array
- * from random import sample as random_sample # <<<<<<<<<<<<<<
- * import sys
- * from copy import copy
- */
- __pyx_t_1 = PyList_New(1); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 24; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __Pyx_GOTREF(((PyObject *)__pyx_t_1));
- __Pyx_INCREF(((PyObject *)__pyx_n_s__sample));
- PyList_SET_ITEM(__pyx_t_1, 0, ((PyObject *)__pyx_n_s__sample));
- __Pyx_GIVEREF(((PyObject *)__pyx_n_s__sample));
- __pyx_t_2 = __Pyx_Import(((PyObject *)__pyx_n_s__random), ((PyObject *)__pyx_t_1)); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 24; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __Pyx_GOTREF(__pyx_t_2);
- __Pyx_DECREF(((PyObject *)__pyx_t_1)); __pyx_t_1 = 0;
- __pyx_t_1 = PyObject_GetAttr(__pyx_t_2, __pyx_n_s__sample); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 24; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __Pyx_GOTREF(__pyx_t_1);
- if (PyObject_SetAttr(__pyx_m, __pyx_n_s__random_sample, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 24; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
- __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
-
/* "MACS2/IO/cFixWidthTrack.pyx":25
* #from array import array
- * from random import sample as random_sample
+ * #from random import sample as random_sample
* import sys # <<<<<<<<<<<<<<
* from copy import copy
* from MACS2.Constants import *
*/
- __pyx_t_2 = __Pyx_Import(((PyObject *)__pyx_n_s__sys), 0); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 25; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __Pyx_GOTREF(__pyx_t_2);
- if (PyObject_SetAttr(__pyx_m, __pyx_n_s__sys, __pyx_t_2) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 25; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
+ __pyx_t_1 = __Pyx_Import(((PyObject *)__pyx_n_s__sys), 0); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 25; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_GOTREF(__pyx_t_1);
+ if (PyObject_SetAttr(__pyx_m, __pyx_n_s__sys, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 25; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
/* "MACS2/IO/cFixWidthTrack.pyx":26
- * from random import sample as random_sample
+ * #from random import sample as random_sample
* import sys
* from copy import copy # <<<<<<<<<<<<<<
* from MACS2.Constants import *
* from libc.stdint cimport uint32_t, uint64_t, int32_t, int64_t
*/
- __pyx_t_2 = PyList_New(1); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 26; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __Pyx_GOTREF(((PyObject *)__pyx_t_2));
+ __pyx_t_1 = PyList_New(1); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 26; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_GOTREF(((PyObject *)__pyx_t_1));
__Pyx_INCREF(((PyObject *)__pyx_n_s__copy));
- PyList_SET_ITEM(__pyx_t_2, 0, ((PyObject *)__pyx_n_s__copy));
+ PyList_SET_ITEM(__pyx_t_1, 0, ((PyObject *)__pyx_n_s__copy));
__Pyx_GIVEREF(((PyObject *)__pyx_n_s__copy));
- __pyx_t_1 = __Pyx_Import(((PyObject *)__pyx_n_s__copy), ((PyObject *)__pyx_t_2)); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 26; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __Pyx_GOTREF(__pyx_t_1);
- __Pyx_DECREF(((PyObject *)__pyx_t_2)); __pyx_t_2 = 0;
- __pyx_t_2 = PyObject_GetAttr(__pyx_t_1, __pyx_n_s__copy); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 26; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __pyx_t_2 = __Pyx_Import(((PyObject *)__pyx_n_s__copy), ((PyObject *)__pyx_t_1)); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 26; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_GOTREF(__pyx_t_2);
- if (PyObject_SetAttr(__pyx_m, __pyx_n_s__copy, __pyx_t_2) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 26; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
+ __Pyx_DECREF(((PyObject *)__pyx_t_1)); __pyx_t_1 = 0;
+ __pyx_t_1 = PyObject_GetAttr(__pyx_t_2, __pyx_n_s__copy); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 26; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_GOTREF(__pyx_t_1);
+ if (PyObject_SetAttr(__pyx_m, __pyx_n_s__copy, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 26; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+ __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
/* "MACS2/IO/cFixWidthTrack.pyx":27
* import sys
@@ -7293,16 +7266,16 @@ PyMODINIT_FUNC PyInit_cFixWidthTrack(void)
* from libc.stdint cimport uint32_t, uint64_t, int32_t, int64_t
* from cpython cimport bool
*/
- __pyx_t_1 = PyList_New(1); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 27; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __Pyx_GOTREF(((PyObject *)__pyx_t_1));
+ __pyx_t_2 = PyList_New(1); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 27; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_GOTREF(((PyObject *)__pyx_t_2));
__Pyx_INCREF(((PyObject *)__pyx_n_s_27));
- PyList_SET_ITEM(__pyx_t_1, 0, ((PyObject *)__pyx_n_s_27));
+ PyList_SET_ITEM(__pyx_t_2, 0, ((PyObject *)__pyx_n_s_27));
__Pyx_GIVEREF(((PyObject *)__pyx_n_s_27));
- __pyx_t_2 = __Pyx_Import(((PyObject *)__pyx_n_s_26), ((PyObject *)__pyx_t_1)); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 27; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __Pyx_GOTREF(__pyx_t_2);
- __Pyx_DECREF(((PyObject *)__pyx_t_1)); __pyx_t_1 = 0;
- if (__pyx_import_star(__pyx_t_2) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 27; __pyx_clineno = __LINE__; goto __pyx_L1_error;};
- __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
+ __pyx_t_1 = __Pyx_Import(((PyObject *)__pyx_n_s_26), ((PyObject *)__pyx_t_2)); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 27; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_GOTREF(__pyx_t_1);
+ __Pyx_DECREF(((PyObject *)__pyx_t_2)); __pyx_t_2 = 0;
+ if (__pyx_import_star(__pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 27; __pyx_clineno = __LINE__; goto __pyx_L1_error;};
+ __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
/* "MACS2/IO/cFixWidthTrack.pyx":31
* from cpython cimport bool
@@ -7311,10 +7284,10 @@ PyMODINIT_FUNC PyInit_cFixWidthTrack(void)
* cimport numpy as np
*
*/
- __pyx_t_2 = __Pyx_Import(((PyObject *)__pyx_n_s__numpy), 0); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 31; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __Pyx_GOTREF(__pyx_t_2);
- if (PyObject_SetAttr(__pyx_m, __pyx_n_s__np, __pyx_t_2) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 31; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
+ __pyx_t_1 = __Pyx_Import(((PyObject *)__pyx_n_s__numpy), 0); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 31; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_GOTREF(__pyx_t_1);
+ if (PyObject_SetAttr(__pyx_m, __pyx_n_s__np, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 31; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
/* "MACS2/IO/cFixWidthTrack.pyx":37
* # constants
@@ -7350,8 +7323,8 @@ PyMODINIT_FUNC PyInit_cFixWidthTrack(void)
* """Fixed Width Locations Track class II along the whole genome
* (commonly with the same annotation type), which are stored in a
*/
- __pyx_t_2 = PyDict_New(); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 49; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __Pyx_GOTREF(((PyObject *)__pyx_t_2));
+ __pyx_t_1 = PyDict_New(); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 49; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_GOTREF(((PyObject *)__pyx_t_1));
/* "MACS2/IO/cFixWidthTrack.pyx":57
* dict. They can be sorted by calling self.sort() function.
@@ -7360,10 +7333,10 @@ PyMODINIT_FUNC PyInit_cFixWidthTrack(void)
* """fw is the fixed-width for all locations.
*
*/
- __pyx_t_1 = __pyx_binding_PyCFunctionType_NewEx(&__pyx_mdef_5MACS2_2IO_14cFixWidthTrack_10FWTrackIII___init__, NULL, __pyx_n_s_31); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 57; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __Pyx_GOTREF(__pyx_t_1);
- if (PyObject_SetItem(__pyx_t_2, __pyx_n_s____init__, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 57; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+ __pyx_t_2 = __pyx_binding_PyCFunctionType_NewEx(&__pyx_mdef_5MACS2_2IO_14cFixWidthTrack_10FWTrackIII___init__, NULL, __pyx_n_s_31); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 57; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_GOTREF(__pyx_t_2);
+ if (PyObject_SetItem(__pyx_t_1, __pyx_n_s____init__, __pyx_t_2) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 57; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
/* "MACS2/IO/cFixWidthTrack.pyx":69
*
@@ -7372,10 +7345,10 @@ PyMODINIT_FUNC PyInit_cFixWidthTrack(void)
* """Add a location to the list according to the sequence name.
*
*/
- __pyx_t_1 = __pyx_binding_PyCFunctionType_NewEx(&__pyx_mdef_5MACS2_2IO_14cFixWidthTrack_10FWTrackIII_1add_loc, NULL, __pyx_n_s_31); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 69; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __Pyx_GOTREF(__pyx_t_1);
- if (PyObject_SetItem(__pyx_t_2, __pyx_n_s__add_loc, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 69; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+ __pyx_t_2 = __pyx_binding_PyCFunctionType_NewEx(&__pyx_mdef_5MACS2_2IO_14cFixWidthTrack_10FWTrackIII_1add_loc, NULL, __pyx_n_s_31); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 69; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_GOTREF(__pyx_t_2);
+ if (PyObject_SetItem(__pyx_t_1, __pyx_n_s__add_loc, __pyx_t_2) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 69; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
/* "MACS2/IO/cFixWidthTrack.pyx":87
* self.__pointer[chromosome][strand] += 1
@@ -7384,10 +7357,10 @@ PyMODINIT_FUNC PyInit_cFixWidthTrack(void)
* arr.resize( arr.size + BUFFER_SIZE, refcheck = False )
* return
*/
- __pyx_t_1 = __pyx_binding_PyCFunctionType_NewEx(&__pyx_mdef_5MACS2_2IO_14cFixWidthTrack_10FWTrackIII_2__expand__, NULL, __pyx_n_s_31); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 87; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __Pyx_GOTREF(__pyx_t_1);
- if (PyObject_SetItem(__pyx_t_2, __pyx_n_s____expand__, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 87; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+ __pyx_t_2 = __pyx_binding_PyCFunctionType_NewEx(&__pyx_mdef_5MACS2_2IO_14cFixWidthTrack_10FWTrackIII_2__expand__, NULL, __pyx_n_s_31); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 87; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_GOTREF(__pyx_t_2);
+ if (PyObject_SetItem(__pyx_t_1, __pyx_n_s____expand__, __pyx_t_2) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 87; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
/* "MACS2/IO/cFixWidthTrack.pyx":91
* return
@@ -7396,10 +7369,10 @@ PyMODINIT_FUNC PyInit_cFixWidthTrack(void)
* """ Resize np arrays for 5' positions and sort them in place """
*
*/
- __pyx_t_1 = __pyx_binding_PyCFunctionType_NewEx(&__pyx_mdef_5MACS2_2IO_14cFixWidthTrack_10FWTrackIII_3finalize, NULL, __pyx_n_s_31); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 91; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __Pyx_GOTREF(__pyx_t_1);
- if (PyObject_SetItem(__pyx_t_2, __pyx_n_s__finalize, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 91; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+ __pyx_t_2 = __pyx_binding_PyCFunctionType_NewEx(&__pyx_mdef_5MACS2_2IO_14cFixWidthTrack_10FWTrackIII_3finalize, NULL, __pyx_n_s_31); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 91; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_GOTREF(__pyx_t_2);
+ if (PyObject_SetItem(__pyx_t_1, __pyx_n_s__finalize, __pyx_t_2) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 91; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
/* "MACS2/IO/cFixWidthTrack.pyx":112
* return
@@ -7408,10 +7381,10 @@ PyMODINIT_FUNC PyInit_cFixWidthTrack(void)
* """Return a tuple of two lists of locations for certain chromosome.
*
*/
- __pyx_t_1 = __pyx_binding_PyCFunctionType_NewEx(&__pyx_mdef_5MACS2_2IO_14cFixWidthTrack_10FWTrackIII_4get_locations_by_chr, NULL, __pyx_n_s_31); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 112; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __Pyx_GOTREF(__pyx_t_1);
- if (PyObject_SetItem(__pyx_t_2, __pyx_n_s_32, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 112; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+ __pyx_t_2 = __pyx_binding_PyCFunctionType_NewEx(&__pyx_mdef_5MACS2_2IO_14cFixWidthTrack_10FWTrackIII_4get_locations_by_chr, NULL, __pyx_n_s_31); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 112; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_GOTREF(__pyx_t_2);
+ if (PyObject_SetItem(__pyx_t_1, __pyx_n_s_32, __pyx_t_2) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 112; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
/* "MACS2/IO/cFixWidthTrack.pyx":121
* raise Exception("No such chromosome name (%s) in TrackI object!\n" % (chromosome))
@@ -7420,10 +7393,10 @@ PyMODINIT_FUNC PyInit_cFixWidthTrack(void)
* """Return all the chromosome names stored in this track object.
* """
*/
- __pyx_t_1 = __pyx_binding_PyCFunctionType_NewEx(&__pyx_mdef_5MACS2_2IO_14cFixWidthTrack_10FWTrackIII_5get_chr_names, NULL, __pyx_n_s_31); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 121; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __Pyx_GOTREF(__pyx_t_1);
- if (PyObject_SetItem(__pyx_t_2, __pyx_n_s__get_chr_names, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 121; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+ __pyx_t_2 = __pyx_binding_PyCFunctionType_NewEx(&__pyx_mdef_5MACS2_2IO_14cFixWidthTrack_10FWTrackIII_5get_chr_names, NULL, __pyx_n_s_31); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 121; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_GOTREF(__pyx_t_2);
+ if (PyObject_SetItem(__pyx_t_1, __pyx_n_s__get_chr_names, __pyx_t_2) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 121; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
/* "MACS2/IO/cFixWidthTrack.pyx":128
* return l
@@ -7432,10 +7405,10 @@ PyMODINIT_FUNC PyInit_cFixWidthTrack(void)
* """Total sequenced length = total number of tags * width of tag
* """
*/
- __pyx_t_1 = __pyx_binding_PyCFunctionType_NewEx(&__pyx_mdef_5MACS2_2IO_14cFixWidthTrack_10FWTrackIII_6length, NULL, __pyx_n_s_31); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 128; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __Pyx_GOTREF(__pyx_t_1);
- if (PyObject_SetItem(__pyx_t_2, __pyx_n_s__length, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 128; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+ __pyx_t_2 = __pyx_binding_PyCFunctionType_NewEx(&__pyx_mdef_5MACS2_2IO_14cFixWidthTrack_10FWTrackIII_6length, NULL, __pyx_n_s_31); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 128; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_GOTREF(__pyx_t_2);
+ if (PyObject_SetItem(__pyx_t_1, __pyx_n_s__length, __pyx_t_2) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 128; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
/* "MACS2/IO/cFixWidthTrack.pyx":133
* return self.total*self.fw
@@ -7444,10 +7417,10 @@ PyMODINIT_FUNC PyInit_cFixWidthTrack(void)
* """Naive sorting for locations.
*
*/
- __pyx_t_1 = __pyx_binding_PyCFunctionType_NewEx(&__pyx_mdef_5MACS2_2IO_14cFixWidthTrack_10FWTrackIII_7sort, NULL, __pyx_n_s_31); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 133; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __Pyx_GOTREF(__pyx_t_1);
- if (PyObject_SetItem(__pyx_t_2, __pyx_n_s__sort, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 133; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+ __pyx_t_2 = __pyx_binding_PyCFunctionType_NewEx(&__pyx_mdef_5MACS2_2IO_14cFixWidthTrack_10FWTrackIII_7sort, NULL, __pyx_n_s_31); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 133; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_GOTREF(__pyx_t_2);
+ if (PyObject_SetItem(__pyx_t_1, __pyx_n_s__sort, __pyx_t_2) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 133; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
/* "MACS2/IO/cFixWidthTrack.pyx":149
* self.__sorted = True
@@ -7456,16 +7429,16 @@ PyMODINIT_FUNC PyInit_cFixWidthTrack(void)
* """Filter the duplicated reads.
*
*/
- __pyx_t_1 = __Pyx_PyBool_FromLong(0); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 149; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __Pyx_GOTREF(__pyx_t_1);
- if (!(likely(__Pyx_TypeTest(__pyx_t_1, __pyx_ptype_7cpython_4bool_bool)))) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 149; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __pyx_k_3 = ((PyBoolObject *)__pyx_t_1);
- __Pyx_GIVEREF(__pyx_t_1);
- __pyx_t_1 = 0;
- __pyx_t_1 = __pyx_binding_PyCFunctionType_NewEx(&__pyx_mdef_5MACS2_2IO_14cFixWidthTrack_10FWTrackIII_8filter_dup, NULL, __pyx_n_s_31); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 149; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __Pyx_GOTREF(__pyx_t_1);
- if (PyObject_SetItem(__pyx_t_2, __pyx_n_s__filter_dup, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 149; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+ __pyx_t_2 = __Pyx_PyBool_FromLong(0); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 149; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_GOTREF(__pyx_t_2);
+ if (!(likely(__Pyx_TypeTest(__pyx_t_2, __pyx_ptype_7cpython_4bool_bool)))) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 149; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __pyx_k_3 = ((PyBoolObject *)__pyx_t_2);
+ __Pyx_GIVEREF(__pyx_t_2);
+ __pyx_t_2 = 0;
+ __pyx_t_2 = __pyx_binding_PyCFunctionType_NewEx(&__pyx_mdef_5MACS2_2IO_14cFixWidthTrack_10FWTrackIII_8filter_dup, NULL, __pyx_n_s_31); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 149; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_GOTREF(__pyx_t_2);
+ if (PyObject_SetItem(__pyx_t_1, __pyx_n_s__filter_dup, __pyx_t_2) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 149; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
/* "MACS2/IO/cFixWidthTrack.pyx":254
* return selfcopy
@@ -7474,10 +7447,10 @@ PyMODINIT_FUNC PyInit_cFixWidthTrack(void)
* """Sample the tags for a given percentage.
*
*/
- __pyx_t_1 = __pyx_binding_PyCFunctionType_NewEx(&__pyx_mdef_5MACS2_2IO_14cFixWidthTrack_10FWTrackIII_9sample_percent, NULL, __pyx_n_s_31); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 254; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __Pyx_GOTREF(__pyx_t_1);
- if (PyObject_SetItem(__pyx_t_2, __pyx_n_s__sample_percent, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 254; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+ __pyx_t_2 = __pyx_binding_PyCFunctionType_NewEx(&__pyx_mdef_5MACS2_2IO_14cFixWidthTrack_10FWTrackIII_9sample_percent, NULL, __pyx_n_s_31); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 254; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_GOTREF(__pyx_t_2);
+ if (PyObject_SetItem(__pyx_t_1, __pyx_n_s__sample_percent, __pyx_t_2) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 254; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
/* "MACS2/IO/cFixWidthTrack.pyx":287
* return
@@ -7486,10 +7459,10 @@ PyMODINIT_FUNC PyInit_cFixWidthTrack(void)
* """Sample the tags for a given percentage.
*
*/
- __pyx_t_1 = __pyx_binding_PyCFunctionType_NewEx(&__pyx_mdef_5MACS2_2IO_14cFixWidthTrack_10FWTrackIII_10sample_num, NULL, __pyx_n_s_31); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 287; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __Pyx_GOTREF(__pyx_t_1);
- if (PyObject_SetItem(__pyx_t_2, __pyx_n_s__sample_num, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 287; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+ __pyx_t_2 = __pyx_binding_PyCFunctionType_NewEx(&__pyx_mdef_5MACS2_2IO_14cFixWidthTrack_10FWTrackIII_10sample_num, NULL, __pyx_n_s_31); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 287; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_GOTREF(__pyx_t_2);
+ if (PyObject_SetItem(__pyx_t_1, __pyx_n_s__sample_num, __pyx_t_2) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 287; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
/* "MACS2/IO/cFixWidthTrack.pyx":298
* return
@@ -7498,10 +7471,10 @@ PyMODINIT_FUNC PyInit_cFixWidthTrack(void)
* """Output FWTrackIII to BED format files. If fhd is given,
* write to a file, otherwise, output to standard output.
*/
- __pyx_t_1 = __pyx_binding_PyCFunctionType_NewEx(&__pyx_mdef_5MACS2_2IO_14cFixWidthTrack_10FWTrackIII_11print_to_bed, NULL, __pyx_n_s_31); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 298; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __Pyx_GOTREF(__pyx_t_1);
- if (PyObject_SetItem(__pyx_t_2, __pyx_n_s__print_to_bed, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 298; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
+ __pyx_t_2 = __pyx_binding_PyCFunctionType_NewEx(&__pyx_mdef_5MACS2_2IO_14cFixWidthTrack_10FWTrackIII_11print_to_bed, NULL, __pyx_n_s_31); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 298; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_GOTREF(__pyx_t_2);
+ if (PyObject_SetItem(__pyx_t_1, __pyx_n_s__print_to_bed, __pyx_t_2) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 298; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
/* "MACS2/IO/cFixWidthTrack.pyx":49
* # ------------------------------------
@@ -7510,22 +7483,22 @@ PyMODINIT_FUNC PyInit_cFixWidthTrack(void)
* """Fixed Width Locations Track class II along the whole genome
* (commonly with the same annotation type), which are stored in a
*/
- if (PyDict_SetItemString(((PyObject *)__pyx_t_2), "__doc__", ((PyObject *)__pyx_kp_s_33)) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 49; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __pyx_t_1 = __Pyx_CreateClass(((PyObject *)__pyx_empty_tuple), ((PyObject *)__pyx_t_2), __pyx_n_s__FWTrackIII, __pyx_n_s_31); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 49; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __Pyx_GOTREF(__pyx_t_1);
- if (PyObject_SetAttr(__pyx_m, __pyx_n_s__FWTrackIII, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 49; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
- __Pyx_DECREF(((PyObject *)__pyx_t_2)); __pyx_t_2 = 0;
+ if (PyDict_SetItemString(((PyObject *)__pyx_t_1), "__doc__", ((PyObject *)__pyx_kp_s_33)) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 49; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __pyx_t_2 = __Pyx_CreateClass(((PyObject *)__pyx_empty_tuple), ((PyObject *)__pyx_t_1), __pyx_n_s__FWTrackIII, __pyx_n_s_31); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 49; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_GOTREF(__pyx_t_2);
+ if (PyObject_SetAttr(__pyx_m, __pyx_n_s__FWTrackIII, __pyx_t_2) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 49; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
+ __Pyx_DECREF(((PyObject *)__pyx_t_1)); __pyx_t_1 = 0;
/* "MACS2/IO/cFixWidthTrack.pyx":1
* # cython: profile=True # <<<<<<<<<<<<<<
- * # Time-stamp: <2012-05-01 21:52:01 Tao Liu>
+ * # Time-stamp: <2012-05-02 01:07:31 Tao Liu>
*
*/
- __pyx_t_2 = PyDict_New(); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __Pyx_GOTREF(((PyObject *)__pyx_t_2));
- if (PyObject_SetAttr(__pyx_m, __pyx_n_s____test__, ((PyObject *)__pyx_t_2)) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- __Pyx_DECREF(((PyObject *)__pyx_t_2)); __pyx_t_2 = 0;
+ __pyx_t_1 = PyDict_New(); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_GOTREF(((PyObject *)__pyx_t_1));
+ if (PyObject_SetAttr(__pyx_m, __pyx_n_s____test__, ((PyObject *)__pyx_t_1)) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __Pyx_DECREF(((PyObject *)__pyx_t_1)); __pyx_t_1 = 0;
/* "cpython/type.pxd":2
*
View
4 MACS2/IO/cFixWidthTrack.pyx
@@ -1,5 +1,5 @@
# cython: profile=True
-# Time-stamp: <2012-05-01 21:52:01 Tao Liu>
+# Time-stamp: <2012-05-02 01:07:31 Tao Liu>
"""Module for FWTrack classes.
@@ -21,7 +21,7 @@ with the distribution).
import logging
#from array import array
-from random import sample as random_sample
+#from random import sample as random_sample
import sys
from copy import copy
from MACS2.Constants import *
View
5,296 MACS2/IO/cScoreTrack.c
2,731 additions, 2,565 deletions not shown because the diff is too large. Please use a local Git client to view these changes.
View
73 MACS2/IO/cScoreTrack.pyx
@@ -1,5 +1,5 @@
# cython: profile=True
-# Time-stamp: <2012-05-01 18:15:11 Tao Liu>
+# Time-stamp: <2012-05-02 16:56:21 Tao Liu>
"""Module for Feature IO classes.
@@ -129,14 +129,21 @@ class scoreTrackI:
this class is 0-indexed and right-open.
"""
- def __init__ (self):
+ def __init__ (self, float effective_depth_in_million = 1.0):
"""Different with bedGraphTrackI, missing values are simply
replaced with 0.
-
+
+ effective_depth_in_million: sequencing depth in million after
+ duplicates being filtered. If
+ treatment is scaled down to
+ control sample size, then this
+ should be control sample size in
+ million. And vice versa.
"""
self.data = {}
self.pointer = {}
self.trackline = False
+ self.effective_depth_in_million = effective_depth_in_million
def enable_trackline(self):
"""Turn on trackline with bedgraph output
@@ -153,7 +160,7 @@ class scoreTrackI:
'100logLR': np.zeros(chrom_max_len, dtype="int32") }
self.pointer[chrom] = 0
- def add (self, str chromosome, int endpos, int sample, float control):
+ def add (self, str chromosome, int endpos, float sample, float control):
"""Add a chr-endpos-sample-control block into data
dictionary. At the mean time, calculate pvalues and log
likelihood.
@@ -166,7 +173,7 @@ class scoreTrackI:
c['pos'][i] = endpos
c['sample'][i] = sample
c['control'][i] = control
- c['-100logp'][i] = get_pscore(sample,control)
+ c['-100logp'][i] = get_pscore(int(sample),control)
c['100logLR'][i] = logLR(sample,control)
self.pointer[chromosome] += 1
@@ -198,51 +205,57 @@ class scoreTrackI:
l = set(self.data.keys())
return l
- def write_bedGraph (self, fhd, str name, str description, str colname):
+ def write_bedGraph ( self, fhd, str name, str description, str colname, bool do_SPMR = False ):
"""Write all data to fhd in Wiggle Format.
fhd: a filehandler to save bedGraph.
+
name/description: the name and description in track line.
- colname: can be 'sample','control','-100logp','-100logq'
+ colname: can be 'sample','control','-100logp','-100logq', '100logLR'
+ do_SPMR: only effective when writing sample/control tracks. When True, save SPMR instead.
+
"""
cdef str chrom
cdef int l, pre, i, p
- cdef double pre_v, v
+ cdef float pre_v, v, scale_factor
if self.trackline:
# this line is REQUIRED by the wiggle format for UCSC browser
- fhd.write("track type=bedGraph name=\"%s\" description=\"%s\"\n" % (name,description))
+ fhd.write( "track type=bedGraph name=\"%s\" description=\"%s\"\n" % ( name,description ) )
+
+ if colname not in [ 'sample', 'control', '-100logp', '-100logq', '100logLR' ]:
+ raise Exception( "%s not supported!" % colname )
+
+ if colname in [ '-100logp', '-100logq', '100logLR' ]:
+ scale_factor = 0.1 # for pvalue or qvalue, divide them by 100 while writing to bedGraph file
+ elif colname in [ 'sample', 'control' ]:
+ if do_SPMR:
+ logging.info( "MACS will save SPMR for fragment pileup using effective depth of %.2f million" % self.effective_depth_in_million )
+ scale_factor = 1.0/self.effective_depth_in_million
+ else:
+ scale_factor = 1
- if colname not in ['sample','control','-100logp','-100logq','100logLR']:
- raise Exception("%s not supported!" % colname)
- if colname in ['-100logp', '-100logq','100logLR']:
- flag100 = True # for pvalue or qvalue, divide them by 100 while writing to bedGraph file
- else:
- flag100 = False
chrs = self.get_chr_names()
for chrom in chrs:
- d = self.data[chrom]
- l = self.pointer[chrom]
+ d = self.data[ chrom ]
+ l = self.pointer[ chrom ]
pre = 0
- pos = d['pos']
- if flag100:
- value = d[colname]/100.0
- else:
- value = d[colname]
- if len(value) == 0: continue
- pre_v = value[0]
+ pos = d[ 'pos' ]
+ value = d[ colname ] * scale_factor
+ if value.shape[ 0 ] == 0: continue # skip if there's no data
+ pre_v = value[ 0 ]
for i in range( 1, l ):
- v = value[i]
- p = pos[i-1]
+ v = value[ i ]
+ p = pos[ i-1 ]
if pre_v != v:
- fhd.write("%s\t%d\t%d\t%.2f\n" % (chrom,pre,p,pre_v))
+ fhd.write( "%s\t%d\t%d\t%.2f\n" % ( chrom, pre, p, pre_v ) )
pre_v = v
pre = p
- p = pos[-1]
+ p = pos[ -1 ]
# last one
- fhd.write("%s\t%d\t%d\t%.2f\n" % (chrom,pre,p,pre_v))
-
+ fhd.write( "%s\t%d\t%d\t%.2f\n" % ( chrom, pre, p, pre_v ) )
+
return True
def make_pq_table ( self ):
View
5 MACS2/OptValidator.py
@@ -1,4 +1,4 @@
-# Time-stamp: <2012-04-10 17:56:07 Tao Liu>
+# Time-stamp: <2012-05-02 16:38:42 Tao Liu>
"""Module Description
@@ -178,6 +178,9 @@ def opt_validate ( options ):
if options.halfext:
options.argtxt += "# MACS will make 1/2d size fragments\n"
+ if options.do_SPMR and options.store_bdg:
+ options.argtxt += "# MACS will save fragment pileup signal per million reads\n"
+
return options
View
6,456 MACS2/cPeakDetect.c
3,271 additions, 3,185 deletions not shown because the diff is too large. Please use a local Git client to view these changes.
View
69 MACS2/cPeakDetect.pyx
@@ -1,5 +1,5 @@
# cython: profile=True
-# Time-stamp: <2012-05-01 21:56:08 Tao Liu>
+# Time-stamp: <2012-05-02 17:05:58 Tao Liu>
"""Module Description
@@ -15,12 +15,13 @@ with the distribution).
@author: Yong Zhang, Tao Liu
@contact: taoliu@jimmy.harvard.edu
"""
-import os
-from array import array
-from copy import deepcopy
+#import os
+#from array import array
+#from copy import deepcopy
from itertools import groupby
from operator import itemgetter
-import subprocess
+import io
+#import subprocess
import gc # use garbage collectior
from MACS2.IO.cPeakIO import PeakIO
@@ -59,6 +60,8 @@ def compare_treatment_vs_control(treat, control, fragment_size, gsize,
Finally, BH process will be applied to adjust pvalue to qvalue.
"""
+ cdef float effective_depth_in_million
+
if PE_MODE:
treat_total = treat[0].total
control_total = control[0].total
@@ -78,8 +81,10 @@ def compare_treatment_vs_control(treat, control, fragment_size, gsize,
# if user want to scale everything to control data
lambda_bg = float(fragment_size)*treat_total/gsize/ratio_treat2control
treat_btrack.apply_func(lambda x:float(x)/ratio_treat2control)
+ effective_depth_in_million = control_total / 1000000.0
else:
lambda_bg = float(fragment_size)*treat_total/gsize
+ effective_depth_in_million = treat_total / 1000000.0
# control data needs multiple steps of calculation
@@ -138,7 +143,7 @@ def compare_treatment_vs_control(treat, control, fragment_size, gsize,
scale_factor_s=scale_factor_s)
# calculate pvalue scores
- score_btrack = treat_btrack.make_scoreTrack_for_macs(control_btrack)
+ score_btrack = treat_btrack.make_scoreTrack_for_macs( control_btrack , effective_depth_in_million = effective_depth_in_million )
#treat_btrack = None # clean them
#control_btrack = None
#gc.collect() # full collect garbage
@@ -316,7 +321,8 @@ class PeakDetect:
for enrichment.
"""
cdef int i
-
+ cdef float lambda_bg, effective_depth_in_million
+
if self.PE_MODE:
treat_total = self.treat[0].total
control_total = self.control[0].total
@@ -327,19 +333,26 @@ class PeakDetect:
self.info("#3 pileup treatment data by extending tags towards 3' to %d length" % self.d)
self.ratio_treat2control = float(treat_total)/control_total
-
if self.opt.tocontrol:
- # if user want to scale everything to control data
+ # if MACS decides to scale treatment to control data because treatment is bigger
+
+ effective_depth_in_million = control_total / 1000000.0
+
lambda_bg = float(self.d)*treat_total/self.gsize/self.ratio_treat2control
+
if self.PE_MODE:
treat_btrack = pileup_frag_bdg(self.treat[0], self.treat[1],
- scale_factor=1/self.ratio_treat2control)
+ scale_factor=1/self.ratio_treat2control)
else:
treat_btrack = pileup_bdg(self.treat, self.d, directional=True,
halfextension=self.opt.halfext,
- scale_factor=1/self.ratio_treat2control)
+ scale_factor=1/self.ratio_treat2control)
else:
+ # if MACS decides to scale control to treatment because control sample is bigger
+ effective_depth_in_million = treat_total / 1000000.0
+
lambda_bg = float(self.d)*treat_total/self.gsize
+
if self.PE_MODE:
treat_btrack = pileup_frag_bdg(self.treat[0], self.treat[1],
scale_factor=1.0)
@@ -393,17 +406,17 @@ class PeakDetect:
# pileup using different extension sizes and scaling factors
if self.PE_MODE:
control_btrack = pileup_frag_w_multiple_d_bdg(self.control[0],
- self.control[1],
- d_s[1:],
- baseline_value=lambda_bg,
- scale_factor_s=scale_factor_s[1:],
- scale_factor_0=scale_factor_s[0])
+ self.control[1],
+ d_s[1:],
+ baseline_value=lambda_bg,
+ scale_factor_s=scale_factor_s[1:],
+ scale_factor_0=scale_factor_s[0])
else:
control_btrack = pileup_w_multiple_d_bdg(self.control, d_s,
- baseline_value=lambda_bg,
- directional=self.shiftcontrol,
- halfextension=self.opt.halfext,
- scale_factor_s=scale_factor_s)
+ baseline_value=lambda_bg,
+ directional=self.shiftcontrol,
+ halfextension=self.opt.halfext,
+ scale_factor_s=scale_factor_s)
# free mem
#self.treat = None
@@ -414,7 +427,7 @@ class PeakDetect:
# calculate pvalue scores
self.info("#3 Build score track ...")
- score_btrack = treat_btrack.make_scoreTrack_for_macs(control_btrack)
+ score_btrack = treat_btrack.make_scoreTrack_for_macs( control_btrack, effective_depth_in_million = effective_depth_in_million )
if self.opt.trackline: score_btrack.enable_trackline()
#treat_btrack = None # clean them
#control_btrack = None
@@ -424,7 +437,7 @@ class PeakDetect:
pqtable = score_btrack.make_pq_table()
self.info("#3 Saving p-value to q-value table ...")
- pqfhd = open(self.opt.pqtable,"w")
+ pqfhd = io.open(self.opt.pqtable,"wb")
pqfhd.write( "-log10pvalue\t-log10qvalue\trank\tbasepairs\n" )
for p in sorted(pqtable.keys(),reverse=True):
#for i in range(pqtable.shape[0]):
@@ -484,9 +497,9 @@ class PeakDetect:
for filename, title, desc, scorecol in tracks:
self.info("#3 save the %s track into bedGraph file..." % desc)
- with open(filename, 'w') as bdgfhd:
+ with io.open(filename, 'wb') as bdgfhd:
score_btrack.write_bedGraph(bdgfhd, title,
- trackdesc % desc, scorecol)
+ trackdesc % desc, scorecol, do_SPMR = self.opt.do_SPMR) # do_SPMR doesn't have effect on p/q/logLR values.
return peaks
def __call_peaks_wo_control (self):
@@ -512,6 +525,10 @@ class PeakDetect:
Finally, a poisson CDF is applied to calculate one-side pvalue
for enrichment.
"""
+ cdef float effective_depth_in_million
+
+ effective_depth_in_million = treat_total / 1000000.0
+
# global lambda
if self.PE_MODE:
# should we support halfext?
@@ -555,7 +572,7 @@ class PeakDetect:
# calculate pvalue scores
self.info("#3 Build score track ...")
- score_btrack = treat_btrack.make_scoreTrack_for_macs(control_btrack)
+ score_btrack = treat_btrack.make_scoreTrack_for_macs( control_btrack, effective_depth_in_million = effective_depth_in_million )
if self.opt.trackline: score_btrack.enable_trackline()
treat_btrack = None # clean them
control_btrack = None
@@ -619,7 +636,7 @@ class PeakDetect:
self.info("#3 save the %s track into bedGraph file..." % desc)
with open(filename, 'w') as bdgfhd:
score_btrack.write_bedGraph(bdgfhd, title,
- trackdesc % desc, scorecol)
+ trackdesc % desc, scorecol, do_SPMR = self.opt.do_SPMR) # do_SPMR doesn't have effect on p/q/logLR values.
return peaks
View
17 bin/macs2
@@ -1,5 +1,5 @@
#!/usr/bin/env python
-# Time-stamp: <2012-04-10 17:56:16 Tao Liu>
+# Time-stamp: <2012-05-02 16:38:16 Tao Liu>
"""Description: MACS 2 main executable
@@ -146,7 +146,11 @@ def add_callpeak_parser( subparsers ):
default = False )
group_output.add_argument( "--verbose", dest = "verbose", type = int, default = 2,
help = "Set verbose level of runtime message. 0: only show critical message, 1: show additional warning message, 2: show process information, 3: show debug messages. DEFAULT:2" )
+ group_output.add_argument( "--no-trackline", dest="trackline", action="store_false", default = True,
+ help = "Tells MACS not to include trackline with bedGraph files. The trackline is required by UCSC. Require -B to be set. Default: keep trackline." )
+ group_output.add_argument( "--SPMR", dest = "do_SPMR", action = "store_true", default = False,
+ help = "If True, MACS will save signal per million reads for fragment pileup profiles. Require -B to be set. Default: False" )
# group for bimodal
group_bimodal = argparser_callpeak.add_argument_group( "Shifting model arguments" )
group_bimodal.add_argument( "-s", "--tsize", dest = "tsize", type = int, default = None,
@@ -172,6 +176,11 @@ def add_callpeak_parser( subparsers ):
help = "Minimum FDR (q-value) cutoff for peak detection. DEFAULT: 0.05. -q and -p are mutually exclusive." )
p_or_q_group.add_argument( "-p", "--pvalue", dest = "pvalue", type = float,
help = "Pvalue cutoff for peak detection. DEFAULT: not set. -q and -p are mutually exclusive." )
+ # about scaling
+ group_callpeak.add_argument( "--to-large", dest = "tolarge", action = "store_true", default = False,
+ help = "When set, scale the small sample up to the bigger sample. By default, the bigger dataset will be scaled down towards the smaller dataset, which will lead to smaller p/qvalues and more specific results. Keep in mind that scaling down will bring down background noise more. DEFAULT: False" )
+ group_callpeak.add_argument( "--down-sample", dest = "downsample", action = "store_true", default = False,
+ help = "When set, random sampling method will scale down the bigger sample. By default, MACS uses linear scaling. Warning: This option will make your result unstable and irreproducible since each time, random reads would be selected. Consider to use 'randsample' script instead. <not implmented>If used together with --SPMR, 1 million unique reads will be randomly picked.</not implemented> Caution: due to the implementation, the final number of selected reads may not be as you expected! DEFAULT: False" )
group_callpeak.add_argument( "--nolambda", dest = "nolambda", action = "store_true",
help = "If True, MACS will use fixed background lambda as local lambda for every peak region. Normally, MACS calculates a dynamic local lambda to reflect the local bias due to potential chromatin structure. ",
@@ -180,10 +189,6 @@ def add_callpeak_parser( subparsers ):
help = "The small nearby region in basepairs to calculate dynamic lambda. This is used to capture the bias near the peak summit region. Invalid if there is no control data. If you set this to 0, MACS will skip slocal lambda calculation. *Note* that MACS will always perform a d-size local lambda calculation. The final local bias should be the maximum of the lambda value from d, slocal, and llocal size windows. DEFAULT: 1000 " )
group_callpeak.add_argument( "--llocal", dest = "largelocal", type = int, default = 10000,
help = "The large nearby region in basepairs to calculate dynamic lambda. This is used to capture the surround bias. If you set this to 0, MACS will skip llocal lambda calculation. *Note* that MACS will always perform a d-size local lambda calculation. The final local bias should be the maximum of the lambda value from d, slocal, and llocal size windows. DEFAULT: 10000." )
- group_callpeak.add_argument( "--to-large", dest = "tolarge", action = "store_true", default = False,
- help = "When set, scale the small sample up to the bigger sample. By default, the bigger dataset will be scaled down towards the smaller dataset, which will lead to smaller p/qvalues and more specific results. Keep in mind that scaling down will bring down background noise more. DEFAULT: False" )
- group_callpeak.add_argument( "--down-sample", dest = "downsample", action = "store_true", default = False,
- help = "When set, random sampling method will scale down the bigger sample. By default, MACS uses linear scaling. Warning: This option will make your result unstable and irreproducible since each time, random reads would be selected. Consider to use 'randsample' script instead. DEFAULT: False" )
group_callpeak.add_argument( "--shift-control", dest = "shiftcontrol", action = "store_true", default = False,
help = "When set, control tags will be shifted just as ChIP tags according to their strand before the extension of d, slocal and llocal. By default, control tags are extended centered at their current positions regardless of strand. You may consider to turn this option on while comparing two ChIP datasets of different condition but the same factor. DEFAULT: False" )
group_callpeak.add_argument( "--half-ext", dest = "halfext", action = "store_true", default = False,
@@ -194,8 +199,6 @@ def add_callpeak_parser( subparsers ):
help = "Cutoff for broad region. This option is not available unless --broad is set. If -p is set, this is a pvalue cutoff, otherwise, it's a qvalue cutoff. DEFAULT: 0.1 " )
group_callpeak.add_argument( "--call-summits", dest="call_summits", action="store_true",
help="If set, MACS will use a more sophisticated approach to find all summits in each enriched peak region. DEFAULT: False",default=False)
- group_callpeak.add_argument("--no-trackline", dest="trackline", action="store_false", default=True,
- help="Tells MACS not to include trackline with bedGraph files. The trackline is required by UCSC.")
return

0 comments on commit 21951a2

Please sign in to comment.