Parse sam files with aligned sequences
429781
def parse_align_file(file_in):
    
    loc_id = 1
    bedfile_clusters = ""
    bamfile = pybedtools.BedTool(file_in)
    bed = pybedtools.BedTool.bam_to_bed(bamfile)
    for c, start, end, name, q, strand in bed:
        loc_id += 1
        bedfile_clusters += "%s\t%s\t%s\t%s\t%s\t%s\n" % \
                            (c, start, end, name, loc_id, strand)
    return bedfile_clusters
239172
def parseArgv():
    
    parser = OptionParser()
    parser.usage = r""" cat <TEXT> | %prog [--unit <UNIT>] [--output <SA_FILE>]

Create the suffix array of TEXT with the processing UNIT and optionally store it in SA_FILE for subsequent use.
UNIT may be set to 'byte', 'character' (given an encoding with the --encoding option) or 'word', which is the default.
"""

    parser.add_option("-i", "--input",
                      action="store", type="string", dest="input",
                      default=False,
                      help="Path of the file containing the input text. When '-' is given, read the standard input (default). If the path ends with '.gz', reads the decompressed file.")

    parser.add_option("-o", "--output",
                      action="store", type="string", dest="output",
                      default=False,
                      help="Store the suffix array of the input to the file OUTPUT. When '-' is given, writes to the standard output. If the filename ends with '.gz', the suffix array will be stored  compressed.")

    parser.add_option("", "--load",
                      action="store", type="string", dest="SAFile",
                      default=False,
                      help="Load a suffix array from SAFILE, this option and --input are mutually exclusive.")

    parser.add_option("-u", "--unit",
                      action="store", type="string", dest="unit",
                      default=DEFAULT_UNIT_STR,
                      help="Processing unit used for the creation of the suffix array." + \
                           'Possible values are "byte", "character" and "word". Default is "%s".' % DEFAULT_UNIT_STR + \
                           "This option is ignored when the suffix array is loaded from SAFILE." + \
                           'For characters, the input is decoded according to the encoding set via the option --encoding.')

    parser.add_option("-e", "--encoding",
                      action="store", type="string", dest="encoding",
                      default=DEFAULT_ENCODING,
                      help="Encoding of the input. This information is required only when processing characters. Default is '%s'." % DEFAULT_ENCODING)

    parser.add_option("-p", "--print",
                      action="store_true", dest="printSA",
                      default=False,
                      help="Prints the suffix array in a human readable format to the standard error output.")

    parser.add_option("", "--verbose",
                      action="store_true", dest="verbose",
                      default=False,
                      help="Prints more information.")

    parser.add_option("", "--no-lcps",
                      action="store_true", dest="noLCPs",
                      default=False,
                      help="Switch off the computation of LCPs. By doing so, the find functions are unusable.")

    (options, args) = parser.parse_args(_argv)
    strings = args[1:]
    return (options, strings)
991145
def bamtag(sam):
    ''' 
    '''
    from pysam import AlignmentFile

    start_time = time.time()

    sam_file = open_bamfile(sam)
    out_file = AlignmentFile("-", "wh", template=sam_file)
    track = sam_file.fetch(until_eof=True)

    # peek at first alignment to determine the annotations
    if is_python3():
        queryalignment = next(track)
    else:
        queryalignment = track.next()
    annotations = detect_alignment_annotations(queryalignment)
    track = itertools.chain([queryalignment], track)

    re_string = construct_transformed_regex(annotations)
    parser_re = re.compile(re_string)

    for count, aln in enumerate(track, start=1):
        if count and not count % 1000000:
            logger.info("Processed %d alignments." % count)

        match = parser_re.match(aln.qname)
        tags = aln.tags

        if "cellular" in annotations:
            aln.tags += [('XC', match.group('CB'))]
        if "molecular" in annotations:
            aln.tags += [('RX', match.group('MB'))]
        if "sample" in annotations:
            aln.tags += [('XS', match.group('SB'))]

        out_file.write(aln)

    total_time = time.time() - start_time
    logger.info('BAM tag conversion done - {:.3}s, {:,} alns/min'.format(total_time, int(60. * count / total_time)))
    logger.info("Processed %d alignments." % count)
