Python pyBigWig.open函数代码示例

def create_new_header(infile, mappings, outfile):
    """Create new header in BigWig, with UCSC chromosome names."""
    with pyBigWig.open(infile) as bw:
        if set(bw.chroms().keys()).issubset(mappings.values()):
            # If chromosome names are already UCSC, just rename input file to output name.
            # Exit with status 0 since this is normal behavior.
            os.rename(infile, outfile)
            sys.exit(0)

        hdr = [(mappings[chrom], length) for chrom, length in bw.chroms().items() if chrom in mappings]

        if not hdr:
            msg = "Neither of the chromosomes in the input file has a valid UCSC pair. No mapping will be done."
            print(warning(msg))
            os.rename(infile, outfile)
            sys.exit(0)

        seq_num = 0
        with pyBigWig.open(outfile, 'w') as bw_output:
            bw_output.addHeader(hdr)
            for chrom, length in bw.chroms().items():
                ints = bw.intervals(chrom, 0, length)
                if ints and chrom in mappings:
                    bw_output.addEntries([mappings[chrom]] * len(ints),
                                         [x[0] for x in ints],
                                         ends=[x[1] for x in ints],
                                         values=[x[2] for x in ints])
                elif chrom not in mappings:
                    seq_num += 1
                    print('UCSC chromosome/conting mapping for {} is missing'.format(chrom))

        if seq_num > 0:
            print(warning("UCSC chromosome/conting mapping for {} sequence(s) is missing. "
                          "This sequence(s) will not be included in the bigWig file.".format(seq_num)))

def getChromSizes(bigwigFilesList):
    """
    Get chromosome sizes from bigWig file with pyBigWig

    Test dataset with two samples covering 200 bp.
    >>> test = Tester()

    Chromosome name(s) and size(s).
    >>> assert(getChromSizes([test.bwFile1, test.bwFile2]) == ([('3R', 200)], set([])))
    """
    # check that the path to USCS bedGraphToBigWig as set in the config
    # is installed and is executable.

    def print_chr_names_and_size(chr_set):
        sys.stderr.write("chromosome\tlength\n")
        for name, size in chr_set:
            sys.stderr.write("{0:>15}\t{1:>10}\n".format(name, size))

    bigwigFilesList = bigwigFilesList[:]

    common_chr = set()
    for fname in bigwigFilesList:
        fh = pyBigWig.open(fname)
        common_chr = common_chr.union(set(fh.chroms().items()))
        fh.close()

    non_common_chr = set()
    for bw in bigwigFilesList:
        _names_and_size = set(pyBigWig.open(bw).chroms().items())
        if len(common_chr & _names_and_size) == 0:
            #  try to add remove 'chr' from the chromosme name
            _corr_names_size = set()
            for chrom_name, size in _names_and_size:
                if chrom_name.startswith('chr'):
                    _corr_names_size.add((chrom_name[3:], size))
                else:
                    _corr_names_size.add(('chr' + chrom_name, size))
            if len(common_chr & _corr_names_size) == 0:
                message = "No common chromosomes found. Are the bigwig files " \
                          "from the same species and same assemblies?\n"
                sys.stderr.write(message)
                print_chr_names_and_size(common_chr)

                sys.stderr.write("\nand the following is the list of the unmatched chromosome and chromosome\n"
                                 "lengths from file\n{}\n".format(bw))
                print_chr_names_and_size(_names_and_size)
                exit(1)
            else:
                _names_and_size = _corr_names_size

        non_common_chr |= common_chr ^ _names_and_size
        common_chr = common_chr & _names_and_size

    if len(non_common_chr) > 0:
        sys.stderr.write("\nThe following chromosome names did not match between the the bigwig files\n")
        print_chr_names_and_size(non_common_chr)

    # get the list of common chromosome names and sizes
    return sorted(common_chr), non_common_chr

