### Navigation

- index
- modules |
- SamSifter 1.0.0 documentation »
- Module code »

# Source code for samsifter.util.papertrail

    #!/usr/bin/env python3
    # -*- coding: utf-8 -*-
    """Annotation of SAM file headers with metadata on applied programs.

    .. moduleauthor:: Florian Aldehoff <samsifter@biohazardous.de>
    """
    import sys
    import fileinput

    def find_last_in_chain(current_id, pg_entries):
        """Recursively identify the last program entry in the chain.

        @PG entries in SAM file headers can be chained by linking each entry to the
        previously applied program using the optional PP tag. The actual last entry
        of the chain can only be found under the assumption that all programs are
        listed within one chain.

        Note
        ----
        May lead to wrong assumptions about the actual order of programs applied to
        dataset if the chain of PG entries is interrupted by missing PP tags. Use
        at own risk!

        Parameters
        ----------
        current_id : str
            Unique ID of program entry that serves as arbitrary start point for the
            recursive search of the chain end.
        pg_entries : dict of dicts
            Dictionary of @PG record tags (key:value) referenced by their unique
            program ID.

        Returns
        -------
        str
            Unique ID of the assumed last entry in the program chain (see note
            above!)
        """
        for entry, tags in pg_entries.items():
            try:
                if current_id == tags['PP']:
                    # another entry follows
                    return find_last_in_chain(entry, pg_entries)
            except KeyError:
                # entry has no PP tag
                continue
        # no other entry follows
        return current_id

    def modify_header(lines, name=None, command_line=None, description=None,
                      version=None):
        """Modify SAM header by inserting a new @PG record into the chain.

        Ensures unique program identifiers and integrity of program chain.
        New entries are always inserted after all current PG entries and their
        order is not changed in case the entries are not chained by optional PP
        tags.

        Note
        ----
        Chaining of entries with PP tags is disabled as the chain of programs is
        often incomplete and may be misleading! The recursion through existing
        entries may not always be sufficient to identify the current last entry in
        the chain.

        Parameters
        ----------
        lines : list of str
            Full set of header lines of SAM file, including the first header line
            with @HD record.
        name : str, optional
            Program name, defaults to None.
        command_line : str, optional
            Program command line excluding program name/call, just the optional and
            positional arguments. Defaults to None.
        description : str, optional
            Short description of program functionality.
        version : str, optional
            Program version number/string.
        """
        first_header = None       # first header line, required if header present
        pg_list = []              # list of PG header lines
        pg_entries = {}           # dict of dict for program records and their tags
        others = []               # list of other header lines

        # parse old header
        for line in lines:
            if line.startswith('@HD'):
                first_header = line
            elif line.startswith('@PG'):
                pg_list.append(line)
                fields = line.lstrip('@PG\t').split('\t')
                tags = {}
                identifier = None
                for field in fields:
                    (key, value) = field.rstrip().split(':')
                    if key == 'ID':
                        identifier = value
                    tags[key] = value

                # store program entry
                pg_entries[identifier] = tags
            else:
                others.append(line)

        # disabled as chain is often incomplete and may be misleading!
    #    # recursively identify last program in chain, start with arbitrary entry
    #    current_id = next(iter(pg_entries.keys()))
    #    last = find_last_in_chain(current_id, pg_entries)

        # generate unique but short and consecutive identifier for new @PG record
        new_id = str(name)
        suffix = 0
        while new_id in pg_entries.keys():
            suffix += 1
            new_id = "%s_%i" % (name, suffix)

        # create new @PG line
        new_pg = create_pg_line(new_id,
                                name=name,
                                command_line=command_line,
                                # deactivated chaining, see above
                                # previous=last,
                                previous=None,
                                description=description,
                                version=version)

        # insert new PG line into the rest of the lines
        modified = [first_header]
        modified.extend(pg_list)
        modified.append(new_pg)
        modified.extend(others)

        return modified

    def create_pg_line(identifier, name=None, command_line=None, previous=None,
                       description=None, version=None):
        """Create a @PG header for SAM files.

        Follows SAMv1 specification, see
        https://samtools.github.io/hts-specs/SAMv1.pdf

        Parameters
        ----------
        identifier : str
            Unique program identifier within the context of the current SAM file
            that is identical to the identifier used in the alignment section. May
            be modified by later processing steps.
        name : str, optional
            Program name, defaults to None.
        command_line : str, optional
            Program command line excluding program name/call, just the optional and
            positional arguments. Defaults to None.
        previous : str, optional
            Program identifier of previous step in program chain.
        description : str, optional
            Short description of program functionality.
        version : str, optional
            Program version number/string.

        Returns
        -------
        str
            Complete standard-compliant @PG header line.
        """

        line = '@PG'

        # Program record identifier. Each @PG line must have a unique ID. The value
        # of ID is used in the alignment PG tag and PP tags of other @PG lines. PG
        # IDs may be modified when merging SAM files in order to handle collisions.
        pg_id = 'ID:%s' % identifier
        line += '\t' + pg_id

        # Program name, optional
        if name is not None:
            pg_pn = 'PN:%s' % name
            line += '\t' + pg_pn

        # Command line, optional
        if command_line is not None:
            pg_cl = 'CL:%s' % command_line
            line += '\t' + pg_cl

        # Previous @PG-ID. Must match another @PG header’s ID tag. @PG records may
        # be chained using PP tag, with the last record in the chain having no PP
        # tag. This chain defines the order of programs that have been applied to
        # the alignment. PP values may be modified when merging SAM files in order
        # to handle collisions of PG IDs. The first PG record in a chain (i.e. the
        # one referred to by the PG tag in a SAM record) describes the most recent
        # program that operated on the SAM record. The next PG record in the chain
        # describes the next most recent program that operated on the SAM record.
        # The PG ID on a SAM record is not required to refer to the newest PG
        # record in a chain. It may refer to any PG record in a chain, implying
        # that the SAM record has been operated on by the program in that PG
        # record, and the program(s) referred to via the PP tag.
        if previous is not None:
            pg_pp = 'PP:%s' % previous
            line += '\t' + pg_pp

        # Description, optional
        if description is not None:
            pg_ds = 'DS:%s' % description
            line += '\t' + pg_ds

        # Program version, optional
        if version is not None:
            pg_vn = 'VN:%s' % version
            line += '\t' + pg_vn

        line += '\n'

        return line

    def main():
        """Simple test of PG header generation."""
        handle = fileinput.input(sys.argv)
        header = []
        for line in handle:
            if line.startswith('@'):
                header.append(line)
        handle.close()

        modified = modify_header(header, name='Test')
        for line in modified:
            print(line.rstrip(), file=sys.stdout)

    if __name__ == "__main__":
        main()

### Navigation

- index
- modules |
- SamSifter 1.0.0 documentation »
- Module code »