Python BigWigFile.get_as_array方法代码示例

# 需要导入模块: from bx.bbi.bigwig_file import BigWigFile [as 别名]
# 或者: from bx.bbi.bigwig_file.BigWigFile import get_as_array [as 别名]
def main():
	usage="%prog [options]"
	parser = OptionParser(usage,version="%prog " + __version__)
	
	parser.add_option("-i","--bwfile1",action="store",type="string",dest="BigWig_File1",help="BigWig files")
	parser.add_option("-j","--bwfile2",action="store",type="string",dest="BigWig_File2",help="BigWig files")
	parser.add_option("-a","--action",action="store",type="string",dest="action",help='After pairwise align two bigwig files, perform the follow actions (Only select one keyword):"Add" = add signals. "Average" = average signals. "Division"= divide bigwig2 from bigwig1. Add 1 to both bigwig. "Max" = pick the signal that is larger. "Min" = pick the signal that is smaller. "Product" = multiply signals. "Subtract" = subtract signals in 2nd bigwig file from the corresponiding ones in the 1st bigwig file. "geometricMean" = take the geometric mean of signals.')
	parser.add_option("-o","--output",action="store",type="string",dest="output_wig",help="Output wig file")
	parser.add_option("-s","--chromSize",action="store",type="string",dest="chromSize",help="Chromosome size file. Tab or space separated text file with 2 columns: first column is chromosome name, second column is size of the chromosome.")
	parser.add_option("-c","--chunk",action="store",type="int",dest="chunk_size",default=100000,help="Chromosome chunk size. Each chomosome will be cut into samll chunks of this size. Decrease chunk size will save more RAM. default=%default (bp)")
	(options,args)=parser.parse_args()
	
	if not (options.BigWig_File1 and options.BigWig_File2  and options.output_wig and options.chromSize):
		parser.print_help()
		sys.exit(0)
	OUT=open(options.output_wig,'w')
	bw1 = BigWigFile( file=open(options.BigWig_File1) )
	bw2 = BigWigFile( file=open(options.BigWig_File2) )
	chrom_sizes = load_chromsize(options.chromSize)
	for chr_name, chr_size in chrom_sizes.items():		#iterate each chrom
		print >>sys.stderr, "Processing " + chr_name + " ..."
		OUT.write('variableStep chrom='+chr_name+'\n')
		for interval in BED.tillingBed(chrName = chr_name,chrSize = chr_size,stepSize = options.chunk_size):
			coord = interval[1]
			bw_signal1 = bw1.get_as_array(chr_name,interval[1],interval[2])
			bw_signal2 = bw2.get_as_array(chr_name,interval[1],interval[2])
			if all_nan(bw_signal1) and all_nan(bw_signal2):
				continue
			bw_signal1 = replace_nan( bw_signal1 )
			bw_signal2 = replace_nan( bw_signal2 )
		
			call_back = getattr(twoList,options.action)
			for v in call_back(bw_signal1,bw_signal2):
				coord +=1
				if v != 0: print >>OUT, "%d\t%.2f" % (coord,v)

# 需要导入模块: from bx.bbi.bigwig_file import BigWigFile [as 别名]
# 或者: from bx.bbi.bigwig_file.BigWigFile import get_as_array [as 别名]
def get_mean_phastcons(bedtool, phastcons_location, sample_size = 1000):
    
    """
    
    Get means phastcons scores for all intervals in a bed tool
    bedtool - bedtool to extract data from
    phastcons_location - location of phastcons file
    
    """
    
    with open(phastcons_location) as bw_file:
        bw = BigWigFile(bw_file)
    
        data = []
        
        for bedline in bedtool.random_subset(min(len(bedtool), sample_size)):
            conservation_values = bw.get_as_array(bedline.chrom, bedline.start, bedline.stop)
            try:
                if len(conservation_values) > 0:
                    mean_phastcons = np.mean(conservation_values)
                else:
                    mean_phastcons = 0
                data.append(mean_phastcons)
            except TypeError:
                pass
    return data

# 需要导入模块: from bx.bbi.bigwig_file import BigWigFile [as 别名]
# 或者: from bx.bbi.bigwig_file.BigWigFile import get_as_array [as 别名]
def get_mean_phastcons(bedtool, phastcons_location):
    
    """
    
    Get means phastcons scores for all intervals in a bed tool
    bedtool - bedtool to extract data from
    phastcons_location - location of phastcons file
    
    """
    
    f = open(phastcons_location, 'r')
    bw = BigWigFile(file=f)

    #if bedtool
    data = np.ndarray(len(bedtool))  
    for i, bedline in enumerate(bedtool):
              
        conservation_values = bw.get_as_array(bedline.chrom, bedline.start, bedline.stop)
        
        if len(conservation_values) > 0:
            mean_phastcons = np.mean(conservation_values)
        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_array [as 别名]