def main():
    usage = 'usage: %prog [options] <in_bw_file> <out_h5_file>'
    parser = OptionParser(usage)
    parser.add_option('-v', dest='verbose', default=False, action='store_true')
    (options,args) = parser.parse_args()

    if len(args) != 2:
        parser.error('Must provide input BigWig and output HDF5.')
    else:
        bw_file = args[0]
        hdf5_file = args[1]

    # open files
    bw_in = pyBigWig.open(bw_file)
    h5_out = h5py.File(hdf5_file, 'w')

    # for each chromosome
    chrom_lengths = bw_in.chroms()
    for chrom in chrom_lengths:
        if options.verbose:
            print(chrom)

        # read values
        x = bw_in.values(chrom, 0, chrom_lengths[chrom], numpy=True).astype('float16')

        # write gzipped into HDF5
        h5_out.create_dataset(chrom, data=x, dtype='float16', compression='gzip', shuffle=True)

    # close files
    h5_out.close()
    bw_in.close()

def mhsmidkernelsmooth(bamfile, bwfile, maxinsert=80, mininsert=1, paired=False, kernelsize=30):
    bamfor = Baminfo.Baminfo(bamfile)

    bw = pyBigWig.open(bwfile, "w")

    bw.addHeader(list(bamfor.chrlen.items()))

    for chromosome in bamfor.chrlen:

        end = bamfor.chrlen[chromosome]

        mhsmidcount = mhsbam.mhsmidcount(bamfile=bamfile, chromosome=chromosome, start=1,
                                         end=end, maxinsert=maxinsert, mininsert=mininsert, paired=paired)

        mhsmidsmoothed = kernelsmooth(mhsmidcount, 1, end, end, kernelsize)

        if mhsmidsmoothed:

            starts = list()

            values = list()

            for start in sorted(mhsmidsmoothed):
                starts.append(start)

                values.append(float(mhsmidsmoothed[start]))

            bw.addEntries(chromosome, starts=starts, values=values,
                          span=1, step=1)

    bw.close()

