Python BigWigFile.get方法代码示例

# 需要导入模块: from bx.bbi.bigwig_file import BigWigFile [as 别名]
# 或者: from bx.bbi.bigwig_file.BigWigFile import get [as 别名]
def get_phastcons(bedtool, phastcons_location, species=None, index=None, ):
    
    """
    
    Get phastcons scores for intervals in a bed tool
    
    """
    
    if species is None and index is None:
        print "Error, must select species or index"
    
    f = open(phastcons_location, 'r')
    bw = BigWigFile(file=f)

    try:
        
        #if its a line
        #for each line fetch bigwig values 
        type(bedtool)
        v = bedtool.chrom #is a single interval
        vals = bw.get(bedtool.chrom, bedtool.start, bedtool.stop)
        consvals = list(v[-1] for v in vals)
        if len(consvals) > 0:
            mean_phastcons = np.mean(consvals)
        else:
            mean_phastcons=0
        data = mean_phastcons
    except:
        
        #if bedtool
        for i, bedline in enumerate(bedtool):
            data = np.ndarray(len(bedtool))        
            vals = bw.get(bedline.chrom, bedline.start, bedline.stop)
            consvals = list(v[-1] for v in vals)
            if len(consvals) > 0:
                mean_phastcons = np.mean(consvals)
            else:
                mean_phastcons=0
            data[i] = mean_phastcons
            
    #returns mean phastcons score for each line 
    #returns inconistant data types, need to convert so it just returns an array 
    return data

# 需要导入模块: from bx.bbi.bigwig_file import BigWigFile [as 别名]
# 或者: from bx.bbi.bigwig_file.BigWigFile import get [as 别名]
 def get_GA_from_bw(self, plus, minus, GTF, filterfxn):
     ##bx-python 'get' method is 0 based, fully closed
     ga = HTSeq.GenomicArray( "auto", typecode='d' , stranded = True)
     with open(plus) as f:
         bw_file = BigWigFile(file=f)
         for GF in GTF:
             if filterfxn( GF ) == False: continue
             window = GF.iv
             chrom, start, stop = window.chrom, window.start, window.end
             vals = bw_file.get(chrom, start, stop)
             for start, stop, value in vals:
                 ga[ HTSeq.GenomicPosition(chrom, start, '+') ] = value
     with open(minus) as f:
         bw_file = BigWigFile(file=f)
         for GF in GTF:
             if filterfxn( GF ) == False: continue
             window = GF.iv
             chrom, start, stop = window.chrom, window.start, window.end
             vals = bw_file.get(chrom, start, stop)
             for start, stop, value in vals:
                 ga[ HTSeq.GenomicPosition(chrom, start, '-') ] = value
     return ga

# 需要导入模块: from bx.bbi.bigwig_file import BigWigFile [as 别名]
# 或者: from bx.bbi.bigwig_file.BigWigFile import get [as 别名]
def get_phastcons(bedtool, species=None, index=None):
    """
    Get phastcons scores for intervals in a bed tool
    """
    if species is None and index is None:
        print "Error, must select species or index"
    if species is not None and index is None:
        if species == "mm9":
            index= basedir + "/yeolab/Conservation/phastCons/mm9_30way/placental/mm9_phastcons.bw"
        elif species == "hg19":
            index = basedir + "/yeolab/Conservation/phastCons/hg19_46way/placentalMammals/reformat/hg19_phastcons.bw"
    f = open(index, 'r')
    bw = BigWigFile(file=f)

    try:
        type(bedtool)
        v = bedtool.chrom #is a single interval
        vals = bw.get(bedtool.chrom, bedtool.start, bedtool.stop)
        consvals = list(v[-1] for v in vals)
        if len(consvals) > 0:
            mean_phastcons = np.mean(consvals)
        else:
            mean_phastcons=0
        data = mean_phastcons


    except:
        for i, bedline in enumerate(bedtool):
            data = np.ndarray(len(bedtool))        
            vals = bw.get(bedline.chrom, bedline.start, bedline.stop)
            consvals = list(v[-1] for v in vals)
            if len(consvals) > 0:
                mean_phastcons = np.mean(consvals)
            else:
                mean_phastcons=0
            data[i] = mean_phastcons
    return data

# 需要导入模块: from bx.bbi.bigwig_file import BigWigFile [as 别名]
# 或者: from bx.bbi.bigwig_file.BigWigFile import get [as 别名]
def main():
    input_filename, output_filename, loc_filename, loc_key, chrom_col, start_col = sys.argv[1:]

    # open input, output, and bigwig files
    location_file = LocationFile( loc_filename )
    bigwig_filename = location_file.get_values( loc_key )
    bwfh = open_or_die( bigwig_filename, message='Error opening BigWig file %s' % bigwig_filename )
    bw = BigWigFile( file=bwfh )
    ifh = open_or_die( input_filename, message='Error opening input file %s' % input_filename )
    ofh = open_or_die( output_filename, mode='w', message='Error opening output file %s' % output_filename )

    # make column numbers 0-based
    chrom_col = int( chrom_col ) - 1
    start_col = int( start_col ) - 1
    min_cols = max( chrom_col, start_col )

    # add score column to imput file
    line_number = 0
    for line in ifh:
        line_number += 1
        line = line.rstrip( '\r\n' )
        elems = line.split( '\t' )
        if len( elems ) > min_cols:
            chrom = elems[chrom_col].strip()
            # base-0 position in chrom
            start = int( elems[start_col] )
            score_list = bw.get( chrom, start, start + 1 )
            score_list_len = len( score_list )
            if score_list_len == 1:
                beg, end, score = score_list[0]
                score_val = '%1.3f' % score
            elif score_list_len == 0:
                score_val = 'NA'
            else:
                die( '%s line %d: chrom=%s, start=%d, score_list_len = %d' % ( input_filename, line_number, chrom, start, score_list_len ) )
            print('\t'.join( [line, score_val] ), file=ofh)
        else:
            print(line, file=ofh)

    bwfh.close()
    ifh.close()
    ofh.close()

