use GATK to split reads with Ns in the CIGAR string, hard clipping regions
    that end up in introns
13322
def gatk_splitreads(data):
    
    broad_runner = broad.runner_from_config(dd.get_config(data))
    ref_file = dd.get_ref_file(data)
    deduped_bam = dd.get_deduped_bam(data)
    base, ext = os.path.splitext(deduped_bam)
    split_bam = base + ".splitN" + ext
    if file_exists(split_bam):
        data = dd.set_split_bam(data, split_bam)
        return data
    gatk_type = broad_runner.gatk_type()
    with file_transaction(data, split_bam) as tx_split_bam:
        params = ["-T", "SplitNCigarReads",
                  "-R", ref_file,
                  "-I", deduped_bam]
        if gatk_type == "gatk4":
            params += ["--output", tx_split_bam]
        else:
            params += ["-rf", "ReassignOneMappingQuality",
                       "-RMQF", "255",
                       "-RMQT", "60",
                       "-rf", "UnmappedRead",
                       "-U", "ALLOW_N_CIGAR_READS",
                       "-o", tx_split_bam]
            if dd.get_quality_format(data) == "illumina":
                params += ["--fix_misencoded_quality_scores", "-fixMisencodedQuals"]
        broad_runner.run_gatk(params)
    bam.index(split_bam, dd.get_config(data))
    data = dd.set_split_bam(data, split_bam)
    return data
767557
def _handle_splits(_str):
    """"""
    _str = _str.replace('/', '-')
    _tmp_dict = {}

    if '-' in _str:
        start, stop = _str.split('-')
        if _check_number(start):
            start = regex.sub(r'[0-9]+\?*', start, stop)
        elif _check_number(stop):
            stop = regex.sub(r'[0-9]+\?*', stop, start)
    else:
        start = _str
        stop = _str
    _tmp_dict['start_raw'] = start
    _tmp_dict['stop_raw'] = stop

    _tmp_dict['start_epoch'] = _get_epoch(start)
    _tmp_dict['stop_epoch'] = _get_epoch(stop)

    return _tmp_dict
1297071
def from_string(string):
        
        lines = list(clean_lines(string.splitlines()))

        def input_mode(line):
            if line[0] == "&":
                return ("sections", line[1:].lower())
            elif "ATOMIC_SPECIES" in line:
                return ("pseudo", )
            elif "K_POINTS" in line:
                return ("kpoints", line.split("{")[1][:-1])
            elif "CELL_PARAMETERS" in line or "ATOMIC_POSITIONS" in line:
                return ("structure", line.split("{")[1][:-1])
            elif line == "/":
                return None
            else:
                return mode

        sections = {"control": {}, "system": {}, "electrons": {}, 
                    "ions": {}, "cell":{}}
        pseudo = {}
        pseudo_index = 0
        lattice = []
        species = []
        coords = []
        structure = None
        site_properties = {"pseudo":[]}
        mode = None
        for line in lines:
            mode = input_mode(line)
            if mode == None:
                pass
            elif mode[0] == "sections":
                section = mode[1]
                m = re.match(r'(\w+)\(?(\d*?)\)?\s*=\s*(.*)', line)
                if m:
                    key = m.group(1).strip()
                    key_ = m.group(2).strip()
                    val = m.group(3).strip()
                    if key_ != "":
                        if sections[section].get(key, None) == None:
                            val_ = [0.0]*20 # MAX NTYP DEFINITION
                            val_[int(key_)-1] = PWInput.proc_val(key, val)
                            sections[section][key] = val_

                            site_properties[key] = []
                        else:
                            sections[section][key][int(key_)-1] = PWInput.proc_val(key, val) 
                    else:
                        sections[section][key] = PWInput.proc_val(key, val)

            elif mode[0] == "pseudo":
                m = re.match(r'(\w+)\s+(\d*.\d*)\s+(.*)', line)
                if m:
                    pseudo[m.group(1).strip()] = {}
                    pseudo[m.group(1).strip()]["index"] = pseudo_index
                    pseudo[m.group(1).strip()]["pseudopot"] = m.group(3).strip()
                    pseudo_index += 1
            elif mode[0] == "kpoints":
                m = re.match(r'(\d+)\s+(\d+)\s+(\d+)\s+(\d+)\s+(\d+)\s+(\d+)', line)
                if m:
                    kpoints_grid = (int(m.group(1)), int(m.group(2)), int(m.group(3)))
                    kpoints_shift = (int(m.group(4)), int(m.group(5)), int(m.group(6)))
                else:
                    kpoints_mode = mode[1]
            elif mode[0] == "structure":
                m_l = re.match(r'(-?\d+\.?\d*)\s+(-?\d+\.?\d*)\s+(-?\d+\.?\d*)', line)
                m_p = re.match(r'(\w+)\s+(-?\d+\.\d*)\s+(-?\d+\.?\d*)\s+(-?\d+\.?\d*)', line)
                if m_l:
                    lattice += [ float(m_l.group(1)), float(m_l.group(2)), float(m_l.group(3)) ]
                elif m_p:
                    site_properties["pseudo"].append(pseudo[m_p.group(1)]["pseudopot"])
                    species += [pseudo[m_p.group(1)]["pseudopot"].split(".")[0]]
                    coords += [[float(m_p.group(2)), float(m_p.group(3)), float(m_p.group(4))]]

                    for k, v in site_properties.items():
                        if k != "pseudo":
                            site_properties[k].append(sections['system'][k][pseudo[m_p.group(1)]["index"]])
                if mode[1] == "angstrom":
                    coords_are_cartesian = True
                elif mode[1] == "crystal":
                    coords_are_cartesian = False

        structure = Structure(Lattice(lattice), species, coords, 
                              coords_are_cartesian=coords_are_cartesian,
                              site_properties=site_properties)
        return PWInput(structure=structure, control=sections["control"],
                       system=sections["system"], electrons=sections["electrons"], 
                       ions=sections["ions"], cell=sections["cell"], kpoints_mode=kpoints_mode,
                       kpoints_grid=kpoints_grid, kpoints_shift=kpoints_shift)
