def ref_build_and_muscle_chunk(data, sample):
<EOL>regions = bedtools_merge(data, sample).strip().split("<STR_LIT:\n>")<EOL>nregions = len(regions)<EOL>chunksize = (nregions / <NUM_LIT:10>) + (nregions % <NUM_LIT:10>)<EOL>LOGGER.debug("<STR_LIT>".format(nregions, chunksize))<EOL>idx = <NUM_LIT:0><EOL>tmpfile = os.path.join(data.tmpdir, sample.name+"<STR_LIT>")<EOL>for i in range(<NUM_LIT:11>):<EOL><INDENT>if os.path.exists(tmpfile.format(i)):<EOL><INDENT>os.remove(tmpfile.format(i))<EOL><DEDENT><DEDENT>fopen = open<EOL>if data.paramsdict["<STR_LIT>"] == "<STR_LIT>":<EOL><INDENT>tmpfile = os.path.join(data.dirs.clusts, sample.name+"<STR_LIT>")<EOL>fopen = gzip.open<EOL><DEDENT>samfile = pysam.AlignmentFile(sample.files.mapped_reads, '<STR_LIT:rb>')<EOL>clusts = []<EOL>nclusts = <NUM_LIT:0><EOL>for region in regions:<EOL><INDENT>chrom, pos1, pos2 = region.split()<EOL>try:<EOL><INDENT>if "<STR_LIT>" in data.paramsdict["<STR_LIT>"]:<EOL><INDENT>clust = fetch_cluster_pairs(data, samfile, chrom, int(pos1), int(pos2))<EOL><DEDENT>else:<EOL><INDENT>clust = fetch_cluster_se(data, samfile, chrom, int(pos1), int(pos2))<EOL><DEDENT><DEDENT>except IndexError as inst:<EOL><INDENT>LOGGER.error("<STR_LIT>".format(chrom, pos1, pos2))<EOL>continue<EOL><DEDENT>if clust:<EOL><INDENT>clusts.append("<STR_LIT:\n>".join(clust))<EOL>nclusts += <NUM_LIT:1><EOL>if nclusts == chunksize:<EOL><INDENT>tmphandle = tmpfile.format(idx)<EOL>with fopen(tmphandle, '<STR_LIT:a>') as tmp:<EOL><INDENT>tmp.write("<STR_LIT>".join(clusts)+"<STR_LIT>")<EOL><DEDENT>idx += <NUM_LIT:1><EOL>nclusts = <NUM_LIT:0><EOL>clusts = []<EOL><DEDENT><DEDENT><DEDENT>if clusts:<EOL><INDENT>with fopen(tmpfile.format(idx), '<STR_LIT:a>') as tmp:<EOL><INDENT>tmp.write("<STR_LIT>".join(clusts)+"<STR_LIT>")<EOL><DEDENT>clusts = []<EOL><DEDENT>if not data.paramsdict["<STR_LIT>"] == "<STR_LIT>":<EOL><INDENT>chunkfiles = glob.glob(os.path.join(data.tmpdir, sample.name+"<STR_LIT>"))<EOL>LOGGER.info("<STR_LIT>", chunkfiles)<EOL><DEDENT>samfile.close()<EOL>
1. Run bedtools to get all overlapping regions
2. Parse out reads from regions using pysam and dump into chunk files. 
   We measure it out to create 10 chunk files per sample. 
3. If we really wanted to speed this up, though it is pretty fast already, 
   we could parallelize it since we can easily break the regions into 
   a list of chunks.