peakcaller.bak.py 文件源码

def peakcaller(tmp_dir, out_dir, gtf_fp, unique_only=False, with_replicates=False, with_control=False, unstranded=False):
    """DOCSTRING
    Args:
    Returns:
    """
    # file handlers
    mbam = pysam.Samfile(os.path.join(out_dir, 'realigned.sorted.bam'),'rb')
    ubam = pysam.Samfile(os.path.join(tmp_dir, 'unique.sorted.bam'),'rb')
    bam_dict = {'ubam.ip':[ubam], 'mbam.ip':[mbam]}
    if unique_only:
        ofile = open(os.path.join(out_dir, 'narrow_peaks.unique.bed'), 'w')
    else:
        ofile = open(os.path.join(out_dir, 'narrow_peaks.combined.bed'), 'w')

    # read in GTF
    gene_annot = read_gtf(gtf_fp)

    # iteratively call peaks in each gene
    peak_counter = 0
    for gene_name in tqdm(gene_annot):
        gene = gene_annot[gene_name]
        BED = call_gene_peak(bam_dict, gene, 
            unique_only=unique_only, with_control=with_control, 
            unstranded=unstranded)
        ofile.write(BED)
        #print BED
        peak_counter += len(BED.split('\n'))
    ofile.close()
    logger.info('called %i peaks'%peak_counter)
    return