445138
def SAMflags(x):
    
    flags=[]

    if x & 1:
        l="1: Read paired"
    else:
        l="0: Read unpaired"
    flags.append(l)

    if x & 2 :
        l="1: Read mapped in proper pair"
    else:
        l="0: Read not mapped in proper pair"
    flags.append(l)

    if x & 4 :
        l="1: Read unmapped"
    else:
        l="0: Read mapped"
    flags.append(l)

    if x & 8 :
        l="1: Mate unmapped"
    else:
        l="0: Mate mapped"
    flags.append(l)

    if x & 16 :
        l="1: Read reverse strand"
    else:
        l="0: Read direct strand"
    flags.append(l)

    if x & 32 :
        l="1: Mate reverse strand"
    else:
        l="0: Mate direct strand"
    flags.append(l)

    if x & 64 :
        l="1: First in pair"
    else:
        l="0: Second in pair"
    flags.append(l)

    if x & 128 :
        l="1: Second in pair"
    else:
        l="0: First in pair"
    flags.append(l)

    if x & 256 :
        l="1: Not primary alignment"
    else:
        l="0: Primary alignment"
    flags.append(l)

    if x & 512 :
        l="1: Read fails platform/vendor quality checks"
    else:
        l="0: Read passes platform/vendor quality checks"
    flags.append(l)

    if x & 1024 :
        l="1: Read is PCR or optical duplicate"
    else:
        l="0: Read is not PCR or optical duplicate"
    flags.append(l)

    if x & 2048 :
        l="1: Supplementary alignment"
    else:
        l="0: Not supplementary alignment"
    flags.append(l)

    return flags
570216
def ascii2h5(dirname, output_fname):
    

    import tarfile
    import sys
    from glob import glob
    from contextlib import closing

    # The datatype that will be used to store extinction, A0
    A0_dtype = 'float16'

    def load_samp_file(f, fname):
        # Parse filename
        fname_chunks = os.path.split(fname)[1].split('_')
        l = float(fname_chunks[0])
        b = float(fname_chunks[1])

        # Load ASCII data
        data_raw = np.loadtxt(f, dtype='float64')

        n_samples = data_raw.shape[1] - 1
        n_dists = data_raw.shape[0]

        # Construct output
        dtype = [
            ('dist', 'int32'),
            ('A0', A0_dtype, (n_samples,))]

        data = np.empty(n_dists, dtype=dtype)
        data['dist'][:] = data_raw[:,0]
        data['A0'][:,:] = data_raw[:,1:]

        return (l,b), data

    def process_tarball(tarball_fname):
        # Write to the progress bar
        print('.', end='')
        sys.stdout.flush()

        with closing(tarfile.open(tarball_fname, mode='r:gz')) as f_tar:
            fnames = f_tar.getnames()

            f = f_tar.extractfile(fnames[0])
            (l,b), data = load_samp_file(f, fnames[0])
            n_dists, n_samples = data['A0'].shape
            n_coords = len(fnames)

            dtype = [
                ('l', 'float32'),
                ('b', 'float32'),
                ('dist', 'int32', (n_dists,)),
                ('A0', A0_dtype, (n_dists, n_samples))]
            data_combined = np.empty(n_coords, dtype=dtype)

            for k,fn in enumerate(fnames):
                # print('File {: >4d} of {:d}'.format(k+1, n_coords))
                f = f_tar.extractfile(fn)
                (l,b), data = load_samp_file(f, fn)
                data_combined['l'][k] = l
                data_combined['b'][k] = b
                data_combined['dist'][k] = data['dist']
                data_combined['A0'][k] = data['A0']

        return data_combined

    def save_data(data, fname):
        with closing(h5py.File(fname, 'w')) as f:
            f.create_dataset(
                'samples',
                data=data,
                chunks=True,
                compression='gzip',
                compression_opts=3)

    print('Progress: ', end='')
    sys.stdout.flush()

    tar_fname_list = glob(os.path.join(dirname, 'A_samp_*.tar.gz'))
    d = np.hstack([process_tarball(fn) for fn in tar_fname_list])

    print('+', end='')
    sys.stdout.flush()

    save_data(d, output_fname)

    print('')