1172850
def remove_read_more_sep(self, raw_content):
        
        if self._read_more_exp is None:
            return raw_content

        sp = self._read_more_exp.split(raw_content, maxsplit=1)
        if len(sp) == 2 and sp[0]:
            result = '\n\n'.join((sp[0].rstrip(), sp[1].lstrip()))
        else:
            result = raw_content
        return result
1001898
def parse(self, contents):
        
        i = 0
        for text in contents.split(self.delim):
            if not len(text.strip()):
                continue
            words = text.split()
            char_offsets = [0] + [
                int(_) for _ in np.cumsum([len(x) + 1 for x in words])[:-1]
            ]
            text = " ".join(words)
            yield {
                "text": text,
                "words": words,
                "pos_tags": [""] * len(words),
                "ner_tags": [""] * len(words),
                "lemmas": [""] * len(words),
                "dep_parents": [0] * len(words),
                "dep_labels": [""] * len(words),
                "char_offsets": char_offsets,
                "abs_char_offsets": char_offsets,
            }
            i += 1
1203383
def split_sequence(seq, n):
    """"""
    tokens = []
    while seq:
        tokens.append(seq[:n])
        seq = seq[n:]
    return tokens
1180194
def get_gene2aart(gene2section2gos, sec2chr):
        """"""
        geneid2str = {}
        for geneid, section2gos_gene in gene2section2gos.items():
            letters = [abc if s in section2gos_gene else "." for s, abc in sec2chr.items()]
            geneid2str[geneid] = "".join(letters)
        return geneid2str
894775
def _sam_to_soft_clipped(self, sam):
        ''''''
        if sam.is_unmapped:
            raise Error('Cannot get soft clip info from an unmapped read')

        if sam.cigar is None or len(sam.cigar) == 0:
            return False, False
        return (sam.cigar[0][0] == 4, sam.cigar[-1][0] == 4)
1245419
function deepSplitString (input, separators) {
  if (input === null || typeof input === 'undefined') {
    return input
  }
  if (separators.length === 0) {
    return input
  }

  const sep = separators.shift()
  let output = input

  if (typeof input === 'string') {
    output = splitString(input, sep)
  } else if (Array.isArray(input)) {
    output = input.map(i => splitString(i, sep))
  } else if (typeof input === 'object') {
    output = {}
    Object.keys(input).forEach(key => {
      output[key] = splitString(input[key], sep)
    })
  }

  if (separators.length > 0) {
    return deepSplitString(output, separators)
  } else {
    return output
  }
}
893073
def _piped_input_cl(data, region, tmp_dir, out_base_file, prep_params):
    
    return data["work_bam"], _gatk_extract_reads_cl(data, region, prep_params, tmp_dir)
