Python pyBigWig.open方法代碼示例

# 需要導入模塊: import pyBigWig [as 別名]
# 或者: from pyBigWig import open [as 別名]
def read_blacklist(blacklist_bed, black_buffer=20):
  """Construct interval trees of blacklist
     regions for each chromosome."""
  black_chr_trees = {}

  if blacklist_bed is not None and os.path.isfile(blacklist_bed):
    for line in open(blacklist_bed):
      a = line.split()
      chrm = a[0]
      start = max(0, int(a[1]) - black_buffer)
      end = int(a[2]) + black_buffer

      if chrm not in black_chr_trees:
        black_chr_trees[chrm] = intervaltree.IntervalTree()

      black_chr_trees[chrm][start:end] = True

  return black_chr_trees 

# 需要導入模塊: import pyBigWig [as 別名]
# 或者: from pyBigWig import open [as 別名]
def are_files_equal(file1, file2):
    equal = True
    with open(file1) as textfile1, open(file2) as textfile2:
        for x, y in zip(textfile1, textfile2):
            if x.startswith('File'):
                continue
            if x != y:
                # handle the case of flipped values
                split_x = x.split('\t')
                split_y = y.split('\t')
                if split_x[0] == split_y[0] and split_x[1] == split_y[1] and split_x[2] == split_y[2]:
                    # to ignore rounding errors after 2th digit
                    if 0 <= abs(abs(float(split_x[3].strip())) - abs(float(split_y[3].strip()))) <= 0.01:
                        continue
                    else:
                        log.debug('split_x {} split_y {}'.format(split_x, split_y))
                equal = False
                break
    return equal 

# 需要導入模塊: import pyBigWig [as 別名]
# 或者: from pyBigWig import open [as 別名]
def extract_fasta_to_file(fasta, output_dir, mode="bcolz", overwrite=False):
    assert mode in _array_writer

    makedirs(output_dir, exist_ok=overwrite)
    fasta_file = FastaFile(fasta)
    file_shapes = {}
    for chrom, size in zip(fasta_file.references, fasta_file.lengths):
        data = np.zeros((size, NUM_SEQ_CHARS), dtype=np.float32)
        seq = fasta_file.fetch(chrom)
        one_hot_encode_sequence(seq, data)
        file_shapes[chrom] = data.shape
        _array_writer[mode](data, os.path.join(output_dir, chrom))

    with open(os.path.join(output_dir, "metadata.json"), "w") as fp:
        json.dump(
            {
                "file_shapes": file_shapes,
                "type": "array_{}".format(mode),
                "source": fasta,
            },
            fp,
        ) 

# 需要導入模塊: import pyBigWig [as 別名]
# 或者: from pyBigWig import open [as 別名]
def skiprows(f):

    try:
        import gzip
        fh = gzip.open(f)
        for i, l in enumerate(fh):
            if l.decode()[0] != "#":
                break
    except (OSError, TypeError):  # not a gzipped file, or StringIO
        fh = open(f)
        for i, l in enumerate(fh):
            if l[0] != "#":
                break

    fh.close()

    return i 

# 需要導入模塊: import pyBigWig [as 別名]
# 或者: from pyBigWig import open [as 別名]
def __get_scores_per_bin(self, genome_range, num_bins, max_try_nums=5):
        # on rare occasions pyBigWig may throw an error, apparently caused by a corruption
        # of the memory. This only occurs when calling trackPlot from different
        # processors. Reloading the file solves the problem.
        num_tries = 0
        while num_tries < max_try_nums:
            num_tries += 1
            try:
                scores_per_bin = np.array(self.bw.stats(genome_range.chrom, genome_range.start,
                                                        genome_range.end, nBins=num_bins)).astype(float)
            except Exception as e:
                import pyBigWig
                self.bw = pyBigWig.open(self.properties['file'])

                log.warning("error found while reading bigwig scores ({}).\nTrying again. Iter num: {}".
                            format(e, num_tries))
                pass
            else:
                if num_tries > 1:
                    log.warning("After {} the scores could be computed".format(num_tries))
                break
        return scores_per_bin 

# 需要導入模塊: import pyBigWig [as 別名]
# 或者: from pyBigWig import open [as 別名]
def getUCSCTracks(infile=PARAMS["filename_ucsc_encode"]):
        '''return a list of UCSC tracks from infile.'''
        tables = []
        with open(infile) as f:
            for line in f:
                if line.startswith("#"):
                    continue
                tablename = line[:-1].strip()
                if tablename == "":
                    continue
                tables.append(tablename)
        return tables

    ############################################################
    ############################################################
    ############################################################
    # import UCSC encode tracks
    ############################################################ 

# 需要導入模塊: import pyBigWig [as 別名]
# 或者: from pyBigWig import open [as 別名]
def exportUCSCEncodeTracks(infile, outfile):

        dbhandle = sqlite3.connect(PARAMS["database_name"])

        outs = open(outfile, "w")
        for tablename in getUCSCTracks():
            outs.write("track name=%s\n" % tablename)

            cc = dbhandle.cursor()
            statement = """SELECT chrom, chrostart, chroend FROM %s
            ORDER by chrom, chrostart""" % (
                tablename)
            cc.executewait(dbhandle, statement)
            for contig, start, end in cc:
                outs.write("%s\t%i\t%i\n" % (contig, start, end))
        outs.close() 