580998
def bam_to_fastq_pair(in_file, target_region, pair):
    
    space, start, end = target_region
    bam_file = pysam.Samfile(in_file, "rb")
    for read in bam_file:
        if (not read.is_unmapped and not read.mate_is_unmapped
                and bam_file.getrname(read.tid) == space
                and bam_file.getrname(read.mrnm) == space
                and read.pos >= start and read.pos <= end
                and read.mpos >= start and read.mpos <= end
                and not read.is_secondary
                and read.is_paired and getattr(read, "is_read%s" % pair)):
            seq = Seq.Seq(read.seq)
            qual = list(read.qual)
            if read.is_reverse:
                seq = seq.reverse_complement()
                qual.reverse()
            yield read.qname, str(seq), "".join(qual)
1208082
def convertReadAlignment(self, read, readGroupSet, readGroupId):
        
        samFile = self.getFileHandle(self._dataUrl)
        # TODO fill out remaining fields
        # TODO refine in tandem with code in converters module
        ret = protocol.ReadAlignment()
        # ret.fragmentId = 'TODO'
        ret.aligned_quality.extend(read.query_qualities)
        ret.aligned_sequence = read.query_sequence
        if SamFlags.isFlagSet(read.flag, SamFlags.READ_UNMAPPED):
            ret.ClearField("alignment")
        else:
            ret.alignment.CopyFrom(protocol.LinearAlignment())
            ret.alignment.mapping_quality = read.mapping_quality
            ret.alignment.position.CopyFrom(protocol.Position())
            ret.alignment.position.reference_name = samFile.getrname(
                read.reference_id)
            ret.alignment.position.position = read.reference_start
            ret.alignment.position.strand = protocol.POS_STRAND
            if SamFlags.isFlagSet(read.flag, SamFlags.READ_REVERSE_STRAND):
                ret.alignment.position.strand = protocol.NEG_STRAND
            for operation, length in read.cigar:
                gaCigarUnit = ret.alignment.cigar.add()
                gaCigarUnit.operation = SamCigar.int2ga(operation)
                gaCigarUnit.operation_length = length
                gaCigarUnit.reference_sequence = ""  # TODO fix this!
        ret.duplicate_fragment = SamFlags.isFlagSet(
            read.flag, SamFlags.DUPLICATE_READ)
        ret.failed_vendor_quality_checks = SamFlags.isFlagSet(
            read.flag, SamFlags.FAILED_QUALITY_CHECK)
        ret.fragment_length = read.template_length
        ret.fragment_name = read.query_name
        for key, value in read.tags:
            # Useful for inspecting the structure of read tags
            # print("{key} {ktype}: {value}, {vtype}".format(
            #     key=key, ktype=type(key), value=value, vtype=type(value)))
            protocol.setAttribute(ret.attributes.attr[key].values, value)

        if SamFlags.isFlagSet(read.flag, SamFlags.MATE_UNMAPPED):
            ret.next_mate_position.Clear()
        else:
            ret.next_mate_position.Clear()
            if read.next_reference_id != -1:
                ret.next_mate_position.reference_name = samFile.getrname(
                    read.next_reference_id)
            else:
                ret.next_mate_position.reference_name = ""
            ret.next_mate_position.position = read.next_reference_start
            ret.next_mate_position.strand = protocol.POS_STRAND
            if SamFlags.isFlagSet(read.flag, SamFlags.MATE_REVERSE_STRAND):
                ret.next_mate_position.strand = protocol.NEG_STRAND
        if SamFlags.isFlagSet(read.flag, SamFlags.READ_PAIRED):
            ret.number_reads = 2
        else:
            ret.number_reads = 1
        ret.read_number = -1
        if SamFlags.isFlagSet(read.flag, SamFlags.FIRST_IN_PAIR):
            if SamFlags.isFlagSet(read.flag, SamFlags.SECOND_IN_PAIR):
                ret.read_number = 2
            else:
                ret.read_number = 0
        elif SamFlags.isFlagSet(read.flag, SamFlags.SECOND_IN_PAIR):
            ret.read_number = 1
        ret.improper_placement = not SamFlags.isFlagSet(
            read.flag, SamFlags.READ_PROPER_PAIR)
        ret.read_group_id = readGroupId
        ret.secondary_alignment = SamFlags.isFlagSet(
            read.flag, SamFlags.SECONDARY_ALIGNMENT)
        ret.supplementary_alignment = SamFlags.isFlagSet(
            read.flag, SamFlags.SUPPLEMENTARY_ALIGNMENT)
        ret.id = readGroupSet.getReadAlignmentId(ret)
        return ret