751949
def clean_readme(fname):
    """"""
    with codecs.open(fname, 'r', 'utf-8') as f:
        return ''.join(
            re.sub(r':\w+:`([^`]+?)( <[^<>]+>)?`', r'``\1``', line)
            for line in f
            if not (line.startswith('.. currentmodule') or line.startswith('.. toctree'))
        )
46697
def _split_tidy(self, string, maxsplit=None):
        """"""

        if maxsplit is None:
            return string.rstrip("\n").split("\t")
        else:
            return string.rstrip("\n").split("\t", maxsplit)
68580
def _split_by_ready_regions(ext, file_key, dir_ext_fn):
    
    def _sort_by_size(region_w_bams):
        region, _ = region_w_bams
        _, start, end = region
        return end - start
    def _assign_bams_to_regions(data):
        """Ensure BAMs aligned with input regions, either global or individual.
        """
        for i, region in enumerate(data["region"]):
            work_bams = []
            for xs in data["region_bams"]:
                if len(xs) == 1:
                    work_bams.append(xs[0])
                else:
                    work_bams.append(xs[i])
            for work_bam in work_bams:
                assert os.path.exists(work_bam), work_bam
            yield region, work_bams
    def _do_work(data):
        if "region" in data:
            name = data["group"][0] if "group" in data else data["description"]
            out_dir = os.path.join(data["dirs"]["work"], dir_ext_fn(data))
            out_file = os.path.join(out_dir, "%s%s" % (name, ext))
            assert isinstance(data["region"], (list, tuple))
            out_parts = []
            for r, work_bams in sorted(_assign_bams_to_regions(data), key=_sort_by_size, reverse=True):
                out_region_dir = os.path.join(out_dir, r[0])
                out_region_file = os.path.join(out_region_dir,
                                               "%s-%s%s" % (name, pregion.to_safestr(r), ext))
                out_parts.append((r, work_bams, out_region_file))
            return out_file, out_parts
        else:
            return None, []
    return _do_work
1150902
def main(args):
   """"""
   txome = Transcriptome()
   #read the reference gpd if one is gven
   if args.reference: 
      rinf = None
      if re.search('\.gz$',args.reference):
         rinf = gzip.open(args.reference)
      else:
         rinf = open(args.reference)
      sys.stderr.write("Reading in reference\n")
      z = 0
      # populate txome with reference transcripts for each chromosome
      for line in rinf:
         z += 1
         gpd = GPD(line)
         gpd.set_payload(z)
         if z%100 == 0:  sys.stderr.write(str(z)+"          \r")
         txome.add_transcript(gpd)
      rinf.close()
      sys.stderr.write(str(z)+"          \r")
      sys.stderr.write("\n")
   txome.sort_transcripts()
   sys.stderr.write("Buffering mappings\n")
   inf = sys.stdin
   if args.input != '-':
      if re.search('\.gz$',args.input):
         inf = gzip.open(args.input)
      else:
         inf = open(args.input)
   tof = gzip.open(args.tempdir+'/reads.gpd.gz','w')
   for line in inf: tof.write(line.rstrip()+"\n")
   tof.close()
  
   sys.stderr.write("1. Process by overlapping locus.\n")
   annotated_reads(txome,gzip.open(args.tempdir+'/reads.gpd.gz'),'initial',args)

   #of.close()

   """Now sequences we should be able to annotate as partial have
      consensus transcripts.  Lets put those to annotation scrutiny again.
   """
   sys.stderr.write("2. Reannotated partially annotated by overlapping locus.\n")
   annotated_reads(txome,gzip.open(args.tempdir+'/initial/partial_annotated_multiexon.sorted.gpd.gz'),'repartial',args)


   txs = {}
   for transcript in txome.transcripts:
      txs[transcript.name] = transcript
   detected = {}
   tinf = gzip.open(args.tempdir+'/initial/annotated.txt.gz')
   for line in tinf:
      f = line.rstrip().split("\t")
      detected[f[1]] = True
   tinf.close()
   tinf = gzip.open(args.tempdir+'/repartial/annotated.txt.gz')
   remove_partial = {}
   for line in tinf:
      f = line.rstrip().split("\t")
      detected[f[1]] = True
      remove_partial[f[1]] = True
   tinf.close()
   tof = gzip.open(args.tempdir+'/candidate.gpd.gz','w')
   for name in detected:
      tof.write(txs[name].get_gpd_line()+"\n")
   tinf = gzip.open(args.tempdir+'/initial/partial_annotated_multiexon.sorted.gpd.gz')
   for line in tinf:
      f = line.rstrip().split("\t")
      if f[1] not in remove_partial:
         tof.write(line)
   tinf = gzip.open(args.tempdir+'/initial/unannotated_singleexon.sorted.gpd.gz')
   for line in tinf:
      f = line.rstrip().split("\t")
      f[0] = str(uuid.uuid4())
      tof.write("\t".join(f)+"\n")
   tinf.close()
   tof.close()
   sort_gpd(args.tempdir+'/candidate.gpd.gz',args.tempdir+'/candidate.sorted.gpd.gz',args)
   """We can ignore the partial annotations that have been detected"""
   ntxome = Transcriptome()
   tinf = gzip.open(args.tempdir+'/candidate.sorted.gpd.gz')
   for line in tinf: ntxome.add_transcript(GPD(line))     
   annotated_reads(ntxome,gzip.open(args.tempdir+'/initial/unannotated_multiexon.sorted.gpd.gz'),'unannotated',args)
   """now we know which unannotated reads actually have annotations"""
   tinf.close()

   tof = gzip.open(args.tempdir+'/final.gpd.gz','w')
   tinf = gzip.open(args.tempdir+'/candidate.sorted.gpd.gz')   
   for line in tinf: tof.write(line)
   tinf.close()
   tinf = gzip.open(args.tempdir+'/unannotated/unannotated_multiexon.sorted.gpd.gz')   
   for line in tinf: tof.write(line)
   tinf.close()
   tinf = gzip.open(args.tempdir+'/unannotated/partial_annotated_multiexon.sorted.gpd.gz')   
   for line in tinf: tof.write(line)
   tinf.close()
   tof.close()
   sort_gpd(args.tempdir+'/final.gpd.gz',args.tempdir+'/final.sorted.gpd.gz',args)
   """Prepare outputs"""
   of = sys.stdout
   if args.output:
      if re.search('\.gz$',args.output):
         of = gzip.open(args.output,'w')
      else:
         of = open(args.output,'w')
   tinf = gzip.open(args.tempdir+'/final.sorted.gpd.gz')
   for line in tinf: of.write(line)
   tinf.close()
   of.close()