# 需要導入模塊: import pyBigWig [as 別名]
# 或者: from pyBigWig import open [as 別名]
def check_files(lst_of_files, action="r"):

	flat_lst = flatten_list(lst_of_files)
	for fil in flat_lst:
		if fil != None:
			if action == "r":
				if os.path.exists(fil):
					try: 
						with open(fil) as f:
							pass
					except:
						sys.exit("ERROR: Could not open file \"{0}\" for reading".format(fil))
				else:
					sys.exit("ERROR: File \"{0}\" does not exists".format(fil))

			elif action == "w":
				if os.path.exists(fil):
					try: 
						with open(fil, "w") as f:
							pass
					except:
						sys.exit("ERROR: Could not open file \"{0}\" for writing. Please check that you do not have the file open and that you have write permission.".format(fil)) 

# 需要導入模塊: import pyBigWig [as 別名]
# 或者: from pyBigWig import open [as 別名]
def dict_to_tab(dict_list, fname, chosen_columns, header=False):

	#Establish header
	if header == True:
		out_str = "\t".join(chosen_columns) + "\n"
	else:
		out_str = ""
	
	#Add lines
	out_str += "\n".join(["\t".join([str(line_dict[column]) for column in chosen_columns]) for line_dict in dict_list])
	
	#Add \n if out_str contains lines
	out_str += "\n" if len(out_str) > 0 else ""

	#Write file
	f = open(fname, "w")
	f.write(out_str)
	f.close()

#Quantile normalization 

# 需要導入模塊: import pyBigWig [as 別名]
# 或者: from pyBigWig import open [as 別名]
def _create_bigwig(bed_column, outpath, genome_size_dict):
    # type: (pd.Series, str, Dict[str, int]) -> None

    logging.info("Creating biwgwig " + outpath)

    bed_column = bed_column.reset_index()

    values = [float(f) for _, _, _, f in bed_column.values]
    unique_chromosomes = list(bed_column.Chromosome.drop_duplicates())
    chromosomes = list(bed_column.Chromosome)
    starts = _to_int(list(bed_column.Bin))
    ends = _to_int(list(bed_column.End + 1))

    header = [(c, int(genome_size_dict[c])) for c in unique_chromosomes]

    bw = pyBigWig.open(outpath, "w")
    bw.addHeader(header)

    bw.addEntries(chromosomes, starts, ends=ends, values=values)
    bw.close() 

# 需要導入模塊: import pyBigWig [as 別名]
# 或者: from pyBigWig import open [as 別名]
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) 

# 需要導入模塊: import pyBigWig [as 別名]
# 或者: from pyBigWig import open [as 別名]
def fetch_from_plainbed(path, chrom, start, end):
    found_chrom = False
    for line in open(path):
        fields = line.strip().split()
        if fields[0] != chrom: continue
        found_chrom = True

        cur_start, cur_end = fields[1:3]
        if int(cur_end) < start or int(cur_start) > end: continue
        yield tx_from_bedfields(fields)

    if not found_chrom:
        warning = "Didn't find chromosome {}; make sure it's formatted correctly (eg 'chr1' vs '1')".format(chrom)
        logging.warn(warning) 

# 需要導入模塊: import pyBigWig [as 別名]
# 或者: from pyBigWig import open [as 別名]
def __init__(self, path, nbins=1000, name=None):
        super().__init__(name)
        
        import pyBigWig
        self.bigwig = pyBigWig.open(path)
        self.nbins = 1000 

# 需要導入模塊: import pyBigWig [as 別名]
# 或者: from pyBigWig import open [as 別名]
def bigwig_open(bw_file, genome_file):
  """ Open the bigwig file for writing and write the header. """

  bw_out = pyBigWig.open(bw_file, 'w')

  chrom_sizes = []
  for line in open(genome_file):
    a = line.split()
    chrom_sizes.append((a[0], int(a[1])))

  bw_out.addHeader(chrom_sizes)

  return bw_out 

# 需要導入模塊: import pyBigWig [as 別名]
# 或者: from pyBigWig import open [as 別名]
def limit_segments(segments, filter_bed):
  """ Limit to segments overlapping the given BED.

    Args
     segments: list of (chrom,start,end) genomic segments
     filter_bed: BED file to filter by

    Returns:
     fsegments: list of (chrom,start,end) genomic segments
    """

  # print segments to BED
  seg_fd, seg_bed_file = tempfile.mkstemp()
  seg_bed_out = open(seg_bed_file, 'w')
  for chrom, seg_start, seg_end in segments:
    print('%s\t%d\t%d' % (chrom, seg_start, seg_end), file=seg_bed_out)
  seg_bed_out.close()

  # intersect w/ filter_bed
  fsegments = []
  p = subprocess.Popen(
      'bedtools intersect -u -a %s -b %s' % (seg_bed_file, filter_bed),
      shell=True,
      stdout=subprocess.PIPE)
  for line in p.stdout:
    a = line.decode('utf-8').split()
    chrom = a[0]
    seg_start = int(a[1])
    seg_end = int(a[2])
    fsegments.append((chrom, seg_start, seg_end))

  p.communicate()

  os.close(seg_fd)
  os.remove(seg_bed_file)

  return fsegments


################################################################################