def dhscutkernelsmooth(bamfile, bwfile, library='Duke', kernelsize=200):
    bamfor = Baminfo.Baminfo(bamfile)

    bw = pyBigWig.open(bwfile, "w")

    bw.addHeader(list(bamfor.chrlen.items()))

    for chromosome in bamfor.chrlen:

        end = bamfor.chrlen[chromosome]

        dhscut = dhsbam.dhcutcount(bamfile=bamfile, chromosome=chromosome, start=1,
                                   end=end, library=library)

        dhscutsmoothed = kernelsmooth(dhscut, 1, end, end, kernelsize)

        if dhscutsmoothed:

            starts = list()

            values = list()

            for start in sorted(dhscutsmoothed):
                starts.append(start)

                values.append(float(dhscutsmoothed[start]))

            bw.addEntries(chromosome, starts=starts, values=values,
                          span=1, step=1)

    bw.close()

    def testBigBed(self):
        fname = "http://www.encodeproject.org/files/ENCFF001JBR/@@download/ENCFF001JBR.bigBed"
        bb = pyBigWig.open(fname)
        assert(bb is not None)
        assert(bb.isBigWig() == 0)
        assert(bb.isBigBed() == 1)
        SQL = """table RnaElements 
"BED6 + 3 scores for RNA Elements data "
    (
    string chrom;      "Reference sequence chromosome or scaffold"
    uint   chromStart; "Start position in chromosome"
    uint   chromEnd;   "End position in chromosome"
    string name;       "Name of item"
    uint   score;      "Normalized score from 0-1000"
    char[1] strand;    "+ or - or . for unknown"
    float level;       "Expression level such as RPKM or FPKM. Set to -1 for no data."
    float signif;      "Statistical significance such as IDR. Set to -1 for no data."
    uint score2;       "Additional measurement/count e.g. number of reads. Set to 0 for no data."
    )
"""
        output = bb.SQL()
        if isinstance(output, bytes):
            output = output.decode('ASCII')
        assert(output == SQL)
        o = bb.entries('chr1',10000000,10020000)
        expected = [(10009333, 10009640, '61035\t130\t-\t0.026\t0.42\t404'), (10014007, 10014289, '61047\t136\t-\t0.029\t0.42\t404'), (10014373, 10024307, '61048\t630\t-\t5.420\t0.00\t2672399')]
        assert(o == expected)
        bb.close()

 def _generate_chunk_output_file(self, i=None):
     records = [
         ("chr1", 1, 2, 1.5),
         ("chr1", 2, 3, 4.5),
         ("chr1", 3, 4, 1.9),
         ("chr1", 4, 5, 0.45),
         ("chr2", 8, 9, 1.0),
         ("chr2", 9, 10, 6.7)
     ]
     fn = tempfile.NamedTemporaryFile(suffix=".bw").name
     _records = records[(i*3):(i*3)+3]
     assert len(_records) == 3
     ranges = {}
     for rec in _records:
         seqid = rec[0]
         pos = rec[1]
         ranges.setdefault(seqid, (sys.maxint, 0))
         ranges[seqid] = (min(ranges[seqid][0], pos),
                          max(ranges[seqid][1], pos))
     bw = pyBigWig.open(fn, "w")
     regions = [ (s, ranges[s][1]+1) for s in sorted(ranges.keys()) ]
     bw.addHeader(regions)
     bw.addEntries([rec[0] for rec in _records],
                   [rec[1]-1 for rec in _records],
                   ends=[rec[2]-1 for rec in _records],
                   values=[rec[3] for rec in _records])
     bw.close()
     return fn

    def coverage_from_bigwig(self, bigwig_file, stepsize=100):

        """Return list of arrays describing the coverage of each genomicRegions from <bigwig_file>.
        
        *Keyword arguments:*
        
        - bigwig_file -- path to bigwig file
        - stepsize -- used stepsize
        
        *Output:*
        
        Class variable <coverage>: a list where the elements correspond to the GenomicRegion. The list elements give
        the number of reads falling into the GenomicRegion.
        
        """

        self.coverage = []

        bwf = pyBigWig.open(bigwig_file)

        for gr in self.genomicRegions:
            steps = int(len(gr) / stepsize)
            try:
                ds = bwf.stats(gr.chrom, gr.initial, gr.final, type="mean", nBins=steps)
                ds = [x if x else 0 for x in ds]
            except:
                ds = [0] * steps
            self.coverage.append(np.array(ds))

        bwf.close()

def coveragetobw(bamfile, bwfile, maxinsert, mininsert, paired=False):
    bamfor = Baminfo.Baminfo(bamfile)

    bw = pyBigWig.open(bwfile, "w")

    bw.addHeader(list(bamfor.chrlen.items()))

    for chromosome in bamfor.chrlen:

        end = bamfor.chrlen[chromosome]

        coveragecount = mhsbam.coveragecount(bamfile=bamfile, chromosome=chromosome, start=1,
                                             end=end, maxinsert=maxinsert, mininsert=mininsert, paired=paired)

        if coveragecount:

            starts = list()

            values = list()

            for start in sorted(coveragecount):
                starts.append(start)

                values.append(float(coveragecount[start]))

            bw.addEntries(chromosome, starts=starts, values=values,
                          span=1, step=1)

    bw.close()