def findInsertions(bwFile, bedData, interval, x):

    if interval =='start':
        sL = int(bedData[x][1])-options.l
        sR = int(bedData[x][1])+options.r
    elif interval == 'end':
        sL = int(bedData[x][2])-options.l
        sR = int(bedData[x][2])+options.r
    else:
        sL = int(bedData[x][1])-options.l
        sR = int(bedData[x][2])+options.r
    
    # get signal data
    f = open(bwFile, "rb")
    bigwig_class = BigWigFile(f)
    try: signal = bigwig_class.get_as_array(bedData[x][0],sL,sR)
    except OverflowError: signal = np.array([np.nan]*(sR-sL))
    f.close()
    
    if signal is not None:
        if np.sum(np.isfinite(signal)) > 0:
            out = np.nanmean(signal)
        else: out = 0
    else: out = 0
    
    out = signal
    return out

# 需要导入模块: from bx.bbi.bigwig_file import BigWigFile [as 别名]
# 或者: from bx.bbi.bigwig_file.BigWigFile import get_as_array [as 别名]
    def summarize(self, interval, bins=None, method='summarize',
                  function='mean'):

        # We may be dividing by zero in some cases, which raises a warning in
        # NumPy based on the IEEE 754 standard (see
        # http://docs.scipy.org/doc/numpy/reference/generated/
        #       numpy.seterr.html)
        #
        # That's OK -- we're expecting that to happen sometimes. So temporarily
        # disable this error reporting for the duration of this method.
        orig = np.geterr()['invalid']
        np.seterr(invalid='ignore')

        if (bins is None) or (method == 'get_as_array'):
            bw = BigWigFile(open(self.fn))
            s = bw.get_as_array(
                interval.chrom,
                interval.start,
                interval.stop,)
            if s is None:
                s = np.zeros((interval.stop - interval.start,))
            else:
                s[np.isnan(s)] = 0

        elif method == 'ucsc_summarize':
            if function in ['mean', 'min', 'max', 'std', 'coverage']:
                return self.ucsc_summarize(interval, bins, function=function)
            else:
                raise ValueError('function "%s" not supported by UCSC\'s'
                                 'bigWigSummary')

        else:
            bw = BigWigFile(open(self.fn))
            s = bw.summarize(
                interval.chrom,
                interval.start,
                interval.stop, bins)
            if s is None:
                s = np.zeros((bins,))
            else:
                if function == 'sum':
                    s = s.sum_data
                if function == 'mean':
                    s = s.sum_data / s.valid_count
                    s[np.isnan(s)] = 0
                if function == 'min':
                    s = s.min_val
                    s[np.isinf(s)] = 0
                if function == 'max':
                    s = s.max_val
                    s[np.isinf(s)] = 0
                if function == 'std':
                    s = (s.sum_squares / s.valid_count)
                    s[np.isnan(s)] = 0

        # Reset NumPy error reporting
        np.seterr(divide=orig)
        return s

# 需要导入模块: from bx.bbi.bigwig_file import BigWigFile [as 别名]
# 或者: from bx.bbi.bigwig_file.BigWigFile import get_as_array [as 别名]
class BigWigWrapper(object):

    """A wrapper for bx-python BigWig file"""

    def __init__(self, filepath):
        self.bw = BigWigFile(open(filepath))

    def __getitem__(self, iv):
        return self.bw.get_as_array(iv.chrom, iv.start, iv.end)

# 需要导入模块: from bx.bbi.bigwig_file import BigWigFile [as 别名]
# 或者: from bx.bbi.bigwig_file.BigWigFile import get_as_array [as 别名]
def profile_bwfile(inbed,bwfile):
	'''retrieve signal from bigwig file for each entry in input bed file'''
	bw = BigWigFile( file=open( bwfile ) )
	
	for line in open(inbed):
		bw_signal=[]
		try:
			if line.startswith('#'):continue
			if line.startswith('track'):continue
			if line.startswith('browser'):continue
			if not line.strip():
				continue
			else:
				line = line.rstrip('\r\n')
				fields = line.split()
				chrom = fields[0]
				start = int(fields[1])
				end = int(fields[2])
		except:
			print >>sys.stderr,"Must be  chrom [space] start [space] end: " + line,
			continue		
		bw_signal.extend(bw.get_as_array(chrom,start,end))
		print chrom +'\t'+ str(start) +'\t'+ str(end) + '\t' + ','.join(str(i) for i in bw_signal)		