958671
def _init_exons(self):
        
        exon_starts = [self.chrom_start + int(s)
                       for s in self.bed_tuple.blockStarts.strip(',').split(',')]
        exon_sizes = list(map(int, self.bed_tuple.blockSizes.strip(',').split(',')))

        # get chromosome intervals
        exons = [(exon_starts[i], exon_starts[i] + exon_sizes[i])
                 for i in range(len(exon_starts))]
        no_utr_exons = self._filter_utr(exons)
        self.exons = no_utr_exons
        self.exon_lens = [e[1] - e[0] for e in self.exons]
        self.num_exons = len(self.exons)
        self.cds_len = sum(self.exon_lens)
        self.five_ss_len = 2*(self.num_exons-1)
        self.three_ss_len = 2*(self.num_exons-1)
        self._init_splice_site_pos()
1107006
def split_denovos(denovo_path, temp_dir):
    
    
    # open the de novos, drop the header line, then sort them (which will be by
    # HGNC symbol, as the first element of each line)
    with open(denovo_path, "r") as handle:
        lines = handle.readlines()
        header = lines.pop(0)
    
    basename = os.path.basename(denovo_path)
    
    # only include genes with 2+ missense SNVs
    counts = count_missense_per_gene(lines)
    counts = dict((k, v) for k, v in counts.items() if v > 1 )
    
    genes = set([])
    for line in sorted(lines):
        gene = line.split("\t")[0]
        
        # open a new output file whenever we hit a different gene
        if gene not in genes and gene in counts:
            genes.add(gene)
            path = os.path.join(temp_dir, "tmp.{}.txt".format(len(genes)))
            output = open(path, "w")
            output.write(header)
        
        if gene in counts:
            output.write(line)
    
    return len(genes)
338920
def _run_split_on_punc(self, text):
        """"""
        chars = list(text)
        i = 0
        start_new_word = True
        output = []
        while i < len(chars):
            char = chars[i]
            if self._is_punctuation(char):
                output.append([char])
                start_new_word = True
            else:
                if start_new_word:
                    output.append([])
                start_new_word = False
                output[-1].append(char)
            i += 1

        return [''.join(x) for x in output]