def gerprunner():

    import pyBigWig

    b = pyBigWig.open("/scratch/ucgd/lustre/u1021864/serial/hg19.gerp.bw")
   # x = list(range(1,23)); x.append("X"), x.append("Y")

    input = sys.argv[1]
    iterator = JimFile(input)
    iterable = windower(iterator, chunker(1))
    cutoff = 1e-3

    def genchunks():
        nsmall = 0
        for i, chunk in enumerate(iterable):
            #if len(chunk) < 5:
            #    continue
            score = b.stats("chr"+chunk[0].chrom, chunk[0].start, chunk[-1].end)
            yield chunk, score[0]
            if i % 100000 == 0:
                print i, chunk[0].chrom, chunk[0].start, score
        print >>sys.stderr, nsmall, "removed for being too short"
        print >>sys.stderr, i, "total chunks"

    vcf_path = "/scratch/ucgd/lustre/u1021864/serial/clinvar-anno.vcf.gz"
    res = eval2(genchunks(), vcf_path,
        "/scratch/ucgd/lustre/u1021864/serial/esp-common.vcf.gz")
    print metrics(res[True], res[False], "gerp.auc.png")

    def doWrite2(self):
        '''
        Test all three modes of storing entries. Also test to ensure that we get error messages when doing something silly

        This is a modified version of the writing example from libBigWig
        '''
        chroms = ["1"]*6
        starts = [0, 100, 125, 200, 220, 230, 500, 600, 625, 700, 800, 850]
        ends = [5, 120, 126, 205, 226, 231]
        values = [0.0, 1.0, 200.0, -2.0, 150.0, 25.0, 0.0, 1.0, 200.0, -2.0, 150.0, 25.0, -5.0, -20.0, 25.0, -5.0, -20.0, 25.0]
        ofile = tempfile.NamedTemporaryFile(delete=False)
        oname = ofile.name
        ofile.close()
        bw = pyBigWig.open(oname, "w")
        bw.addHeader([("1", 1000000), ("2", 1500000)])

        #Intervals
        bw.addEntries(chroms[0:3], starts[0:3], ends=ends[0:3], values=values[0:3])
        bw.addEntries(chroms[3:6], starts[3:6], ends=ends[3:6], values=values[3:6])

        #IntervalSpans
        bw.addEntries("1", starts[6:9], values=values[6:9], span=20)
        bw.addEntries("1", starts[9:12], values=values[9:12], span=20)

        #IntervalSpanSteps, this should instead take an int
        bw.addEntries("1", 900, values=values[12:15], span=20, step=30)
        bw.addEntries("1", 990, values=values[15:18], span=20, step=30)

        #Attempt to add incorrect values. These MUST raise an exception
        try:
            bw.addEntries(chroms[0:3], starts[0:3], ends=ends[0:3], values=values[0:3])
            assert(1==0)
        except RuntimeError:
            pass
        try:
            bw.addEntries("1", starts[6:9], values=values[6:9], span=20)
            assert(1==0)
        except RuntimeError:
            pass
        try:
            bw.addEntries("3", starts[6:9], values=values[6:9], span=20)
            assert(1==0)
        except RuntimeError:
            pass
        try:
            bw.addEntries("1", 900, values=values[12:15], span=20, step=30)
            assert(1==0)
        except RuntimeError:
            pass

        #Add a few intervals on a new chromosome
        bw.addEntries(["2"]*3, starts[0:3], ends=ends[0:3], values=values[0:3])
        bw.close()
        #check md5sum, this is the simplest method to check correctness
        h = hashlib.md5(open(oname, "rb").read()).hexdigest()
        assert(h=="b1ca91d2ff42afdd2efa19a007c1ded4")
        #Clean up
        os.remove(oname)

def fetch_from_bigbed(path, chrom, start, end):
    import pyBigWig

    bed = pyBigWig.open(path)
    assert bed.isBigBed(), "Oops, for some reason I was expecting a bed file: {}".format(path)

    chrom = match_chrom_format(chrom, bed.chroms().keys())
    for cur_start, cur_end, bed_line in bed.entries(chrom, start, end):
        bed_line = bed_line.split()
        yield tx_from_bedfields([chrom, cur_start, cur_end] + bed_line)

    def __init__(self, wig_location):
        """
        Arguments
        ---------
        wig_location: Path to bigwig

        """
        self.wig_location = wig_location
        try:
            self.wig = pyBigWig.open(self.wig_location)
        except Exception as e:
            raise MocaException('Error reading wig file: {}'.format(e))

