Python BedTool.subtract方法代码示例

# 需要导入模块: from pybedtools import BedTool [as 别名]
# 或者: from pybedtools.BedTool import subtract [as 别名]
def subtract_bed_sd(bed_name, bed_filter):
    """REMOVES regions of annotation of interest that overlap with
    segmental duplications
    """
    pybedtools.set_tempdir('/sc/orga/scratch/richtf01')
    if not os.path.isfile(bed_name + '.noRmsk.noSD.bed'):
        bed = BedTool(bed_name + '.noRmsk.bed')
        print "Removing calls in seg dup from " + bed_name + "..."
        bed_no_overlap = bed.subtract(bed_filter)
        bed_no_overlap.saveas(bed_name + '.noRmsk.noSD.bed')
        print bed_name + " done!"
    else:
        print bed_name + " Seg dup calls already removed"

# 需要导入模块: from pybedtools import BedTool [as 别名]
# 或者: from pybedtools.BedTool import subtract [as 别名]
def subtract_bed_rmsk(bed_name, bed_filter):
    """REMOVES regions of annotation of interest that overlap with
    repeat-masked regions
    """
    pybedtools.set_tempdir('/sc/orga/scratch/richtf01')
    if not os.path.isfile(bed_name + '.noRmsk.bed'):
        bed = BedTool(bed_name + '.bed') # .merged.sorted
        print "Removing calls in rmsk from " + bed_name + "..."
        bed_no_overlap = bed.subtract(bed_filter)
        bed_no_overlap.saveas(bed_name + '.noRmsk.bed')
        print bed_name + " done!"
    else:
        print bed_name + " rmsk calls already removed"

# 需要导入模块: from pybedtools import BedTool [as 别名]
# 或者: from pybedtools.BedTool import subtract [as 别名]
def remainder_size(bed_file, original_bed):
    """
    Calculate the number of bases in the small panel not included in the truth regions.

    The amount of the panel that overlaps with the truth regions is calculated. This gives an idea of how well the
    panel is represented in the reference sample and whether it is an accurate reflection of the accuracy of the
    process.

    :param bed_file: File path to the specific panel bed file
    :type bed_file: String
    :param original_bed: Path to the original BED file that was given as an argument
    :type original_bed: String
    :return: Total length of regions in remainder BED file
    :rtype: Int
    """
    print('Calculating remainder')
    print(bed_file)
    print(original_bed)
    truth_region_bed = bed_file.replace('_1based', '')
    print(truth_region_bed)

    truth_region_tool = BedTool(truth_region_bed)

    original_tool = BedTool(original_bed)

    remainder_name = bed_file.replace('_truth_regions_1based', '_remainder')
    remainder_bed = original_tool.subtract(truth_region_tool)
    remainder_bed.saveas(remainder_name)

    f = open(remainder_name, 'r')
    regions = [line.strip('\n') for line in f.readlines()]
    f.close()

    total_length = 0

    for region in regions:
        fields = region.split('\t')
        start = int(fields[1])
        end = int(fields[2])

        length = end - start
        total_length += length

    return total_length

# 需要导入模块: from pybedtools import BedTool [as 别名]
# 或者: from pybedtools.BedTool import subtract [as 别名]
    def compare_bed(self, bed_file, orig_from_string, pp_bed):
        d = '/home/bioinfo/Natalie/wc/genes/compare_beds/'

        original_bed = BedTool(bed_file, from_string=orig_from_string)
        original_name = os.path.basename(bed_file)
        saveas = d + original_name
        original_bed.saveas(saveas)
        saveas_pp = d + 'pp_' + original_name
        pp_bed.saveas(saveas_pp)

        missing_from_pp = original_bed.subtract(pp_bed)
        print 'Missing from exported BED'
        print missing_from_pp
        saveas_missing_pp = d + 'missing_pp_' + original_name
        missing_from_pp.saveas(saveas_missing_pp)
        missing_from_orig = pp_bed.subtract(original_bed)
        print 'Missing from original BED'
        print missing_from_orig
        saveas_missing_orig = d+ 'missing_orig_' + original_name
        missing_from_orig.saveas(saveas_missing_orig)

# 需要导入模块: from pybedtools import BedTool [as 别名]
# 或者: from pybedtools.BedTool import subtract [as 别名]
#        line=f.readline()
#        if line=='':
#            break
#        lineSplit=line.split('\t')
#        features=lineSplit[8].split(';')
#        if lineSplit[2]=='exon':
#            print >>cuffToBed, lineSplit[0]+'\t'+str(int(lineSplit[3])-1)+'\t'+\
#                  lineSplit[4]+'\t'+features[1][16:len(features[1])-1]\
#                  +'\t'+lineSplit[5]+'\t'+lineSplit[6]
#    f.close()


miRNA=BedTool('hsa.bed')
exome=BedTool('exome.bed')
print 'bed files loaded'
intergenic=miRNA.subtract(exome)

windowed=intergenic.slop(b=50000,g='human.hg19.genome')


cuffToBed('/raid6/ogan/cuffLinksOutput/'+str(fileno)+'/transcripts.gtf','/raid6/ogan/cuffBedT/'+str(fileno))
#cuffToBed2('/raid6/ogan/cuffLinksOutput/'+str(fileno)+'/transcripts.gtf','/raid6/ogan/cuffBedE/'+str(fileno))

cuff=BedTool('/raid6/ogan/cuffBedT/'+str(fileno))
#cuff2=BedTool('/raid6/ogan/cuffBedT/'+str(fileno))
cuffmiRNA=cuff.intersect(intergenic,wa=True,s=True,wb=True)
#cuff2miRNA=cuff2.intersect(intergenic,wa=True,s=True,wb=True)



cuffmiRNA.saveas('/raid6/ogan/cuffInt/cuffmiRNA'+str(fileno))

# 需要导入模块: from pybedtools import BedTool [as 别名]
# 或者: from pybedtools.BedTool import subtract [as 别名]
                         # -1 to account for the center nucleotide!
                         r=int(args.core_length / 2) +
                         (args.core_length % 2) - 1,
                         genome=args.genome_id).each(offset_zero_by_one).saveas(pos_core_bed_fn)

flanks_upstream, flanks_downstream = get_flanks(cores)
get_seqs(cores, flanks_upstream, flanks_downstream, pos_seq_fa_fn)

# prepare negative sites if requested
if args.negative_site_candidate_regions_fn:
    # get negative candidate regions
    negative_site_candidate_regions = BedTool(
        args.negative_site_candidate_regions_fn)
    # remove input binding sites from negative candidate regions
    processed_negative_site_candidate_regions = negative_site_candidate_regions.subtract(
        bsites,
        s=True).saveas()

    # create negative core sites by placing within candidate regions
    logging.info("preparing negative instances")
    logging.info("starting from " + str(cores.count()) + " positive cores")
    if args.chromosome_limits:
        logging.debug("using chromosome_limits " + args.chromosome_limits)
        neg_cores = cores.shuffle(
            g=args.chromosome_limits,
            chrom=True,
            incl=processed_negative_site_candidate_regions.fn,
            noOverlapping=True).each(prefix_neg).saveas(neg_core_bed_fn)
        logging.info("derived negative cores: " + str(neg_cores.count()))
        neg_fup, neg_fdown = get_flanks(neg_cores)
        get_seqs(neg_cores, neg_fup, neg_fdown, neg_seq_fa_fn)