# 需要导入模块: from bx.bbi.bigwig_file import BigWigFile [as 别名]
# 或者: from bx.bbi.bigwig_file.BigWigFile import get_as_array [as 别名]
def findInsertions(bwFile, bedData, x):
    if options.tn5 is not None:
        bwFile = options.b + options.tn5 + "." + bedData[x][0] + ".Scores.bw"

    sL = int(bedData[x][1]) - options.l
    sR = int(bedData[x][2]) + options.r

    # get signal data
    f = open(bwFile, "rb")
    bw = BigWigFile(f)
    try:
        signal = bw.get_as_array(bedData[x][0], sL, sR)
    except OverflowError:
        signal = np.array([np.nan] * (sR - sL))
    f.close()

    out = signal
    try:
        if bedData[x][3] == "-":
            out = out[::-1]
    except IndexError:
        pass
    return out

# 需要导入模块: from bx.bbi.bigwig_file import BigWigFile [as 别名]
# 或者: from bx.bbi.bigwig_file.BigWigFile import get_as_array [as 别名]
 or left_or_right == 'r' and strand == '-'):
 # 5' site
 if (chrom, coordinate) in annotated_5p:
     fivep_splice_site_counts = (
             annotated_fivep_splice_site_counts
         )
     line_counts = annotated_line_counts
 else:
     fivep_splice_site_counts = (
             unannotated_fivep_splice_site_counts
         )
     line_counts = unannotated_line_counts
 if strand == '+':
     bwvals = bw.get_as_array(
                 chrom,
                 coordinate - args.extension,
                 coordinate + args.extension
             )
     if bwvals is None:
         continue
     for i, j in enumerate(
                     xrange(-args.extension, args.extension)
                 ):
         if not math.isnan(bwvals[i]):
             fivep_splice_site_counts[j] += bwvals[i]
             line_counts[j] += 1
 elif strand == '-':
     bwvals = bw.get_as_array(
                 chrom,
                 coordinate - (args.extension - 1),
                 coordinate + (args.extension + 1)

# 需要导入模块: from bx.bbi.bigwig_file import BigWigFile [as 别名]
# 或者: from bx.bbi.bigwig_file.BigWigFile import get_as_array [as 别名]
Create a site profile vector showing the average signal accumulated from a
bigwig file around the center of each interval from a BED file.

Output is the average signal value at that relative position across the 
intervals.

usage: %prog bigwig_file.bw padding < bed_file.bed 
"""

import sys
from numpy import *

from bx.intervals.io import GenomicIntervalReader
from bx.bbi.bigwig_file import BigWigFile

bw = BigWigFile( open( sys.argv[1] ) )
padding = int( sys.argv[2] )
totals = zeros( padding*2, dtype=float64 )
valid = zeros( padding*2, dtype=int32 )

for interval in GenomicIntervalReader( sys.stdin ):
    center = floor( ( interval.start + interval.end ) / 2 )
    values = bw.get_as_array( interval.chrom, center - padding, center + padding )
    # Determine which positions had data and mask the rest for totalling
    invalid = isnan( values )
    values[ invalid ] = 0
    totals += values
    valid += ( ~ invalid )

savetxt( sys.stdout, totals/valid )

# 需要导入模块: from bx.bbi.bigwig_file import BigWigFile [as 别名]
# 或者: from bx.bbi.bigwig_file.BigWigFile import get_as_array [as 别名]
	if args.inFile2 is not None:
		inFile2 = BedReader(open(args.inFile2))
else:	
	raise Exception("Unrecognized format!");



for locus in scanThese:
	if args.verbose>0:
		print("Scanning %s"%(locus[GENOMEDATA.NAME]))
	stF = max(locus[GENOMEDATA.ST] - padding,0);
	enF = locus[GENOMEDATA.EN] + padding + inclusive;
	try:
		if (locus[GENOMEDATA.STR]=="-" and args.inFile2 is not None):
			if args.format=="BIGWIG" or args.format=="BW" or args.format=="BIGBED" or args.format=="BB":
				values = inFile2.get_as_array( locus[GENOMEDATA.CHR], stF, enF )
		else:
			if args.format=="BIGWIG" or args.format=="BW" or args.format=="BIGBED" or args.format=="BB":
				values = inFile1.get_as_array( locus[GENOMEDATA.CHR], stF, enF )
	except OverflowError as e:
		sys.stderr.write("OverflowError at '%s'; st=%d, en=%d\n"%(locus[GENOMEDATA.NAME],locus[GENOMEDATA.ST],locus[GENOMEDATA.EN]));
		raise(e);
	if values is None and args.correctChr>0:
		#try again adding chr or taking it away
		if locus[GENOMEDATA.CHR][:3]=="chr":
			 locus[GENOMEDATA.CHR] =  locus[GENOMEDATA.CHR][3:]
		else:
			 locus[GENOMEDATA.CHR]="chr"+ locus[GENOMEDATA.CHR];
		try:
			if (locus[GENOMEDATA.STR]=="-" and args.inFile2 is not None):
				if args.format=="BIGWIG" or args.format=="BW" or args.format=="BIGBED" or args.format=="BB":

# 需要导入模块: from bx.bbi.bigwig_file import BigWigFile [as 别名]
# 或者: from bx.bbi.bigwig_file.BigWigFile import get_as_array [as 别名]
def main(args):
    bw_file = BigWigFile( open(args.bigWigFile) )
    bw_file.get_as_array(chrom, st, end)

# 需要导入模块: from bx.bbi.bigwig_file import BigWigFile [as 别名]
# 或者: from bx.bbi.bigwig_file.BigWigFile import get_as_array [as 别名]
inFile.close();

chromSizesFile = MYUTILS.smartGZOpen(args.chrsFile,'r');
chromSizes = {};
for line in chromSizesFile:
	if line is None or line == "" or line[0]=="#": continue
	data=line.rstrip().split("\t");
	chromSizes[data[0]]=int(data[1]);

curBW = BigWigFile(open(args.inBW))

outStream = MYUTILS.smartGZOpen("%s.wig.gz"%(args.outFPre),"w");
outStream.write("track type=wiggle_0\n")

for chr in transChrs:
	values = curBW.get_as_array( chr, 0, chromSizes[oldToNew[chr]] )
	#print(chr);
	if values is not None:
		sys.stderr.write("Adding %s -> %s\n"%(chr, oldToNew[chr]));
		outStream.write("fixedStep chrom=%s start=1 step=1\n"%(oldToNew[chr]))
		outStream.write("\n".join(map(str,values)));
		outStream.write("\n");

toBW = subprocess.Popen(["wigToBigWig","%s.wig.gz"%(args.outFPre),args.chrsFile,"%s.bw"%(args.outFPre)])
temp = toBW.communicate()
if temp[0] is not None:
	sys.stderr.write("wigToBigWig: %s"%(temp[0]));
if temp[1] is not None:
	sys.stderr.write("wigToBigWig: %s"%(temp[1]));
if temp[0] is None and temp[1] is None and os.path.isfile("%s.bw"%(args.outFPre)): # if no errors, delete the original
	os.remove("%s.wig.gz"%(args.outFPre))

# 需要导入模块: from bx.bbi.bigwig_file import BigWigFile [as 别名]
# 或者: from bx.bbi.bigwig_file.BigWigFile import get_as_array [as 别名]
wSize = options.w
wSmooth = np.ones(wSize)
step = options.s

#### SCRIPT #####
# split genome into chunks
for i in range(0,len(gSizes)):
    # break chrs into pieces
    chrN = gSizes[i][0]
    chrLen = int(gSizes[i][1])
    sVals = np.arange(1,int(gSizes[i][1]),chunkSize)
    
    # read in bigWig
    for j in range(0,len(sVals)):
        # get data, pass if not available
        signal = bw.get_as_array(chrN,sVals[j],sVals[j]+chunkSize+padLen)
        try: signal.any()
        except: continue
        
        # smooth data
        print chrN, sVals[j]
        signal[np.isnan(signal)] = 0
        convM = np.convolve(signal,wSmooth,'same')
            
        # save data
        sList = np.arange(sVals[j],sVals[j]+chunkSize+padLen,step)
        eList = sList+step
        chrList = np.array([chrN]*len(sList))
        meanSig = convM[range(step/2,chunkSize+padLen,step)]
        
        # save out

# 需要导入模块: from bx.bbi.bigwig_file import BigWigFile [as 别名]
# 或者: from bx.bbi.bigwig_file.BigWigFile import get_as_array [as 别名]
curBW = BigWigFile(open(args.inBW))

outStream = MYUTILS.smartGZOpen("%s.wig.gz"%(args.outFPre),"w");
outStream.write("track type=wiggle_0\n")


for chr in chromSizes.keys():
	last = 0;
	final = chromSizes[chr];
	sys.stderr.write("Outputting data for %s:\n"%(chr));
	while last!=final: # this breaks it up into chunks so that I'm not piping entire (human) chromosomes at once
		if args.verbose>0: sys.stderr.write("  Section %i - %i:\n"%(last,curLast));
		curLast = np.min([last+args.chunks,final]);
		curEnd = np.min([curLast+additionalFlankSize, final]);
		curSt = np.max([last-additionalFlankSize,0]);
		values = curBW.get_as_array( chr, curSt, curEnd )
	#print(chr);
		if values is not None:
			for f in allFunctions:
				values = applyFunction(values,f);
			values = values[(last - curSt):(curLast-last + (last-curSt))];# set them only to the middle part of this data so that the additionalFlankSize regions are not output.
			#print(values.shape);
			if last==0:
				outStream.write("fixedStep chrom=%s start=1 step=1\n"%(chr))
			outStream.write("\n".join(map(str,values)));
			outStream.write("\n");
			outStream.flush();
		last=curLast;

outStream.close();