Python BedTool.closest方法代码示例

# 需要导入模块: from pybedtools import BedTool [as 别名]
# 或者: from pybedtools.BedTool import closest [as 别名]
def feat_dist(vf, af, name):
    print "inside feat_dist"
    v = BedTool(vf)
    a = BedTool(af)
    closest = v.closest(a, D="b")
    results = dict([ (r.name, int(r[len(r.fields)-1])) for r in closest ])
    print "exiting feat_dist"
    return Series(results, name=name)

# 需要导入模块: from pybedtools import BedTool [as 别名]
# 或者: from pybedtools.BedTool import closest [as 别名]
def main():
    """
    Runs Python example from the manuscript
    """
    bedtools_dir = path.split(__file__)[0]
    snps = BedTool(path.join(bedtools_dir, '../test/data/snps.bed.gz'))
    genes = BedTool(path.join(bedtools_dir, '../test/data/hg19.gff'))

    intergenic_snps = (snps - genes)

    nearby = genes.closest(intergenic_snps, d=True, stream=True)

    for gene in nearby:
        if int(gene[-1]) < 5000:
            print gene.name

# 需要导入模块: from pybedtools import BedTool [as 别名]
# 或者: from pybedtools.BedTool import closest [as 别名]
def add_closest(aname, bname):
    a, b = BedTool(aname), BedTool(bname)

    afields = a.field_count()
    c = a.closest(b, d=True)
    get_name = gen_get_name(b, afields)

    dbed = open(BedTool._tmp(), "w")
    # keep the name and distance
    seen_by_line = collections.defaultdict(list)
    for feat in c:
        key = "\t".join(feat[:afields])
        seen_by_line[key].append([feat[-1], get_name(feat)])

    for key, dist_names in seen_by_line.iteritems():
        if len(dist_names) > 0:
            assert len(set([d[0] for d in dist_names])) == 1
        names = ",".join(sorted(set(d[1] for d in dist_names)))
        new_line = "\t".join([key] + [names] + [dist_names[0][0]])
        dbed.write(new_line + "\n")
    dbed.close()
    d = BedTool(dbed.name)
    assert len(d) == len(a)
    return d

# 需要导入模块: from pybedtools import BedTool [as 别名]
# 或者: from pybedtools.BedTool import closest [as 别名]
#!/usr/bin/python
"""
Example from the manuscript to print the names of genes that are <5000 bp away
from intergenic SNPs.  See sh_ms_example.sh for the shell script equivalent.
"""
from pybedtools import BedTool

snps = BedTool('../test/data/snps.bed.gz')
genes = BedTool('../test/data/hg19.gff')

intergenic_snps = (snps - genes)

nearby = genes.closest(intergenic_snps, d=True, stream=True)

for gene in nearby:
    if int(gene[-1]) < 5000:
        print gene.name