1301748
def stats(args):
    
    logger.info("Reading sequeces")
    data = parse_ma_file(args.ma)
    logger.info("Get sequences from sam")
    is_align = _read_sam(args.sam)
    is_json, is_db = _read_json(args.json)
    res = _summarise_sam(data, is_align, is_json, is_db)
    _write_suma(res, os.path.join(args.out, "stats_align.dat"))
    logger.info("Done")
1308240
def extract_sequences(self, ids, new_path=None):
        """"""
        # Temporary path #
        if new_path is None: new_fasta = self.__class__(new_temp_path())
        elif isinstance(new_path, FASTA): new_fasta = new_path
        else:                new_fasta = self.__class__(new_path)
        # Do it #
        new_fasta.create()
        for seq in self:
            if seq.id in ids: new_fasta.add_seq(seq)
        new_fasta.close()
        return new_fasta
1062248
def bam_to_fastq(bam_file, is_paired):
    
    out_files, out_handles = _get_fastq_handles(bam_file,
            is_paired)
    if len(out_handles) > 0:
        in_bam = pysam.Samfile(bam_file, mode='rb')
        for read in in_bam:
            num = 1 if (not read.is_paired or read.is_read1) else 2
            # reverse the sequence and quality if mapped to opposite strand
            if read.is_reverse:
                seq = str(Seq.reverse_complement(Seq.Seq(read.seq)))
                qual = "".join(reversed(read.qual))
            else:
                seq = read.seq
                qual = read.qual
            out_handles[num].write("@%s\n%s\n+\n%s\n" % (read.qname,
                seq, qual))
    [h.close() for h in out_handles.values()]
    return out_files
1013855
def get_sequences(options, fasta_file_contents):
    ''' 
    '''
    errors = []
    fasta_files_str = ", ".join(fasta_file_contents.keys())
    fasta_records = None
    reverse_mapping = {}

    try:
        fasta_records, reverse_mapping = parse_FASTA_files(options, fasta_file_contents)
        if not fasta_records:
            errors.append("No protein sequences found in the FASTA file(s) %s." % fasta_files_str)
    except Exception, e:
        e = '\n'.join([l for l in traceback.format_exc(), str('e') if l.strip()])
        errors.append("Error parsing FASTA file(s) %s:\n%s" % (fasta_files_str, str(e)))

    if not fasta_records:
        return None, {}, errors

    colorprinter.message('Found %d protein sequence(s).' % len(fasta_records))
    return fasta_records, reverse_mapping, errors
749136
def _parse_seqs(self, LOS):
        
        self.nseqs = len(LOS)
        self.width = len(LOS[0])
        for i in range(self.width):
            Dc = {'A': 0, 'C': 0, 'T': 0, 'G': 0, 'N': 0}
            for seq in LOS:
                key = seq[i]
                Dc[key] = Dc[key] + 1
            del(Dc['N'])
            self.counts.append(Dc)
1185087
def get_reads(sam, \
        contigs = False, mismatches = False, mm_option = False, \
        sort_sam = True, req_map = False, region = False, sbuffer = False):
    
    tempdir = '%s/' % (os.path.abspath(sam).rsplit('/', 1)[0])
    if sort_sam is True:
        mapping = '%s.sorted.sam' % (sam.rsplit('.', 1)[0])
        if sam != '-':
            if os.path.exists(mapping) is False:
                os.system("\
                    sort -k1 --buffer-size=%sG -T %s -o %s %s\
                    " % (sbuffer, tempdir, mapping, sam))
        else:
            mapping = 'stdin-sam.sorted.sam'
            p = Popen("sort -k1 --buffer-size=%sG -T %s -o %s" \
                    % (sbuffer, tempdir, mapping), stdin = sys.stdin, shell = True)
            p.communicate()
        mapping = open(mapping)
    else:
        if sam == '-':
            mapping = sys.stdin
        else:
            mapping = open(sam)
    for read in reads_from_mapping(mapping, contigs, mismatches, mm_option, req_map, region):
        yield read