# 需要导入模块: from bx.bbi.bigwig_file import BigWigFile [as 别名]
# 或者: from bx.bbi.bigwig_file.BigWigFile import get [as 别名]
def getNumberOfFragmentsPerRegionFromBigWig(bw, chromSizes):
    """
    Get the number of all mapped fragments per region in all chromosomes
    from a bigWig. Utilizing bx-python.

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

    Get number of fragments in sample.
    >>> getNumberOfFragmentsPerRegionFromBigWig(test.bwFile1, [('3R', 200)])
    3.0
    >>> getNumberOfFragmentsPerRegionFromBigWig(test.bwFile2, [('3R', 200)])
    4.0
    """
    bwh = BigWigFile(open(bw, "rb"))
    mapped = 0
    for cname, csize in chromSizes:
        regions = bwh.get(cname, 0, csize) # region = bwh.get(chrom_name, start, end)
        for region in regions:
            mapped += region[2]
    return mapped

# 需要导入模块: from bx.bbi.bigwig_file import BigWigFile [as 别名]
# 或者: from bx.bbi.bigwig_file.BigWigFile import get [as 别名]
import numpy as np
fl=sys.argv[1]
dist=int(sys.argv[2])
from bx.bbi.bigwig_file import BigWigFile

genes=read.dat("/home/ssaberi/resources/list.genes.txt",'\t')
table=read.dat("/projects/epigenomics/MarcoJuliaPon/peaks.txt",'\t')
mygenes=read.dat("/projects/epigenomics/MarcoJuliaPon/mygenes.txt",'\t')
ens=[]
for i in mygenes:
	for gn in genes:
		if i in gn[0]:
			ens.append(gn[1])
			break

genespos=read.read_gene_pos('/home/ssaberi/resources/hg19v69_genes.TSS_2000.pc.A03480.H3K27me3.GE02.coverage')
genesbed=bedtools.makebed_genpos(ens,genespos,100000)
              


f = open(fl)
bw = BigWigFile(file=f)
mat=[]
for bed_i in genesbed:
   vals = bw.get( bed_i[0], bed_i[1], bed_i[2])
   mat.append(np.array(vals))
mat=np.array(mat)
plt.matshow(mat,aspect='auto',cmap='YlOrBr')
fl=fl[-fl[::-1].index('/'):-fl[::-1].index('.')]
plt.save(fl+".pdf")

# 需要导入模块: from bx.bbi.bigwig_file import BigWigFile [as 别名]
# 或者: from bx.bbi.bigwig_file.BigWigFile import get [as 别名]
class BigWig(object):
    def __init__(self, filename):
        self.filename = filename
        self.determine_sizes()
        self.bwf = BigWigFile(open(filename))

    def determine_sizes(self):
        self.sizes = {}
        fh = open(self.filename, "rb")
        # read magic number to guess endianness
        magic = fh.read(4)
        if magic == '&\xfc\x8f\x88':
            endianness = '<'
        elif magic == '\x88\x8f\xfc&':
            endianness = '>'
        else:
            raise IOError("The file is not in bigwig format")

        # read the header
        info = struct.unpack(endianness + 'HHQQQHHQQIQ', fh.read(60))
        self.version = info[0]
        self.zoom_levels = info[1]
        self.chromosome_tree_offset = info[2]
        self.full_data_offset = info[3]
        self.full_index_offset = info[4]
        self.field_count = info[5]
        self.defined_field_count = info[6]
        self.auto_SQL_offset = info[7]
        self.total_summary_offset = info[8]
        self.uncompress_buf_size = info[9]
        
        # go to the data
        fh.seek(self.chromosome_tree_offset)
        # read magic again
        magic = fh.read(4)
        if magic == '\x91\x8c\xcax':
            endianness = '<'
        elif magic == 'x\xca\x8c\x91':
            endianness = '>'
        else:
            raise ValueError("Wrong magic for this bigwig data file")

        info2 = struct.unpack(endianness + 'IIIQQ', fh.read(28))
        self.block_size = info2[0]
        self.key_size = info2[1]
        self.val_size = info2[2]
        self.item_count = info2[3]

        info3 = struct.unpack(endianness + 'BBH', fh.read(4))
        self.is_leaf = info3[0]
        self.count = info3[2]

        for n in range(self.count):
            format_code = endianness + str(self.key_size) + 'sII'
            info = struct.unpack(format_code, fh.read(self.key_size + 2 * 4))
            key, chrom_id, chrom_size = info

            key = key.replace('\x00', '')
            self.sizes[key] = chrom_size

    def get_as_array(self, chrom, start, end):
        return self.bwf.get_as_array(chrom, start, end)

    def get(self, chrom, start, end):
        return self.bwf.get(chrom, start, end)

    def query(self, chrom, start, end, number):
        return self.bwf.query(chrom, start, end, number)