def bedGraphToBigWig(chromSizes, bedGraphPath, bigWigPath, sort=True):
    """
    takes a bedgraph file, orders it and converts it to
    a bigwig file using pyBigWig.
    """

    from tempfile import NamedTemporaryFile
    from os import remove, system

    # Make a list of tuples for the bigWig header, this MUST be sorted identically to the bedGraph file
    sort_cmd = cfg.config.get('external_tools', 'sort')
    _file = NamedTemporaryFile(delete=False)
    for chrom, size in chromSizes:
        _file.write(toBytes("{}\t{}\n".format(chrom, size)))
    _file.close()
    system("LC_ALL=C {} -k1,1 -k2,2n {} > {}.sorted".format(sort_cmd, _file.name, _file.name))
    cl = []
    f = open("{}.sorted".format(_file.name))
    for line in f:
        chrom, chromLen = line.split()
        cl.append((chrom, int(chromLen)))
    f.close()
    remove(_file.name)
    remove("{}.sorted".format(_file.name))

    # check if the file is empty
    if os.stat(bedGraphPath).st_size < 10:
        import sys
        sys.stderr.write(
            "Error: The generated bedGraphFile was empty. Please adjust\n"
            "your deepTools settings and check your input files.\n")
        exit(1)

    if sort:
        # temporary file to store sorted bedgraph file
        _file = NamedTemporaryFile(delete=False)
        tempfilename1 = _file.name
        system("LC_ALL=C {} -k1,1 -k2,2n {} > {}".format(sort_cmd, bedGraphPath, tempfilename1))
        bedGraphPath = tempfilename1

    bw = pyBigWig.open(bigWigPath, "w")
    assert(bw is not None)
    # The lack of maxZooms will change the results a bit, perhaps the defaults are better
    bw.addHeader(cl, maxZooms=10)
    f = open(bedGraphPath)
    for line in f:
        interval = line.split()
        bw.addEntries([interval[0]], [int(interval[1])], ends=[int(interval[2])], values=[float(interval[3])])
    f.close()
    bw.close()

    if sort:
        remove(tempfilename1)

    def readValuesPyBigWig(self, reference, start, end):
        """
        Use pyBigWig package to read a BigWig file for the
        given range and return a protocol object.

        pyBigWig returns an array of values that fill the query range.
        Not sure if it is possible to get the step and span.

        This method trims NaN values from the start and end.

        pyBigWig throws an exception if end is outside of the
        reference range. This function checks the query range
        and throws its own exceptions to avoid the ones thrown
        by pyBigWig.
        """
        if not self.checkReference(reference):
            raise exceptions.ReferenceNameNotFoundException(reference)
        if start < 0:
            start = 0
        bw = pyBigWig.open(self._sourceFile)
        referenceLen = bw.chroms(reference)
        if referenceLen is None:
            raise exceptions.ReferenceNameNotFoundException(reference)
        if end > referenceLen:
            end = referenceLen
        if start >= end:
            raise exceptions.ReferenceRangeErrorException(
                reference, start, end)

        data = protocol.Continuous()
        curStart = start
        curEnd = curStart + self._INCREMENT
        while curStart < end:
            if curEnd > end:
                curEnd = end
            for i, val in enumerate(bw.values(reference, curStart, curEnd)):
                if not math.isnan(val):
                    if len(data.values) == 0:
                        data.start = curStart + i
                    data.values.append(val)
                    if len(data.values) == self._MAX_VALUES:
                        yield data
                        data = protocol.Continuous()
                elif len(data.values) > 0:
                    # data.values.append(float('NaN'))
                    yield data
                    data = protocol.Continuous()
            curStart = curEnd
            curEnd = curStart + self._INCREMENT

        bw.close()
        if len(data.values) > 0:
            yield data