1181625
def get_total_alignment_score(bam):
    ''''''
    sam_reader = pysam.Samfile(bam, "rb")
    total = 0
    for sam in sam_reader.fetch(until_eof=True):
        try:
            total += sam.opt('AS')
        except:
            pass
    return total
1162886
def select_unaligned_read_pairs(in_bam, extra, out_dir, config):
    
    runner = broad.runner_from_config(config)
    base, ext = os.path.splitext(os.path.basename(in_bam))
    nomap_bam = os.path.join(out_dir, "{}-{}{}".format(base, extra, ext))
    if not utils.file_exists(nomap_bam):
        with file_transaction(nomap_bam) as tx_out:
            runner.run("FilterSamReads", [("INPUT", in_bam),
                                          ("OUTPUT", tx_out),
                                          ("EXCLUDE_ALIGNED", "true"),
                                          ("WRITE_READS_FILES", "false"),
                                          ("SORT_ORDER", "queryname")])
    has_reads = False
    with pysam.Samfile(nomap_bam, "rb") as in_pysam:
        for read in in_pysam:
            if read.is_paired:
                has_reads = True
                break
    if has_reads:
        out_fq1, out_fq2 = ["{}-{}.fq".format(os.path.splitext(nomap_bam)[0], i) for i in [1, 2]]
        runner.run_fn("picard_bam_to_fastq", nomap_bam, out_fq1, out_fq2)
        return out_fq1, out_fq2
    else:
        return None, None
73234
def _par_read(dirname, compressed=True):
    
    templates = []
    if compressed:
        arc = tarfile.open(dirname, "r:*")
        members = arc.getmembers()
        _parfile = [member for member in members
                    if member.name.split(os.sep)[-1] ==
                    'template_parameters.csv']
        if len(_parfile) == 0:
            arc.close()
            raise MatchFilterError(
                'No template parameter file in archive')
        parfile = arc.extractfile(_parfile[0])
    else:
        parfile = open(dirname + '/' + 'template_parameters.csv', 'r')
    for line in parfile:
        t_in = Template()
        for key_pair in line.rstrip().split(','):
            if key_pair.split(':')[0].strip() == 'name':
                t_in.__dict__[key_pair.split(':')[0].strip()] = \
                    key_pair.split(':')[-1].strip()
            elif key_pair.split(':')[0].strip() == 'filt_order':
                try:
                    t_in.__dict__[key_pair.split(':')[0].strip()] = \
                        int(key_pair.split(':')[-1])
                except ValueError:
                    pass
            else:
                try:
                    t_in.__dict__[key_pair.split(':')[0].strip()] = \
                        float(key_pair.split(':')[-1])
                except ValueError:
                    pass
        templates.append(t_in)
    parfile.close()
    if compressed:
        arc.close()
    return templates
431769
def parse_tmhmm_long(tmhmm_results):
    
    with open(tmhmm_results) as f:
        lines = f.read().splitlines()

    infodict = defaultdict(dict)

    for l in lines:
        if 'Number of predicted TMHs:' in l:
            gene = l.split(' Number')[0].strip('# ')
            infodict[gene]['num_tm_helices'] = int(l.split(': ')[1])

        if 'WARNING' in l:
            log.warning('{}: no TMHMM predictions'.format(l))
            continue

        # TODO: POSSIBLE N-term signal sequence - parse this

        # Look for the lines without #, these are the TM predicted regions
        if '#' not in l:
            stuff = l.split()
            if stuff[1] == 'TMHMM2.0':

                gene = stuff[0]
                region = stuff[2]
                region_start = stuff[3]
                region_end = stuff[4]

                if 'sequence' in infodict[gene]:
                    tm_seq = infodict[gene]['sequence']
                else:
                    tm_seq = ''

                if region == 'outside':
                    info = 'O'
                elif region == 'inside':
                    info = 'I'
                elif region == 'TMhelix':
                    info = 'T'
                else:
                    log.error('{}: unknown region type'.format(info))
                    info = '-'

                for r in range(int(region_start), int(region_end) + 1):
                    tm_seq += info

                infodict[gene]['sequence'] = tm_seq

    return infodict