# 需要导入模块: from pybedtools import BedTool [as 别名]
# 或者: from pybedtools.BedTool import closest [as 别名]
def run(chipdir, refseq, filedir, DMSO, CA):
    TSS = (-200, 1000)

    a = BedTool(chipdir)
    b = a.closest(refseq, d=True)
    b.cut([9, 10, 11, 12, 13, 14, 21]).saveas(filedir + "/SRF_closest.bed")
    d = dict()
    with open(filedir + "/SRF_closest.bed") as F:
        for line in F:
            line = line.strip().split()
            chrom, start, stop = line[0:3]
            d[chrom + "\t" + start + "\t" + stop + "\t"] = "\t".join(line[3:])
    outfile = open(filedir + "/SRF_closest.rmdup.bed", "w")
    for key in d:
        if "." not in key.split():
            outfile.write(key + d[key] + "\n")
    outfile.close()
    # os.system("sort -k1,1 -k2,2n " + filedir + "/SRF_closest.rmdup.bed > " + filedir + "/SRF_closest.rmdup.sorted.bed")
    a = BedTool(filedir + "/SRF_closest.rmdup.bed")
    a.sort().saveas(filedir + "/SRF_closest.rmdup.sorted.bed")
    outfile = open(filedir + "/SRF.TSS.bed", "w")
    outfile2 = open(filedir + "/SRF.gene.bed", "w")
    with open(filedir + "/SRF_closest.rmdup.sorted.bed") as F:
        for line in F:
            chrom, start, stop, gene, number, strand, distance = line.strip().split()
            if int(stop) - int(start) > 2000 and int(distance) > 10000:
                if strand is "+":
                    outfile.write(chrom + "\t" + str(int(start) + TSS[0]) + "\t" + str(int(start) + TSS[1]) + "\n")
                    outfile2.write(chrom + "\t" + str(int(start) + TSS[1]) + "\t" + stop + "\n")
                else:
                    outfile.write(chrom + "\t" + str(int(stop) - TSS[1]) + "\t" + str(int(stop) - TSS[0]) + "\n")
                    outfile2.write(chrom + "\t" + start + "\t" + str(int(stop) - TSS[1]) + "\n")
    outfile.close()
    outfile2.close()
    a = BedTool(filedir + "/SRF.TSS.bed")
    a.sort().saveas(filedir + "/SRF.TSS.bed")
    a = BedTool(filedir + "/SRF.gene.bed")
    a.sort().saveas(filedir + "/SRF.gene.bed")

    TSS = filedir + "/SRF.TSS.bed"
    genes = filedir + "/SRF.gene.bed"

    os.system("bedtools map -a " + genes + " -b " + DMSO + " -c 4 -o sum > " + filedir + "/DMSO.genes.bed")
    os.system("bedtools map -a " + TSS + " -b " + DMSO + " -c 4 -o sum > " + filedir + "/DMSO.TSS.bed")
    os.system("bedtools map -a " + genes + " -b " + CA + " -c 4 -o sum > " + filedir + "/CA.genes.bed")
    os.system("bedtools map -a " + TSS + " -b " + CA + " -c 4 -o sum > " + filedir + "/CA.TSS.bed")

    TRx = list()
    TRy = list()
    expressionlist = list()

    with open(filedir + "/DMSO.genes.bed") as a, open(filedir + "/DMSO.TSS.bed") as b, open(
        filedir + "/CA.genes.bed"
    ) as c, open(filedir + "/CA.TSS.bed") as d:
        for line in a:
            bline = b.readline().strip().split()[-1]
            cline = c.readline().strip().split()[-1]
            dline = d.readline().strip().split()[-1]
            if line.strip().split()[-1] is ".":
                DMSOgene = 0.0
            else:
                DMSOgene = float(line.strip().split()[-1])
            if bline is ".":
                DMSOTSS = 0.0
            else:
                DMSOTSS = float(bline)
            if cline is ".":
                CAgene = 0.0
            else:
                CAgene = float(cline)
            if dline is ".":
                CATSS = 0.0
            else:
                CATSS = float(dline)
            if DMSOgene == 0.0:
                TRx.append(0.0)
            else:
                TRx.append((DMSOTSS / DMSOgene))
            if CAgene == 0.0:
                TRy.append(0.0)
            else:
                TRy.append((CATSS / CAgene))
            expressionlist.append((np.log2(DMSOgene) + np.log2(CAgene)) / 2.0)

    F6 = plt.figure()
    ax1 = F6.add_subplot(111)
    xy = np.vstack([TRx, TRy])
    z = gaussian_kde(xy)(xy)
    ax1.scatter(TRx, TRy, c=z, edgecolor="")
    # ax1.scatter(TRx2,TRy2,c='red',edgecolor="",s=expressionlist2)
    ax1.set_title("Pausing Index")
    ax1.set_ylabel("CA")
    ax1.set_xlabel("DMSO")
    ax1.get_xaxis().tick_bottom()
    ax1.get_yaxis().tick_left()
    # ax1.plot([0,1/slope1],[intercept1,1],color = 'r')
    ax1.set_xlim([0, 20])
    ax1.set_ylim([0, 20])
    ax1.plot([0, 50.0], [0, 50.0], color="k")
    # ax1.text(8,18, "Pearson = " + str(pearsons)[0:5])
#.........这里部分代码省略.........