API Reference

This page dynamically pulls docstrings from the GuideMaker source files.

Core Modules

Core classes and functions for GuideMaker.

Annotation

Annotation class for data and methods on targets and gene annotations.

Source code in guidemaker/core.py

class Annotation:

    """
    Annotation class for data and methods on targets and gene annotations.

    """

    def __init__(self, annotation_list: List[str], annotation_type: str, target_bed_df: object) -> None:
        """
        Annotation class for data and methods on targets and gene annotations

        Args:
            annotation_list (List[str]): A list of genbank files from a single genome
            annotation_type (str): "genbank" | "gff"
            target_bed_df (object): A pandas dataframe in Bed format with the
                locations of targets in the genome

        Returns:
            None
        """
        self.annotation_list: List[str] = annotation_list
        self.annotation_type = annotation_type
        self.target_bed_df: object = target_bed_df
        self.genbank_bed_df: object = None
        self.feature_dict: Dict = None
        self.nearby: object = None
        self.filtered_df: object = None
        self.qualifiers: object = None

    def check_annotation_type(self):
        """open GTF/GFF and determine if the file provided by the GFF argument is a GFF or GTF file

            Args: None

            Returns (str): ["gff" | "gtf"]
        """
        def search(f):
            line1 = f.readline()
            gffmatch = re.search("gff-version", line1)
            if gffmatch is not None:
                return "gff"
            gtfmatch = re.search("gtf-version", line1)
            if gtfmatch is not None:
                return "gtf"
            else:
                logger.error("Could not verify the GFF/GTF file type. Please make sure your GFF/GTF file starts with '#gtf-version' or '##gff-version'")
                raise ValueError
        testfile = self.annotation_list[0]
        if is_gzip(testfile):
            with gzip.open(testfile, 'rt') as f:
                return search(f)
        else:
            with open(testfile, 'r') as f:
                return search(f)

    def get_annotation_features(self, feature_types: List[str] = None) -> None:
        """
        Parse annotation records into pandas DF/Bed format and dict format saving to self

        Args:
            feature_types (List[str]): a list of Genbank feature types to use

        Returns:
            None
        """
        if feature_types is None:
            feature_types = ["CDS"]
        feature_dict = {}
        pddict = dict(chrom=[], chromStart=[], chromEnd=[], name=[], strand=[])
        if self.annotation_type == "genbank":
            for gbfile in self.annotation_list:
                try:
                    if is_gzip(gbfile):
                        f = gzip.open(gbfile, mode='rt')
                    else:
                        f = open(gbfile, mode='r')
                except IOError as e:
                    logger.error("The genbank file %s could not be opened" % gbfile)
                    raise e
                genbank_file = SeqIO.parse(f, "genbank")
                for entry in genbank_file:
                    for record in entry.features:
                        if record.type in feature_types:
                            strand_val = getattr(record, 'strand', getattr(record.location, 'strand', None))
                            if strand_val in [1, -1, "+", "-"]:
                                pddict["strand"].append("-" if str(strand_val) in ['-1', '-'] else "+")
                            featid = hashlib.md5(str(record).encode()).hexdigest()
                            pddict['chrom'].append(entry.id)
                            pddict["chromStart"].append(int(record.location.start))
                            pddict["chromEnd"].append(int(record.location.end))
                            pddict["name"].append(featid)
                            for qualifier_key, qualifier_val in record.qualifiers.items():
                                if not qualifier_key in feature_dict:
                                    feature_dict[qualifier_key] = {}
                                feature_dict[qualifier_key][featid] = qualifier_val
            genbankbed = pd.DataFrame.from_dict(pddict)
            self.genbank_bed_df = genbankbed
            self.feature_dict = feature_dict
            f.close()
        elif self.annotation_type == "gff":
            anno_format = self.check_annotation_type()
            for gff in self.annotation_list:
                gr = pr.read_gff3(gff) if anno_format == 'gff' else pr.read_gtf(gff)
                df = gr.df
                filtered_df = df[df['Feature'].isin(feature_types)]
                for idx, row in filtered_df.iterrows():
                    pddict['chrom'].append(str(row['Chromosome']))
                    pddict['chromStart'].append(int(row['Start']))
                    pddict['chromEnd'].append(int(row['End']))
                    pddict['strand'].append(str(row['Strand']))
                    rec_str = f"{row['Chromosome']}_{row['Start']}_{row['End']}_{row['Feature']}_{idx}"
                    featid = hashlib.md5(rec_str.encode()).hexdigest()
                    pddict['name'].append(featid)

                    for col in df.columns:
                        if col not in ['Chromosome', 'Source', 'Feature', 'Start', 'End', 'Score', 'Strand', 'Frame', 'Attribute']:
                            val = row[col]
                            if pd.notna(val) and str(val).strip() != '':
                                if col not in feature_dict:
                                    feature_dict[col] = {}
                                feature_dict[col][featid] = str(val)
            genbankbed = pd.DataFrame.from_dict(pddict)
            self.genbank_bed_df = genbankbed
            self.feature_dict = feature_dict


    def _get_qualifiers(self, configpath, excluded: List[str] = None) -> object:
        """
        Create a dataframe with features and their qualifier values

        Create a dataframe with features and their qualifier values for
        all qualifiers over the minimum threshold (except 'translation'). Add
        to self.qualifiers

        Args:
            min_prop (float): A float between 0-1 representing the fraction of
            features the qualifier must be present in to be included in the dataframe
            excluded (List(str)): A list of genbank qualifiers to exclude, Default ["translation"]

        Returns:
            None
        """
        with open(configpath) as cf:
            config = yaml.safe_load(cf)

        min_prop = config['MINIMUM_PROPORTION']

        if excluded is None:
            excluded = ["translation"]
        final_quals = []
        qual_df = pd.DataFrame(data={"Feature id": []})
        for featkey, quals in self.feature_dict.items():
            if len(quals) / len(self.feature_dict[featkey]) > min_prop:
                final_quals.append(featkey)
        for qualifier in final_quals:
            if qualifier not in excluded:
                featlist = []
                quallist = []
                for feat, qual in self.feature_dict[qualifier].items():
                    featlist.append(feat)
                    if isinstance(qual, list):
                        quallist.append(";".join([str(i) for i in qual]))
                    else:
                        quallist.append(qual)
                tempdf = pd.DataFrame({'Feature id': featlist, qualifier: quallist})
                qual_df = qual_df.merge(tempdf, how="outer", on="Feature id")
        self.qualifiers = qual_df

    def _get_nearby_features(self) -> None:
        """
        Adds downstream and upstream feature information to target sequences using NumPy regional window search.

        Args:
            None

        Returns:
            None

        Note:
            Writes a dataframe of nearby features to self.nearby
        """
        target_sorted = self.target_bed_df.sort_values(by=['chrom', 'chromstart']).reset_index(drop=True)
        feature_sorted = self.genbank_bed_df.sort_values(by=['chrom', 'chromStart']).reset_index(drop=True)

        rows_down = []
        rows_up = []

        for chrom, t_sub in target_sorted.groupby('chrom', observed=True):
            f_sub = feature_sorted[feature_sorted['chrom'] == chrom]
            if f_sub.empty:
                continue

            f_starts = f_sub['chromStart'].to_numpy(dtype=np.int64)
            f_ends = f_sub['chromEnd'].to_numpy(dtype=np.int64)
            f_strands = f_sub['strand'].to_numpy()
            f_ids = f_sub['name'].to_numpy()
            max_feat_len = int(np.max(f_ends - f_starts)) if len(f_starts) > 0 else 0

            t_starts = t_sub['chromstart'].to_numpy(dtype=np.int64)
            t_ends = t_sub['chromend'].to_numpy(dtype=np.int64)
            t_seqs = t_sub['name'].to_numpy()
            t_strands = t_sub['strand'].to_numpy()

            n_targets = len(t_starts)
            batch_size = 5000

            for b_i in range(0, n_targets, batch_size):
                b_s = t_starts[b_i:b_i + batch_size]
                b_e = t_ends[b_i:b_i + batch_size]
                b_seq = t_seqs[b_i:b_i + batch_size]
                b_str = t_strands[b_i:b_i + batch_size]
                B = len(b_s)

                w_start = max(0, b_s[0] - 100000)
                w_end = b_e[-1] + 100000

                f_idx_start = np.searchsorted(f_starts, max(0, w_start - max_feat_len), side='left')
                f_idx_end = np.searchsorted(f_starts, w_end, side='right')

                s_f_starts = f_starts[f_idx_start:f_idx_end]
                s_f_ends = f_ends[f_idx_start:f_idx_end]
                s_f_strands = f_strands[f_idx_start:f_idx_end]
                s_f_ids = f_ids[f_idx_start:f_idx_end]

                if len(s_f_starts) == 0:
                    s_f_starts = f_starts
                    s_f_ends = f_ends
                    s_f_strands = f_strands
                    s_f_ids = f_ids

                K = len(s_f_starts)

                fs_2d = s_f_starts[None, :]  # 1 x K
                fe_2d = s_f_ends[None, :]    # 1 x K
                bs_2d = b_s[:, None]         # B x 1
                be_2d = b_e[:, None]         # B x 1

                # Downstream features: f_end > g_start
                down_valid = fe_2d > bs_2d
                overlaps_down = (fs_2d < be_2d) & (fe_2d > bs_2d)
                dists_down = np.where(overlaps_down, 0, np.where(fs_2d >= be_2d, fs_2d - be_2d, bs_2d - fe_2d))
                dists_down = np.where(down_valid, dists_down, np.inf)

                best_down_idx = np.argmin(dists_down, axis=1)
                has_down = np.min(dists_down, axis=1) < np.inf

                for i in range(B):
                    if has_down[i]:
                        idx = best_down_idx[i]
                        rows_down.append({
                            'Accession': chrom, 'Guide start': b_s[i], 'Guide end': b_e[i],
                            'Guide sequence': b_seq[i], 'Guide strand': b_str[i],
                            'Feature Accession': chrom, 'Feature start': s_f_starts[idx],
                            'Feature end': s_f_ends[idx], 'Feature id': s_f_ids[idx],
                            'Feature strand': s_f_strands[idx], 'Feature distance': int(dists_down[i, idx]),
                            'direction': 'downstream'
                        })

                # Upstream features: f_start < g_end
                up_valid = fs_2d < be_2d
                overlaps_up = (fs_2d < be_2d) & (fe_2d > bs_2d)
                dists_up = np.where(overlaps_up, 0, np.where(fs_2d >= be_2d, fs_2d - be_2d, bs_2d - fe_2d))
                dists_up = np.where(up_valid, dists_up, np.inf)

                best_up_idx = np.argmin(dists_up, axis=1)
                has_up = np.min(dists_up, axis=1) < np.inf

                for i in range(B):
                    if has_up[i]:
                        idx = best_up_idx[i]
                        rows_up.append({
                            'Accession': chrom, 'Guide start': b_s[i], 'Guide end': b_e[i],
                            'Guide sequence': b_seq[i], 'Guide strand': b_str[i],
                            'Feature Accession': chrom, 'Feature start': s_f_starts[idx],
                            'Feature end': s_f_ends[idx], 'Feature id': s_f_ids[idx],
                            'Feature strand': s_f_strands[idx], 'Feature distance': int(dists_up[i, idx]),
                            'direction': 'upstream'
                        })

        df_down = pd.DataFrame(rows_down)
        df_up = pd.DataFrame(rows_up)
        nearby = pd.concat([df_down, df_up], axis=0, ignore_index=True)
        self.nearby = nearby


    def _filter_features(self, before_feat: int = 100, after_feat: int = 200 ) -> None:
        """
        Merge targets with Feature list and filter for guides close enough to interact.

        Args:
            before_feat (int): The maximum distance before the start of a feature measured from closest point to guide
            after_feat (int): The maximum distance after the start codon (into the gene)

        Returns:
            None
        """
        # for guides in the same orientation as the targets ( +/+ or -/-) select guides that are within
        #  before_feat of the gene start
        filtered_df = self.nearby.query(
            '`Guide strand` == `Feature strand` and 0 < `Feature distance` < @before_feat')
        # for guides in the +/+ orientation select guides where the end is within [before_feat] of the gene start
        p1 = (self.nearby.query('`Guide strand` == "+" and `Feature strand` == "+" \
                                             and `Feature distance` == 0 and \
                                             `Guide end` - `Feature start` < @after_feat'))
        # for guides in the -/- orientation select guides where the end is within [before_feat] of the gene start
        p2 = (self.nearby.query('`Guide strand` == "-" and `Feature strand` == "-" \
                                                     and `Feature distance` == 0 \
                                                     and `Feature end` - `Guide start` < @after_feat'))
        # Select guides where target is + and guide is - and the guide is infront of the gene
        p3 = (self.nearby.query('`Guide strand` == "-" and `Feature strand` == "+" and \
                                                     0 <`Feature start` - `Guide end` < @before_feat'))
        # Select guides where target is - and guide is + and the guide is infront of the gene
        p4 = (self.nearby.query('`Guide strand` == "+" and `Feature strand` == "-" and \
                                                     0 <`Guide start` - `Feature end` < @before_feat'))
        # Select guides where target is + and guide is - and the guide is is within [before_feat] of the gene start
        p5 = (self.nearby.query('`Guide strand` == "-" and `Feature strand` == "+" and \
                                                             0 <`Guide end` -`Feature start`  < @after_feat'))
        # Select guides where target is - and guide is + and the guide is is within [before_feat] of the gene start
        p6 = (self.nearby.query('`Guide strand` == "+" and `Feature strand` == "-" and \
                                                             0 <`Feature end` - `Guide start` < @after_feat'))
        self.filtered_df = pd.concat([filtered_df, p1, p2, p3, p4, p5, p6], axis=0)

    def _format_guide_table(self, targetprocessor_object) -> pd.DataFrame:
        """
        Create guide table for output

        Args:
            target- a dataframe with targets from targetclass

        Returns:
            (PandasDataFrame): A formated pandas dataframe
        """
        def gc(seq):
            cnt = 0
            for letter in seq:
                if letter in ["G", "C"]:
                    cnt += 1
            return cnt / len(seq)

        def get_guide_hash(seq):
            return hashlib.md5(seq.encode()).hexdigest()

        def checklen30(seq):
            if len(seq) == 30:
                return True
            return False

        def get_off_target_score(seq):
            dlist = targetprocessor_object.neighbors[seq]["neighbors"]["dist"]
            s = [str(i) for i in dlist]
            return ";".join(s)

        def get_off_target_seqs(seq):
            slist = targetprocessor_object.neighbors[seq]["neighbors"]["seqs"]
            return ";".join(slist)
        pretty_df = deepcopy(self.filtered_df)  # anno class object
        # retrive only guides that are in neighbors keys.
        pretty_df = pretty_df[pretty_df["Guide sequence"].isin(
            list(targetprocessor_object.neighbors.keys()))]
        pretty_df['GC'] = pretty_df['Guide sequence'].apply(gc)
        pretty_df['Guide name'] = pretty_df['Guide sequence'].apply(get_guide_hash)
        pretty_df['Target strand'] = np.where(
            pretty_df['Guide strand'] == pretty_df['Feature strand'], 'coding', 'non-coding')
        pretty_df['Similar guide distances'] = pretty_df['Guide sequence'].apply(
            get_off_target_score)
        pretty_df['Similar guides'] = pretty_df['Guide sequence'].apply(get_off_target_seqs)
        pretty_df = pd.merge(pretty_df, targetprocessor_object.targets, how="left",
         left_on=['Guide sequence', 'Guide start', 'Guide end', 'Accession'],
            right_on=['target', 'start', 'stop', 'seqid'])
        # rename exact_pam to PAM
        pretty_df = pretty_df.rename(columns={"exact_pam": "PAM"})

        pretty_df = pretty_df[['Guide name', 'Guide sequence', 'GC', 'dtype', 'Accession', 'Guide start', 'Guide end',
                    'Guide strand', 'PAM', 'Feature id',
                    'Feature start', 'Feature end', 'Feature strand',
                    'Feature distance', 'Similar guides', 'Similar guide distances','target_seq30']]
        pretty_df = pretty_df.merge(self.qualifiers, how="left", on="Feature id")
        pretty_df = pretty_df.sort_values(by=['Accession', 'Feature start'])
        # to match with the numbering with other tools- offset
        pretty_df['Guide start'] = pretty_df['Guide start'] + 1
        pretty_df['Feature start'] = pretty_df['Feature start'] + 1
        pretty_df=pretty_df.loc[pretty_df['target_seq30'].apply(checklen30)==True].reset_index(drop=True)
        self.pretty_df = pretty_df

    def _filterlocus(self, attribute:str , filter_by_locus:list = []) -> pd.DataFrame:
        """
        Create guide table for output for a selected attribute type

        Args:
            attribute: The key in the attributes column (column 9) of the GFF/GTF file to filter on
            filter_by_locus: A list of Identifiers to filter the full data frame by

        Returns:
            (PandasDataFrame): A formated pandas dataframe
        """

        df = deepcopy(self.pretty_df)  # anno class object
        if len (filter_by_locus) > 0:
            df = df[df[attribute].isin(filter_by_locus)]
        return df

    def locuslen(self) -> int:
        """
        Count the number of locus tag in the genebank file

        Args:
            None

        Returns:
            (int): Number of locus tag
        """
        da_keys = self.feature_dict.keys()
        firsttag = (list(da_keys)[0])
        if firsttag:
            locus_count = len(self.feature_dict[firsttag].keys())
            return firsttag, locus_count
        else:
            logger.warning("A locus key could not be found.")
            return "notag", 0

__init__(annotation_list, annotation_type, target_bed_df)

Annotation class for data and methods on targets and gene annotations

Parameters:

Name	Type	Description	Default
annotation_list	List[str]	A list of genbank files from a single genome	required
annotation_type	str	"genbank" | "gff"	required
target_bed_df	object	A pandas dataframe in Bed format with the locations of targets in the genome	required

Returns:

Type	Description
None	None

Source code in guidemaker/core.py

def __init__(self, annotation_list: List[str], annotation_type: str, target_bed_df: object) -> None:
    """
    Annotation class for data and methods on targets and gene annotations

    Args:
        annotation_list (List[str]): A list of genbank files from a single genome
        annotation_type (str): "genbank" | "gff"
        target_bed_df (object): A pandas dataframe in Bed format with the
            locations of targets in the genome

    Returns:
        None
    """
    self.annotation_list: List[str] = annotation_list
    self.annotation_type = annotation_type
    self.target_bed_df: object = target_bed_df
    self.genbank_bed_df: object = None
    self.feature_dict: Dict = None
    self.nearby: object = None
    self.filtered_df: object = None
    self.qualifiers: object = None

check_annotation_type()

open GTF/GFF and determine if the file provided by the GFF argument is a GFF or GTF file

Args: None

Returns (str): ["gff" | "gtf"]

Source code in guidemaker/core.py

def check_annotation_type(self):
    """open GTF/GFF and determine if the file provided by the GFF argument is a GFF or GTF file

        Args: None

        Returns (str): ["gff" | "gtf"]
    """
    def search(f):
        line1 = f.readline()
        gffmatch = re.search("gff-version", line1)
        if gffmatch is not None:
            return "gff"
        gtfmatch = re.search("gtf-version", line1)
        if gtfmatch is not None:
            return "gtf"
        else:
            logger.error("Could not verify the GFF/GTF file type. Please make sure your GFF/GTF file starts with '#gtf-version' or '##gff-version'")
            raise ValueError
    testfile = self.annotation_list[0]
    if is_gzip(testfile):
        with gzip.open(testfile, 'rt') as f:
            return search(f)
    else:
        with open(testfile, 'r') as f:
            return search(f)

get_annotation_features(feature_types=None)

Parse annotation records into pandas DF/Bed format and dict format saving to self

Parameters:

Name	Type	Description	Default
feature_types	List[str]	a list of Genbank feature types to use	None

Returns:

Type	Description
None	None

Source code in guidemaker/core.py

def get_annotation_features(self, feature_types: List[str] = None) -> None:
    """
    Parse annotation records into pandas DF/Bed format and dict format saving to self

    Args:
        feature_types (List[str]): a list of Genbank feature types to use

    Returns:
        None
    """
    if feature_types is None:
        feature_types = ["CDS"]
    feature_dict = {}
    pddict = dict(chrom=[], chromStart=[], chromEnd=[], name=[], strand=[])
    if self.annotation_type == "genbank":
        for gbfile in self.annotation_list:
            try:
                if is_gzip(gbfile):
                    f = gzip.open(gbfile, mode='rt')
                else:
                    f = open(gbfile, mode='r')
            except IOError as e:
                logger.error("The genbank file %s could not be opened" % gbfile)
                raise e
            genbank_file = SeqIO.parse(f, "genbank")
            for entry in genbank_file:
                for record in entry.features:
                    if record.type in feature_types:
                        strand_val = getattr(record, 'strand', getattr(record.location, 'strand', None))
                        if strand_val in [1, -1, "+", "-"]:
                            pddict["strand"].append("-" if str(strand_val) in ['-1', '-'] else "+")
                        featid = hashlib.md5(str(record).encode()).hexdigest()
                        pddict['chrom'].append(entry.id)
                        pddict["chromStart"].append(int(record.location.start))
                        pddict["chromEnd"].append(int(record.location.end))
                        pddict["name"].append(featid)
                        for qualifier_key, qualifier_val in record.qualifiers.items():
                            if not qualifier_key in feature_dict:
                                feature_dict[qualifier_key] = {}
                            feature_dict[qualifier_key][featid] = qualifier_val
        genbankbed = pd.DataFrame.from_dict(pddict)
        self.genbank_bed_df = genbankbed
        self.feature_dict = feature_dict
        f.close()
    elif self.annotation_type == "gff":
        anno_format = self.check_annotation_type()
        for gff in self.annotation_list:
            gr = pr.read_gff3(gff) if anno_format == 'gff' else pr.read_gtf(gff)
            df = gr.df
            filtered_df = df[df['Feature'].isin(feature_types)]
            for idx, row in filtered_df.iterrows():
                pddict['chrom'].append(str(row['Chromosome']))
                pddict['chromStart'].append(int(row['Start']))
                pddict['chromEnd'].append(int(row['End']))
                pddict['strand'].append(str(row['Strand']))
                rec_str = f"{row['Chromosome']}_{row['Start']}_{row['End']}_{row['Feature']}_{idx}"
                featid = hashlib.md5(rec_str.encode()).hexdigest()
                pddict['name'].append(featid)

                for col in df.columns:
                    if col not in ['Chromosome', 'Source', 'Feature', 'Start', 'End', 'Score', 'Strand', 'Frame', 'Attribute']:
                        val = row[col]
                        if pd.notna(val) and str(val).strip() != '':
                            if col not in feature_dict:
                                feature_dict[col] = {}
                            feature_dict[col][featid] = str(val)
        genbankbed = pd.DataFrame.from_dict(pddict)
        self.genbank_bed_df = genbankbed
        self.feature_dict = feature_dict

locuslen()

Count the number of locus tag in the genebank file

Returns:

Type	Description
int	Number of locus tag

Source code in guidemaker/core.py

def locuslen(self) -> int:
    """
    Count the number of locus tag in the genebank file

    Args:
        None

    Returns:
        (int): Number of locus tag
    """
    da_keys = self.feature_dict.keys()
    firsttag = (list(da_keys)[0])
    if firsttag:
        locus_count = len(self.feature_dict[firsttag].keys())
        return firsttag, locus_count
    else:
        logger.warning("A locus key could not be found.")
        return "notag", 0

GuideMakerPlot

A class with functions to plot guides over genome cooridinates.

Source code in guidemaker/core.py

class GuideMakerPlot:

    """
    A class with functions to plot guides over genome cooridinates.

    """

    def __init__(self, prettydf: pd.DataFrame, outdir: str) -> None:
        """
        GuideMakerPlot class for visualizing distrubution of gRNA, features, and locus.

        Args:
            prettydf (PandasDataFrame): Final output from GuideMaker
            outdir (str): Output Directory

        Returns:
            None
        """
        self.prettydf = prettydf
        self.accession = list(set(self.prettydf['Accession']))

        def _singleplot(df):
            """
            Returns guidemaker plot describing PAM targets

            Args:
                df(PandasDataFrame): Final output from GuideMaker for a single accession

            Return:
                None
            """
            source = df
            brush = alt.selection(type='interval', encodings=['x'])
            binNum = int(round(source['Feature end'].max() / 200, 0))
            display_info = source.columns.tolist()

            # Feature density
            densityF = alt.Chart(source).transform_density(
            'Feature start',
            as_=['Feature start', 'Feature Density'],
            extent=[1, source['Feature end'].max()],
            bandwidth=binNum,
            ).mark_area(color='black', opacity=0.6).encode(
            x=alt.X('Feature start', axis=alt.Axis(title='Genome Coordinates (bp)', tickCount=5)),
            y='Feature Density:Q',
            ).properties(height=50, width=500)

            # Guide density
            densityG = alt.Chart(source).transform_density(
            'Guide start',
            as_=['Guide start', 'Guide Density'],
            extent=[1, source['Feature end'].max()],
            bandwidth=binNum,
            ).mark_area(color='pink', opacity=0.6).encode(
            x=alt.X('Guide start', axis=alt.Axis(title='Genome Coordinates (bp)', tickCount=5)),
            y='Guide Density:Q',
            ).properties(height=50, width=500).add_selection(brush)

            # locus bar
            locus = alt.Chart(source).mark_bar(cornerRadiusTopLeft=3, cornerRadiusTopRight=3).encode(
            x='count(locus_tag):Q',
            y=alt.Y('locus_tag', axis=alt.Axis(title='Locus')),
            color='PAM:N',
            tooltip=display_info
            ).transform_filter(
            brush
            ).interactive().properties(height=500, width=500)
            guidemakerChart = (densityF & densityG & locus)
            return(guidemakerChart)

        outdir_real = os.path.realpath(outdir)
        for accession in self.accession:
            df = self.prettydf[self.prettydf['Accession'] == accession]
            accession_plot = _singleplot(df)
            clean_accession = re.sub(r'[^A-Za-z0-9._-]', '_', str(accession))
            plot_file_name = os.path.join(outdir, f"{clean_accession}.html")
            plot_file_real = os.path.realpath(plot_file_name)
            if not plot_file_real.startswith(outdir_real):
                raise ValueError(f"Path traversal detected in plot filename for accession: {accession}")
            accession_plot.save(plot_file_name)

__init__(prettydf, outdir)

GuideMakerPlot class for visualizing distrubution of gRNA, features, and locus.

Parameters:

Name	Type	Description	Default
prettydf	PandasDataFrame	Final output from GuideMaker	required
outdir	str	Output Directory	required

Returns:

Type	Description
None	None

Source code in guidemaker/core.py

def __init__(self, prettydf: pd.DataFrame, outdir: str) -> None:
    """
    GuideMakerPlot class for visualizing distrubution of gRNA, features, and locus.

    Args:
        prettydf (PandasDataFrame): Final output from GuideMaker
        outdir (str): Output Directory

    Returns:
        None
    """
    self.prettydf = prettydf
    self.accession = list(set(self.prettydf['Accession']))

    def _singleplot(df):
        """
        Returns guidemaker plot describing PAM targets

        Args:
            df(PandasDataFrame): Final output from GuideMaker for a single accession

        Return:
            None
        """
        source = df
        brush = alt.selection(type='interval', encodings=['x'])
        binNum = int(round(source['Feature end'].max() / 200, 0))
        display_info = source.columns.tolist()

        # Feature density
        densityF = alt.Chart(source).transform_density(
        'Feature start',
        as_=['Feature start', 'Feature Density'],
        extent=[1, source['Feature end'].max()],
        bandwidth=binNum,
        ).mark_area(color='black', opacity=0.6).encode(
        x=alt.X('Feature start', axis=alt.Axis(title='Genome Coordinates (bp)', tickCount=5)),
        y='Feature Density:Q',
        ).properties(height=50, width=500)

        # Guide density
        densityG = alt.Chart(source).transform_density(
        'Guide start',
        as_=['Guide start', 'Guide Density'],
        extent=[1, source['Feature end'].max()],
        bandwidth=binNum,
        ).mark_area(color='pink', opacity=0.6).encode(
        x=alt.X('Guide start', axis=alt.Axis(title='Genome Coordinates (bp)', tickCount=5)),
        y='Guide Density:Q',
        ).properties(height=50, width=500).add_selection(brush)

        # locus bar
        locus = alt.Chart(source).mark_bar(cornerRadiusTopLeft=3, cornerRadiusTopRight=3).encode(
        x='count(locus_tag):Q',
        y=alt.Y('locus_tag', axis=alt.Axis(title='Locus')),
        color='PAM:N',
        tooltip=display_info
        ).transform_filter(
        brush
        ).interactive().properties(height=500, width=500)
        guidemakerChart = (densityF & densityG & locus)
        return(guidemakerChart)

    outdir_real = os.path.realpath(outdir)
    for accession in self.accession:
        df = self.prettydf[self.prettydf['Accession'] == accession]
        accession_plot = _singleplot(df)
        clean_accession = re.sub(r'[^A-Za-z0-9._-]', '_', str(accession))
        plot_file_name = os.path.join(outdir, f"{clean_accession}.html")
        plot_file_real = os.path.realpath(plot_file_name)
        if not plot_file_real.startswith(outdir_real):
            raise ValueError(f"Path traversal detected in plot filename for accession: {accession}")
        accession_plot.save(plot_file_name)

PamTarget

A Class representing a Protospacer Adjacent Motif (PAM) and targets.

The classincludes all targets for given PAM as a dataframe,PAM and target attributes, and methods to find target and control sequences.

Source code in guidemaker/core.py

class PamTarget:

    """
    A Class representing a Protospacer Adjacent Motif (PAM) and targets.

    The classincludes all targets for given PAM as a dataframe,PAM and target attributes,
    and methods to find target and control sequences.

    """

    def __init__(self, pam: str, pam_orientation: str, dtype: str) -> None:
        """
        Pam __init__

        Args:
            pam (str): A DNA string in ambiguous IUPAC format
            pam_orientation (str): [5prime | 3prime ]
                5prime means the order is 5'-[pam][target]-3'
                3prime means the order is 5'-[target][pam]-3'
            dtype (str): hamming or leven

        Returns:
            None
        """
        for letter in pam.upper():
            assert letter in ['A', 'C', 'G', 'T', 'M', 'R', 'W',
                'S', 'Y', 'K', 'V', 'H', 'D', 'B', 'X', 'N']
        assert pam_orientation in ["3prime", "5prime"]
        self.pam: str = pam.upper()
        self.pam_orientation: str = pam_orientation
        self.dtype: str = dtype

    def __str__(self) -> str:
        """
        str __init__

        Args:
            self

        Returns:
            self(str)
        """
        return "A PAM object: {self.pam}".format(self=self)

    def find_targets(self, seq_record_iter: object, target_len: int) -> pd.DataFrame:
        """
        Find all targets on a sequence that match for the PAM on both strand(s)

        Args:
            seq_record_iter (object): A Biopython SeqRecord iterator from SeqIO.parse
            target_len (int): The length of the target sequence

        Returns:
            PandasDataFrame: A pandas dataframe with of matching targets
        """

        def reverse_complement(seq: str) -> str:
            """
            Reverse complement of the PAM sequence

            Args:
                seq (str): A DNA string

            Returns:
                str: A reverse complement of DNA string
            """
            bpseq = Seq.Seq(seq)
            return str(bpseq.reverse_complement())

        def pam2re(pam: str) -> str:
            """
            Convert an IUPAC ambiguous PAM to a Regex expression

            Args:
                pam (str): A DNA string

            Returns:
                str: A Regex expression
            """
            dnaval = {'A': 'A', 'C': 'C', 'G': 'G', 'T': 'T',
                      'M': '[A|C]', 'R': '[A|G]', 'W': '[A|T]', 'S': '[C|G]',
                      'Y': '[C|T]', 'K': '[G|T]', 'V': '[A|C|G]', 'H': '[A|C|T]',
                      'D': '[A|G|T]', 'B': '[C|G|T]', 'X': '[G|A|T|C]', 'N': '[G|A|T|C]'}
            return "".join([dnaval[base] for base in pam])

        #                5prime means the order is 5'-[pam][target]-3'
        #                3prime means the order is 5'-[target][pam]-3'

        def check_target(seq: str, target_len: int) -> bool:
            """
            Check targets for guidelength and DNA bases

            Args:
                seq (str): A DNA string
                target_len(int): Guide length

            Returns:
                bool: True or False
            """
            if len(seq) == target_len and all(letters in ['A', 'T', 'C', 'G'] for letters in seq):  # if not ATCG in the target then ignore those targets
                return True
            return False

        def run_for_5p(pam_pattern: str, dnaseq: str, target_len: int) -> Generator:
            """
            Search for guides with 5prime pam orientation in the forward strand

            Args:
                pam_pattern (str): A DNA string representing PAM
                dnaseq (str): A DNA string representing genome
                target_len (int): Guide length

            Returns:
                (Generator): A generator with target_seq, exact_pam, start, stop, strand, and pam_orientation
            """
            for match_obj in regex.finditer(pattern=pam_pattern, string=dnaseq, overlapped=True):
                target_seq = dnaseq[match_obj.end(): match_obj.end() + target_len]
                target_seq30 = dnaseq[match_obj.start()-3: match_obj.start()+27]
                ## 5'-[guide of 25 nt][exact pam, 3nt][next two]-3'
                if check_target(target_seq, target_len):
                    exact_pam = match_obj.group(0)
                    start = match_obj.end()
                    stop = match_obj.end() + target_len
                    # 5prime =True, 3prime = False
                    pam_orientation = True
                    # forward =True, reverse = False
                    strand = True
                    yield target_seq, exact_pam, start, stop, strand, pam_orientation, target_seq30



        def run_for_3p(pam_pattern, dnaseq, target_len) -> Generator:
            """
            Search for guides with 3prime pam orientation in the reverse strand

            Args:
                pam_pattern (str): A DNA string representing PAM
                dnaseq (str): A DNA string representing genome
                target_len (int): Guide length

            Returns:
                (Generator): A generator with target_seq, exact_pam, start, stop, strand, and pam_orientation
            """
            for match_obj in regex.finditer(pattern=pam_pattern, string=dnaseq, overlapped=True):
                target_seq = dnaseq[match_obj.start() - target_len: match_obj.start()]
                target_seq30 = dnaseq[match_obj.end()-27 :match_obj.end()+3]
                if check_target(target_seq, target_len):
                    exact_pam = match_obj.group(0)
                    start = match_obj.start() - target_len
                    stop = match_obj.start()
                    # 5prime =True, 3prime = False
                    pam_orientation = False
                    # forward =True, reverse = False
                    strand = True
                    yield target_seq, exact_pam, start, stop, strand, pam_orientation, target_seq30

        def run_rev_5p(pam_pattern, dnaseq, target_len) -> Generator:
            """
            Search for guides with 5prime pam orientation in the reverse strand

            Args:
                pam_pattern (str): A DNA string representing PAM
                dnaseq (str): A DNA string representing genome
                target_len (int): Guide length

            Returns:
                (Generator): A generator with target_seq, exact_pam, start, stop, strand, and pam_orientation
            """
            for match_obj in regex.finditer(pattern=pam_pattern, string=dnaseq, overlapped=True):
                target_seq = reverse_complement(
                    dnaseq[match_obj.start() - target_len: match_obj.start()])
                target_seq30 = reverse_complement(
                    dnaseq[match_obj.end()-27:match_obj.end()+3])
                if check_target(target_seq, target_len):
                    exact_pam = reverse_complement(match_obj.group(0))
                    start = match_obj.start() - target_len
                    stop = match_obj.start()
                    # 5prime =True, 3prime = False
                    pam_orientation = True
                    # forward =True, reverse = False
                    strand = False
                    yield target_seq, exact_pam, start, stop, strand, pam_orientation, target_seq30

        def run_rev_3p(pam_pattern, dnaseq, target_len) -> Generator:
            """
            Search for guides with 3prime pam orientation in the reverse strand

            Args:
                pam_pattern (str): A DNA string representing PAM
                dnaseq (str): A DNA string representing genome
                target_len (int): Guide length

            Returns:
                (Generator): A generator with target_seq, exact_pam, start, stop, strand, and pam_orientation
            """
            for match_obj in regex.finditer(pattern=pam_pattern, string=dnaseq, overlapped=True):
                target_seq = reverse_complement(
                    dnaseq[match_obj.end(): match_obj.end() + target_len])
                target_seq30 = reverse_complement(dnaseq[match_obj.start()-3:match_obj.start()+27])
                if check_target(target_seq, target_len):
                    exact_pam = reverse_complement(match_obj.group(0))
                    start = match_obj.end()
                    stop = match_obj.end() + target_len
                    # 5prime =True, 3prime = False
                    pam_orientation = False
                    # forward =True, reverse = False
                    strand = False
                    yield target_seq, exact_pam, start, stop, strand, pam_orientation, target_seq30

        target_list = []
        for record in seq_record_iter:
            record_id = record.id
            seq = str(record.seq)
            if self.pam_orientation == "5prime":
                # forward
                for5p = pd.DataFrame(run_for_5p(pam2re(self.pam), seq, target_len), columns=[
                                     "target", "exact_pam", "start", "stop", "strand", "pam_orientation", "target_seq30"])
                for5p["seqid"] = record_id
                # string to boolean conversion is not straight - as all string were set to Trues- so change the encoding in functions above.
                # https://stackoverflow.com/questions/715417/converting-from-a-string-to-boolean-in-python/715455#715455
                for5p = for5p.astype({"target": 'str', "exact_pam": 'category', "start": 'uint32',
                                     "stop": 'uint32', "strand": 'bool', "pam_orientation": 'bool', "seqid": 'category'})
                target_list.append(for5p)
                # reverse
                rev5p = pd.DataFrame(run_rev_5p(pam2re(reverse_complement(self.pam)), seq, target_len), columns=[
                                     "target", "exact_pam", "start", "stop", "strand", "pam_orientation","target_seq30"])
                rev5p["seqid"] = record_id
                rev5p = rev5p.astype({"target": 'str', "exact_pam": 'category', "start": 'uint32',
                                     "stop": 'uint32', "strand": 'bool', "pam_orientation": 'bool', "seqid": 'category'})
                target_list.append(rev5p)
                # Question? Append directly vs. concat then append? https://ravinpoudel.github.io/AppendVsConcat/
            elif self.pam_orientation == "3prime":
                # forward
                for3p = pd.DataFrame(run_for_3p(pam2re(self.pam), seq, target_len), columns=[
                                     "target", "exact_pam", "start", "stop", "strand", "pam_orientation","target_seq30"])
                for3p["seqid"] = record_id
                for3p = for3p.astype({"target": 'str', "exact_pam": 'category', "start": 'uint32',
                                     "stop": 'uint32', "strand": 'bool', "pam_orientation": 'bool', "seqid": 'category'})
                target_list.append(for3p)
                # reverse
                rev3p = pd.DataFrame(run_rev_3p(pam2re(reverse_complement(self.pam)), seq, target_len), columns=[
                                     "target", "exact_pam", "start", "stop", "strand", "pam_orientation","target_seq30"])
                rev3p["seqid"] = record_id
                rev3p = rev3p.astype({"target": 'str', "exact_pam": 'category', "start": 'uint32',
                                     "stop": 'uint32', "strand": 'bool', "pam_orientation": 'bool', "seqid": 'category'})
                target_list.append(rev3p)
            gc.collect()  # clear memory after each chromosome
        target_list = [item for item in target_list if not item.empty]
        df_targets = pd.concat(target_list, ignore_index=True)
        df_targets = df_targets.assign(seedseq=np.nan, hasrestrictionsite=np.nan, isseedduplicated=np.nan)
        df_targets = df_targets.astype({"seedseq": 'str', "isseedduplicated": 'bool'})
        df_targets = df_targets.assign(dtype=self.dtype)
        df_targets = df_targets.astype({"dtype": 'category'})
        return df_targets

__init__(pam, pam_orientation, dtype)

Pam init

Parameters:

Name	Type	Description	Default
pam	str	A DNA string in ambiguous IUPAC format	required
pam_orientation	str	[5prime | 3prime ] 5prime means the order is 5'-[pam][target]-3' 3prime means the order is 5'-[target][pam]-3'	required
dtype	str	hamming or leven	required

Returns:

Type	Description
None	None

Source code in guidemaker/core.py

def __init__(self, pam: str, pam_orientation: str, dtype: str) -> None:
    """
    Pam __init__

    Args:
        pam (str): A DNA string in ambiguous IUPAC format
        pam_orientation (str): [5prime | 3prime ]
            5prime means the order is 5'-[pam][target]-3'
            3prime means the order is 5'-[target][pam]-3'
        dtype (str): hamming or leven

    Returns:
        None
    """
    for letter in pam.upper():
        assert letter in ['A', 'C', 'G', 'T', 'M', 'R', 'W',
            'S', 'Y', 'K', 'V', 'H', 'D', 'B', 'X', 'N']
    assert pam_orientation in ["3prime", "5prime"]
    self.pam: str = pam.upper()
    self.pam_orientation: str = pam_orientation
    self.dtype: str = dtype

__str__()

str init

Returns:

Type	Description
str	self(str)

Source code in guidemaker/core.py

def __str__(self) -> str:
    """
    str __init__

    Args:
        self

    Returns:
        self(str)
    """
    return "A PAM object: {self.pam}".format(self=self)

find_targets(seq_record_iter, target_len)

Find all targets on a sequence that match for the PAM on both strand(s)

Parameters:

Name	Type	Description	Default
seq_record_iter	object	A Biopython SeqRecord iterator from SeqIO.parse	required
target_len	int	The length of the target sequence	required

Returns:

Name	Type	Description
PandasDataFrame	DataFrame	A pandas dataframe with of matching targets

Source code in guidemaker/core.py

def find_targets(self, seq_record_iter: object, target_len: int) -> pd.DataFrame:
    """
    Find all targets on a sequence that match for the PAM on both strand(s)

    Args:
        seq_record_iter (object): A Biopython SeqRecord iterator from SeqIO.parse
        target_len (int): The length of the target sequence

    Returns:
        PandasDataFrame: A pandas dataframe with of matching targets
    """

    def reverse_complement(seq: str) -> str:
        """
        Reverse complement of the PAM sequence

        Args:
            seq (str): A DNA string

        Returns:
            str: A reverse complement of DNA string
        """
        bpseq = Seq.Seq(seq)
        return str(bpseq.reverse_complement())

    def pam2re(pam: str) -> str:
        """
        Convert an IUPAC ambiguous PAM to a Regex expression

        Args:
            pam (str): A DNA string

        Returns:
            str: A Regex expression
        """
        dnaval = {'A': 'A', 'C': 'C', 'G': 'G', 'T': 'T',
                  'M': '[A|C]', 'R': '[A|G]', 'W': '[A|T]', 'S': '[C|G]',
                  'Y': '[C|T]', 'K': '[G|T]', 'V': '[A|C|G]', 'H': '[A|C|T]',
                  'D': '[A|G|T]', 'B': '[C|G|T]', 'X': '[G|A|T|C]', 'N': '[G|A|T|C]'}
        return "".join([dnaval[base] for base in pam])

    #                5prime means the order is 5'-[pam][target]-3'
    #                3prime means the order is 5'-[target][pam]-3'

    def check_target(seq: str, target_len: int) -> bool:
        """
        Check targets for guidelength and DNA bases

        Args:
            seq (str): A DNA string
            target_len(int): Guide length

        Returns:
            bool: True or False
        """
        if len(seq) == target_len and all(letters in ['A', 'T', 'C', 'G'] for letters in seq):  # if not ATCG in the target then ignore those targets
            return True
        return False

    def run_for_5p(pam_pattern: str, dnaseq: str, target_len: int) -> Generator:
        """
        Search for guides with 5prime pam orientation in the forward strand

        Args:
            pam_pattern (str): A DNA string representing PAM
            dnaseq (str): A DNA string representing genome
            target_len (int): Guide length

        Returns:
            (Generator): A generator with target_seq, exact_pam, start, stop, strand, and pam_orientation
        """
        for match_obj in regex.finditer(pattern=pam_pattern, string=dnaseq, overlapped=True):
            target_seq = dnaseq[match_obj.end(): match_obj.end() + target_len]
            target_seq30 = dnaseq[match_obj.start()-3: match_obj.start()+27]
            ## 5'-[guide of 25 nt][exact pam, 3nt][next two]-3'
            if check_target(target_seq, target_len):
                exact_pam = match_obj.group(0)
                start = match_obj.end()
                stop = match_obj.end() + target_len
                # 5prime =True, 3prime = False
                pam_orientation = True
                # forward =True, reverse = False
                strand = True
                yield target_seq, exact_pam, start, stop, strand, pam_orientation, target_seq30



    def run_for_3p(pam_pattern, dnaseq, target_len) -> Generator:
        """
        Search for guides with 3prime pam orientation in the reverse strand

        Args:
            pam_pattern (str): A DNA string representing PAM
            dnaseq (str): A DNA string representing genome
            target_len (int): Guide length

        Returns:
            (Generator): A generator with target_seq, exact_pam, start, stop, strand, and pam_orientation
        """
        for match_obj in regex.finditer(pattern=pam_pattern, string=dnaseq, overlapped=True):
            target_seq = dnaseq[match_obj.start() - target_len: match_obj.start()]
            target_seq30 = dnaseq[match_obj.end()-27 :match_obj.end()+3]
            if check_target(target_seq, target_len):
                exact_pam = match_obj.group(0)
                start = match_obj.start() - target_len
                stop = match_obj.start()
                # 5prime =True, 3prime = False
                pam_orientation = False
                # forward =True, reverse = False
                strand = True
                yield target_seq, exact_pam, start, stop, strand, pam_orientation, target_seq30

    def run_rev_5p(pam_pattern, dnaseq, target_len) -> Generator:
        """
        Search for guides with 5prime pam orientation in the reverse strand

        Args:
            pam_pattern (str): A DNA string representing PAM
            dnaseq (str): A DNA string representing genome
            target_len (int): Guide length

        Returns:
            (Generator): A generator with target_seq, exact_pam, start, stop, strand, and pam_orientation
        """
        for match_obj in regex.finditer(pattern=pam_pattern, string=dnaseq, overlapped=True):
            target_seq = reverse_complement(
                dnaseq[match_obj.start() - target_len: match_obj.start()])
            target_seq30 = reverse_complement(
                dnaseq[match_obj.end()-27:match_obj.end()+3])
            if check_target(target_seq, target_len):
                exact_pam = reverse_complement(match_obj.group(0))
                start = match_obj.start() - target_len
                stop = match_obj.start()
                # 5prime =True, 3prime = False
                pam_orientation = True
                # forward =True, reverse = False
                strand = False
                yield target_seq, exact_pam, start, stop, strand, pam_orientation, target_seq30

    def run_rev_3p(pam_pattern, dnaseq, target_len) -> Generator:
        """
        Search for guides with 3prime pam orientation in the reverse strand

        Args:
            pam_pattern (str): A DNA string representing PAM
            dnaseq (str): A DNA string representing genome
            target_len (int): Guide length

        Returns:
            (Generator): A generator with target_seq, exact_pam, start, stop, strand, and pam_orientation
        """
        for match_obj in regex.finditer(pattern=pam_pattern, string=dnaseq, overlapped=True):
            target_seq = reverse_complement(
                dnaseq[match_obj.end(): match_obj.end() + target_len])
            target_seq30 = reverse_complement(dnaseq[match_obj.start()-3:match_obj.start()+27])
            if check_target(target_seq, target_len):
                exact_pam = reverse_complement(match_obj.group(0))
                start = match_obj.end()
                stop = match_obj.end() + target_len
                # 5prime =True, 3prime = False
                pam_orientation = False
                # forward =True, reverse = False
                strand = False
                yield target_seq, exact_pam, start, stop, strand, pam_orientation, target_seq30

    target_list = []
    for record in seq_record_iter:
        record_id = record.id
        seq = str(record.seq)
        if self.pam_orientation == "5prime":
            # forward
            for5p = pd.DataFrame(run_for_5p(pam2re(self.pam), seq, target_len), columns=[
                                 "target", "exact_pam", "start", "stop", "strand", "pam_orientation", "target_seq30"])
            for5p["seqid"] = record_id
            # string to boolean conversion is not straight - as all string were set to Trues- so change the encoding in functions above.
            # https://stackoverflow.com/questions/715417/converting-from-a-string-to-boolean-in-python/715455#715455
            for5p = for5p.astype({"target": 'str', "exact_pam": 'category', "start": 'uint32',
                                 "stop": 'uint32', "strand": 'bool', "pam_orientation": 'bool', "seqid": 'category'})
            target_list.append(for5p)
            # reverse
            rev5p = pd.DataFrame(run_rev_5p(pam2re(reverse_complement(self.pam)), seq, target_len), columns=[
                                 "target", "exact_pam", "start", "stop", "strand", "pam_orientation","target_seq30"])
            rev5p["seqid"] = record_id
            rev5p = rev5p.astype({"target": 'str', "exact_pam": 'category', "start": 'uint32',
                                 "stop": 'uint32', "strand": 'bool', "pam_orientation": 'bool', "seqid": 'category'})
            target_list.append(rev5p)
            # Question? Append directly vs. concat then append? https://ravinpoudel.github.io/AppendVsConcat/
        elif self.pam_orientation == "3prime":
            # forward
            for3p = pd.DataFrame(run_for_3p(pam2re(self.pam), seq, target_len), columns=[
                                 "target", "exact_pam", "start", "stop", "strand", "pam_orientation","target_seq30"])
            for3p["seqid"] = record_id
            for3p = for3p.astype({"target": 'str', "exact_pam": 'category', "start": 'uint32',
                                 "stop": 'uint32', "strand": 'bool', "pam_orientation": 'bool', "seqid": 'category'})
            target_list.append(for3p)
            # reverse
            rev3p = pd.DataFrame(run_rev_3p(pam2re(reverse_complement(self.pam)), seq, target_len), columns=[
                                 "target", "exact_pam", "start", "stop", "strand", "pam_orientation","target_seq30"])
            rev3p["seqid"] = record_id
            rev3p = rev3p.astype({"target": 'str', "exact_pam": 'category', "start": 'uint32',
                                 "stop": 'uint32', "strand": 'bool', "pam_orientation": 'bool', "seqid": 'category'})
            target_list.append(rev3p)
        gc.collect()  # clear memory after each chromosome
    target_list = [item for item in target_list if not item.empty]
    df_targets = pd.concat(target_list, ignore_index=True)
    df_targets = df_targets.assign(seedseq=np.nan, hasrestrictionsite=np.nan, isseedduplicated=np.nan)
    df_targets = df_targets.astype({"seedseq": 'str', "isseedduplicated": 'bool'})
    df_targets = df_targets.assign(dtype=self.dtype)
    df_targets = df_targets.astype({"dtype": 'category'})
    return df_targets

TargetProcessor

A Class representing a set of guide RNA targets.

The class includes all targets in a dataframe, methods to process target and a dict with edit distances for sequences.

Source code in guidemaker/core.py

class TargetProcessor:

    """
    A Class representing a set of guide RNA targets.

    The class includes all targets in a dataframe, methods to process target and a dict with edit distances for sequences.

    """

    def __init__(self, targets: pd.DataFrame, lsr: int, editdist: int = 2, knum: int = 2) -> None:
        """
        TargetProcessor __init__

        Args:
            targets (PandasDataFrame): Dataframe with output from class PamTarget
            lsr (int): Length of seed region
            editdist (int): Edit distance
            knum (int): Number of negative controls

        Returns:
            None
        """
        self.targets = targets  # pandas dataframe
        self.lsr: int = lsr  # length of seed region
        self.editdist: int = editdist
        self.knum: int = knum
        self.nmslib_index: object = None
        self.neighbors: dict = {}
        self.closest_neighbor_df: pd.DataFrame = None
        self.ncontrolsearched: int = None
        self.gc_percent: float = None
        self.genomesize: float = None
        self.pam_orientation: bool = targets['pam_orientation'].iat[0]

    def __str__(self) -> None:
        """
        str __init__

        Args:
            self

        Return:
            None
        """
        info = "TargetList: contains a set of {} potential PAM targets".format(len(self.targets))
        return info

    def __len__(self) -> int:
        """
        len __init__ to display length of self.targets

        Args:
            self.targets

        Return:
            (int): Length of the self.targets
        """
        return len(self.targets)

    def check_restriction_enzymes(self, restriction_enzyme_list: list = []) -> None:
        """
        Check for restriction enzymes and its reverse complement within gRNA sequence

        Args:
            restriction_enzyme_list (list): A list with sequence for restriction enzymes

        Returns:
            None
        """
        element_to_exclude = []
        for record in set(restriction_enzyme_list):
            for letter in record.upper():
                assert letter in ['A', 'C', 'G', 'T', 'M', 'R', 'W',
                    'S', 'Y', 'K', 'V', 'H', 'D', 'B', 'X', 'N']
            record_seq = Seq.Seq(record.upper())
            element_to_exclude.append(extend_ambiguous_dna(str(record_seq)))
            element_to_exclude.append(extend_ambiguous_dna(
                str(record_seq.reverse_complement())))  # reverse complement
        element_to_exclude = sum(element_to_exclude, [])  # flatout list of list to list with restriction enzyme sites
        if len(element_to_exclude) > 0:
            self.targets['hasrestrictionsite'] = self.targets['target'].str.contains('|'.join(element_to_exclude))
        else:
            self.targets['hasrestrictionsite'] = False

    def _one_hot_encode(self, seq_list: List[object]) -> List[str]:
        """One hot encode Target DNA as a binary string representation for NMSLIB."""
        charmap = {'A': '1 0 0 0', 'C': '0 1 0 0', 'G': '0 0 1 0', 'T': '0 0 0 1'}

        def seq_to_bin(seq):
            charlist = [charmap[letter] for letter in seq]
            return " ".join(charlist)
        return list(map(seq_to_bin, seq_list))

    def find_unique_near_pam(self) -> None:
        """
        Identify unique sequences in the target list

        The function filters a list of Target objects for targets that
        are unique in the region closest to the PAM. The region length is defined
        by the lsr (length of seed region that need to be unique).

        Args:
            lsr (int): Length of seed region that is close to PAM

        Returns:
            None
        """
        def _get_prox(tseq):  # get target sequence as input
            if self.pam_orientation == True:  # 5prime = True 3prime=False
                if self.lsr == 0:
                    return tseq
                else:
                    return tseq[0:self.lsr]
            elif self.pam_orientation == False:  # 5prime = True 3prime=False
                if self.lsr == 0:
                    return tseq
                else:
                    return tseq[(len(tseq) - self.lsr):]
        # https://stackoverflow.com/questions/12555323/adding-new-column-to-existing-dataframe-in-python-pandas
        self.targets = deepcopy(self.targets)
        self.targets['seedseq'] = self.targets['target'].apply(_get_prox)
        self.targets['isseedduplicated'] = self.targets['seedseq'].duplicated()

    def create_index(self, configpath: str, num_threads=2):
        """
        Create nmslib index

        Converts self.targets to binary one hot encoding and returns NMSLIB index

        Args:
            num_threads (int): cpu threads
            configpath (str): Path to config file which contains hyper parameters for NMSLIB

                M (int): Controls the number of bi-directional links created for each element
                during index construction. Higher values lead to better results at the expense
                of memory consumption. Typical values are 2 -100, but for most datasets a
                range of 12 -48 is suitable. Can’t be smaller than 2.

                efC (int): Size of the dynamic list used during construction. A larger value means
                   a better quality index, but increases build time. Should be an integer value
                   between 1 and the size of the dataset.

        Returns:
            None (but writes NMSLIB index to self)
        """
        with open(configpath) as cf:
            config = yaml.safe_load(cf)

        M, efC, post = config['NMSLIB']['M'], config['NMSLIB']['efc'], config['NMSLIB']['post']

        # index everything but not duplicates
        self.notduplicated_targets = list(set(self.targets['target'].tolist()))
        notduplicated_targets = self.notduplicated_targets
        #mod_logger.info("unique targets for index: %s" % len(notduplicated_targets))
        if self.targets['dtype'].iat[0] == "hamming":
            bintargets = self._one_hot_encode(notduplicated_targets)
            index_params = {'M': M, 'indexThreadQty': num_threads, 'efConstruction': efC, 'post': post}
            index = nmslib.init(space='bit_hamming',
                            dtype=nmslib.DistType.INT,
                            data_type=nmslib.DataType.OBJECT_AS_STRING,
                            method='hnsw')
            index.addDataPointBatch(bintargets) # notduplicated_targets
            index.createIndex(index_params, print_progress=True)
            self.nmslib_index = index
        else:
            bintargets = notduplicated_targets
            index_params = {'M': M, 'indexThreadQty': num_threads, 'efConstruction': efC, 'post': post}
            index = nmslib.init(space='leven',
                            dtype=nmslib.DistType.INT,
                            data_type=nmslib.DataType.OBJECT_AS_STRING,
                            method='hnsw')
            index.addDataPointBatch(bintargets) # notduplicated_targets
            index.createIndex(index_params, print_progress=True)
            self.nmslib_index = index



    def get_neighbors(self, configpath, num_threads=2) -> None:
        """
        Get nearest neighbors for sequences removing sequences that
        have neighbors less than the Hamming distance threshold.
        For the list of all targets calculate the (knum) nearest neighbors.
        filter out targets with close neighbors and
        Writes a dictionary to self.neighbors:
        self.neighbors[seq]{target: seq_obj, neighbors: {seqs:[s1, s1, ...], dist:[d1, d1,...]}}

        Args:
            configpath (str): Path to a parameter config file
            num_threads (int): Number of threads

        Returns:
            None
        """
        with open(configpath) as cf:
            config = yaml.safe_load(cf)

        ef = config['NMSLIB']['ef']

        # unique_targets = self.targets.loc[self.targets['isseedduplicated']
        #     == False]['target'].tolist()
        # For indexing we need to use all targets -- for checking off-targets. For searching neighbours remove seed duplicated and one wiht restriction site.
        unique_targets = self.targets.loc[(self.targets['isseedduplicated']==False) & (self.targets['hasrestrictionsite']==False)]['target'].tolist()
        if self.targets['dtype'].iat[0] == "hamming":
            unique_bintargets = self._one_hot_encode(unique_targets)  # search unique seed one
        else:
            unique_bintargets = unique_targets

        self.nmslib_index.setQueryTimeParams({'efSearch': ef})
        results_list = self.nmslib_index.knnQueryBatch(unique_bintargets,
                                               k=self.knum, num_threads=num_threads)
        neighbor_dict = {}
        target_lookup = getattr(self, 'notduplicated_targets', self.targets['target'].values)
        for i, entry in enumerate(results_list):
            queryseq = unique_targets[i]
            hitseqidx = entry[0].tolist()
            editdist = entry[1].tolist()
            if self.targets['dtype'].iat[0] == "hamming":
                # check that the closest sequence meets the min. dist. requirment. We multiply by 2 b/c each 
                # base is one hot encoded. e.g. 1000 vs 0100 = 2 differences
                if editdist[1] >= 2 * self.editdist:
                    neighbors = {"seqs": [target_lookup[x] for x in hitseqidx],
                                "dist": [int(x / 2) for x in editdist]} 
                    neighbor_dict[queryseq] = {"target": unique_targets[i],
                                            "neighbors": neighbors}
            else:
               if editdist[1] >= self.editdist: 
                    neighbors = {"seqs": [target_lookup[x] for x in hitseqidx],
                                "dist": [int(x) for x in editdist]}
                    neighbor_dict[queryseq] = {"target": unique_targets[i],
                                            "neighbors": neighbors}
        self.neighbors = neighbor_dict

    def export_bed(self) -> object:
        """
        Export the targets in self.neighbors to a bed format file

        Args:
            file (str): the name and location of file to export

        Returns:
            (obj): A Pandas Dataframe in Bed format
        """
        # df = self.targets.copy()
        # why deepcopy - https://stackoverflow.com/questions/55745948/why-doesnt-deepcopy-of-a-pandas-dataframe-affect-memory-usage
        # select only guides that are not duplecated in the seedseq
        df = deepcopy(self.targets.loc[self.targets['isseedduplicated'] == False])
        df = df[["seqid", "start", "stop", "target", "strand"]]
        strand_col = df['strand'].apply(lambda x: '+' if x in [True, '+', 'true', 'True'] else '-')
        df = df.assign(strand=strand_col)
        df.columns = ["chrom", "chromstart", "chromend", "name", "strand"]
        df.sort_values(by=['chrom', 'chromstart'], inplace=True)
        return df

    def get_control_seqs(self, seq_record_iter: object, configpath, length: int = 20, n: int = 10,
                         num_threads: int = 2) -> pd.DataFrame:
        """
        Create random sequences with a specified GC probability and find seqs with the greatest
        distance to any sequence flanking a PAM site

        Args:
            seq_record_iter (Bio.SeqIO): An iterator of fastas
            length (int): Length of the sequence, must match the index
            n (int): Number of sequences to  return
            num_threads (int): Number of processor threads

        Returns:
            (PandasDataFrame): A pandas dataframe with control sequence
        """

        with open(configpath) as cf:
            config = yaml.safe_load(cf)

        MINIMUM_HMDIST = config['CONTROL']['MINIMUM_HMDIST']

        MAX_CONTROL_SEARCH_MULTIPLE = max(config['CONTROL']['CONTROL_SEARCH_MULTIPLE'])

        #  search_mult (int): search this times n sequences
        CONTROL_SEARCH_MULTIPLE = config['CONTROL']['CONTROL_SEARCH_MULTIPLE']

        # get GC percent
        totlen = 0
        gccnt = 0
        for record in seq_record_iter:
            gccnt += gc_fraction(record.seq) * len(record)
            totlen += len(record)
        gc = gccnt / (totlen)
        self.gc_percent = gc * 100
        self.genomesize = totlen / (1024 * 1024)

        minimum_hmdist = 0
        sm_count = 0
        search_mult = 0

        try:
            while minimum_hmdist < MINIMUM_HMDIST or search_mult == MAX_CONTROL_SEARCH_MULTIPLE:
                # generate random sequences
                seqs = []
                search_mult = CONTROL_SEARCH_MULTIPLE[sm_count]
                for i in range(n * search_mult):
                    seqs.append("".join(np.random.choice(a=["G", "C", "A", "T"], size=length,
                                                         replace=True, p=[gc / 2, gc / 2, (1 - gc) / 2, (1 - gc) / 2])))
                # one hot encode sequences
                binseq = []
                charmap = {'A': '1 0 0 0', 'C': '0 1 0 0', 'G': '0 0 1 0', 'T': '0 0 0 1'}
                for seq in seqs:
                    if self.targets['dtype'].iat[0] == "hamming":
                        charlist = [charmap[letter] for letter in seq]
                        binseq.append(" ".join(charlist))
                    else: # leven
                        binseq.append(seq)

                rand_seqs = self.nmslib_index.knnQueryBatch(binseq, k=2, num_threads=num_threads)
                distlist = []
                for i in rand_seqs:
                    distlist.append(i[1][0])
                zipped = list(zip(seqs, distlist))
                dist_seqs = sorted(zipped, reverse=True, key=lambda x: x[1])
                sort_seq = [item[0] for item in dist_seqs][0:n]

                #sort_dist
                if self.targets['dtype'].iat[0] == "hamming":
                    sort_dist = [item[1] / 2 for item in dist_seqs][0:n]  ### ? does divide by 2 holds for leven???
                else:
                    sort_dist = [item[1] for item in dist_seqs][0:n]  ### ? does divide by 2 holds for leven???

                minimum_hmdist = int(min(sort_dist))
                sm_count += 1
        except IndexError as e:
            raise e

        total_ncontrolsearched = search_mult * n
        self.ncontrolsearched = total_ncontrolsearched
        randomdf = pd.DataFrame(data={"Sequences": sort_seq, "Hamming distance": sort_dist})

        def create_name(seq):
            return "Cont-" + hashlib.md5(seq.encode()).hexdigest()
        randomdf['name'] = randomdf["Sequences"].apply(create_name)
        randomdf = randomdf[["name", "Sequences", "Hamming distance"]]
        randomdf.head()
        return (min(sort_dist),
                statistics.median(sort_dist),
                randomdf)

__init__(targets, lsr, editdist=2, knum=2)

TargetProcessor init

Parameters:

Name	Type	Description	Default
targets	PandasDataFrame	Dataframe with output from class PamTarget	required
lsr	int	Length of seed region	required
editdist	int	Edit distance	2
knum	int	Number of negative controls	2

Returns:

Type	Description
None	None

Source code in guidemaker/core.py

def __init__(self, targets: pd.DataFrame, lsr: int, editdist: int = 2, knum: int = 2) -> None:
    """
    TargetProcessor __init__

    Args:
        targets (PandasDataFrame): Dataframe with output from class PamTarget
        lsr (int): Length of seed region
        editdist (int): Edit distance
        knum (int): Number of negative controls

    Returns:
        None
    """
    self.targets = targets  # pandas dataframe
    self.lsr: int = lsr  # length of seed region
    self.editdist: int = editdist
    self.knum: int = knum
    self.nmslib_index: object = None
    self.neighbors: dict = {}
    self.closest_neighbor_df: pd.DataFrame = None
    self.ncontrolsearched: int = None
    self.gc_percent: float = None
    self.genomesize: float = None
    self.pam_orientation: bool = targets['pam_orientation'].iat[0]

__len__()

len init to display length of self.targets

Return

(int): Length of the self.targets

Source code in guidemaker/core.py

def __len__(self) -> int:
    """
    len __init__ to display length of self.targets

    Args:
        self.targets

    Return:
        (int): Length of the self.targets
    """
    return len(self.targets)

__str__()

str init

Return

None

Source code in guidemaker/core.py

def __str__(self) -> None:
    """
    str __init__

    Args:
        self

    Return:
        None
    """
    info = "TargetList: contains a set of {} potential PAM targets".format(len(self.targets))
    return info

check_restriction_enzymes(restriction_enzyme_list=[])

Check for restriction enzymes and its reverse complement within gRNA sequence

Parameters:

Name	Type	Description	Default
restriction_enzyme_list	list	A list with sequence for restriction enzymes	[]

Returns:

Type	Description
None	None

Source code in guidemaker/core.py

def check_restriction_enzymes(self, restriction_enzyme_list: list = []) -> None:
    """
    Check for restriction enzymes and its reverse complement within gRNA sequence

    Args:
        restriction_enzyme_list (list): A list with sequence for restriction enzymes

    Returns:
        None
    """
    element_to_exclude = []
    for record in set(restriction_enzyme_list):
        for letter in record.upper():
            assert letter in ['A', 'C', 'G', 'T', 'M', 'R', 'W',
                'S', 'Y', 'K', 'V', 'H', 'D', 'B', 'X', 'N']
        record_seq = Seq.Seq(record.upper())
        element_to_exclude.append(extend_ambiguous_dna(str(record_seq)))
        element_to_exclude.append(extend_ambiguous_dna(
            str(record_seq.reverse_complement())))  # reverse complement
    element_to_exclude = sum(element_to_exclude, [])  # flatout list of list to list with restriction enzyme sites
    if len(element_to_exclude) > 0:
        self.targets['hasrestrictionsite'] = self.targets['target'].str.contains('|'.join(element_to_exclude))
    else:
        self.targets['hasrestrictionsite'] = False

create_index(configpath, num_threads=2)

Create nmslib index

Converts self.targets to binary one hot encoding and returns NMSLIB index

Parameters:

Name	Type	Description	Default
num_threads	int	cpu threads	2
configpath	str	Path to config file which contains hyper parameters for NMSLIB M (int): Controls the number of bi-directional links created for each element during index construction. Higher values lead to better results at the expense of memory consumption. Typical values are 2 -100, but for most datasets a range of 12 -48 is suitable. Can’t be smaller than 2. efC (int): Size of the dynamic list used during construction. A larger value means a better quality index, but increases build time. Should be an integer value between 1 and the size of the dataset.	required

Returns:

Type	Description
	None (but writes NMSLIB index to self)

Source code in guidemaker/core.py

def create_index(self, configpath: str, num_threads=2):
    """
    Create nmslib index

    Converts self.targets to binary one hot encoding and returns NMSLIB index

    Args:
        num_threads (int): cpu threads
        configpath (str): Path to config file which contains hyper parameters for NMSLIB

            M (int): Controls the number of bi-directional links created for each element
            during index construction. Higher values lead to better results at the expense
            of memory consumption. Typical values are 2 -100, but for most datasets a
            range of 12 -48 is suitable. Can’t be smaller than 2.

            efC (int): Size of the dynamic list used during construction. A larger value means
               a better quality index, but increases build time. Should be an integer value
               between 1 and the size of the dataset.

    Returns:
        None (but writes NMSLIB index to self)
    """
    with open(configpath) as cf:
        config = yaml.safe_load(cf)

    M, efC, post = config['NMSLIB']['M'], config['NMSLIB']['efc'], config['NMSLIB']['post']

    # index everything but not duplicates
    self.notduplicated_targets = list(set(self.targets['target'].tolist()))
    notduplicated_targets = self.notduplicated_targets
    #mod_logger.info("unique targets for index: %s" % len(notduplicated_targets))
    if self.targets['dtype'].iat[0] == "hamming":
        bintargets = self._one_hot_encode(notduplicated_targets)
        index_params = {'M': M, 'indexThreadQty': num_threads, 'efConstruction': efC, 'post': post}
        index = nmslib.init(space='bit_hamming',
                        dtype=nmslib.DistType.INT,
                        data_type=nmslib.DataType.OBJECT_AS_STRING,
                        method='hnsw')
        index.addDataPointBatch(bintargets) # notduplicated_targets
        index.createIndex(index_params, print_progress=True)
        self.nmslib_index = index
    else:
        bintargets = notduplicated_targets
        index_params = {'M': M, 'indexThreadQty': num_threads, 'efConstruction': efC, 'post': post}
        index = nmslib.init(space='leven',
                        dtype=nmslib.DistType.INT,
                        data_type=nmslib.DataType.OBJECT_AS_STRING,
                        method='hnsw')
        index.addDataPointBatch(bintargets) # notduplicated_targets
        index.createIndex(index_params, print_progress=True)
        self.nmslib_index = index

export_bed()

Export the targets in self.neighbors to a bed format file

Parameters:

Name	Type	Description	Default
file	str	the name and location of file to export	required

Returns:

Type	Description
obj	A Pandas Dataframe in Bed format

Source code in guidemaker/core.py

def export_bed(self) -> object:
    """
    Export the targets in self.neighbors to a bed format file

    Args:
        file (str): the name and location of file to export

    Returns:
        (obj): A Pandas Dataframe in Bed format
    """
    # df = self.targets.copy()
    # why deepcopy - https://stackoverflow.com/questions/55745948/why-doesnt-deepcopy-of-a-pandas-dataframe-affect-memory-usage
    # select only guides that are not duplecated in the seedseq
    df = deepcopy(self.targets.loc[self.targets['isseedduplicated'] == False])
    df = df[["seqid", "start", "stop", "target", "strand"]]
    strand_col = df['strand'].apply(lambda x: '+' if x in [True, '+', 'true', 'True'] else '-')
    df = df.assign(strand=strand_col)
    df.columns = ["chrom", "chromstart", "chromend", "name", "strand"]
    df.sort_values(by=['chrom', 'chromstart'], inplace=True)
    return df

find_unique_near_pam()

Identify unique sequences in the target list

The function filters a list of Target objects for targets that are unique in the region closest to the PAM. The region length is defined by the lsr (length of seed region that need to be unique).

Parameters:

Name	Type	Description	Default
lsr	int	Length of seed region that is close to PAM	required

Returns:

Type	Description
None	None

Source code in guidemaker/core.py

def find_unique_near_pam(self) -> None:
    """
    Identify unique sequences in the target list

    The function filters a list of Target objects for targets that
    are unique in the region closest to the PAM. The region length is defined
    by the lsr (length of seed region that need to be unique).

    Args:
        lsr (int): Length of seed region that is close to PAM

    Returns:
        None
    """
    def _get_prox(tseq):  # get target sequence as input
        if self.pam_orientation == True:  # 5prime = True 3prime=False
            if self.lsr == 0:
                return tseq
            else:
                return tseq[0:self.lsr]
        elif self.pam_orientation == False:  # 5prime = True 3prime=False
            if self.lsr == 0:
                return tseq
            else:
                return tseq[(len(tseq) - self.lsr):]
    # https://stackoverflow.com/questions/12555323/adding-new-column-to-existing-dataframe-in-python-pandas
    self.targets = deepcopy(self.targets)
    self.targets['seedseq'] = self.targets['target'].apply(_get_prox)
    self.targets['isseedduplicated'] = self.targets['seedseq'].duplicated()

get_control_seqs(seq_record_iter, configpath, length=20, n=10, num_threads=2)

Create random sequences with a specified GC probability and find seqs with the greatest distance to any sequence flanking a PAM site

Parameters:

Name	Type	Description	Default
seq_record_iter	SeqIO	An iterator of fastas	required
length	int	Length of the sequence, must match the index	20
n	int	Number of sequences to return	10
num_threads	int	Number of processor threads	2

Returns:

Type	Description
PandasDataFrame	A pandas dataframe with control sequence

Source code in guidemaker/core.py

def get_control_seqs(self, seq_record_iter: object, configpath, length: int = 20, n: int = 10,
                     num_threads: int = 2) -> pd.DataFrame:
    """
    Create random sequences with a specified GC probability and find seqs with the greatest
    distance to any sequence flanking a PAM site

    Args:
        seq_record_iter (Bio.SeqIO): An iterator of fastas
        length (int): Length of the sequence, must match the index
        n (int): Number of sequences to  return
        num_threads (int): Number of processor threads

    Returns:
        (PandasDataFrame): A pandas dataframe with control sequence
    """

    with open(configpath) as cf:
        config = yaml.safe_load(cf)

    MINIMUM_HMDIST = config['CONTROL']['MINIMUM_HMDIST']

    MAX_CONTROL_SEARCH_MULTIPLE = max(config['CONTROL']['CONTROL_SEARCH_MULTIPLE'])

    #  search_mult (int): search this times n sequences
    CONTROL_SEARCH_MULTIPLE = config['CONTROL']['CONTROL_SEARCH_MULTIPLE']

    # get GC percent
    totlen = 0
    gccnt = 0
    for record in seq_record_iter:
        gccnt += gc_fraction(record.seq) * len(record)
        totlen += len(record)
    gc = gccnt / (totlen)
    self.gc_percent = gc * 100
    self.genomesize = totlen / (1024 * 1024)

    minimum_hmdist = 0
    sm_count = 0
    search_mult = 0

    try:
        while minimum_hmdist < MINIMUM_HMDIST or search_mult == MAX_CONTROL_SEARCH_MULTIPLE:
            # generate random sequences
            seqs = []
            search_mult = CONTROL_SEARCH_MULTIPLE[sm_count]
            for i in range(n * search_mult):
                seqs.append("".join(np.random.choice(a=["G", "C", "A", "T"], size=length,
                                                     replace=True, p=[gc / 2, gc / 2, (1 - gc) / 2, (1 - gc) / 2])))
            # one hot encode sequences
            binseq = []
            charmap = {'A': '1 0 0 0', 'C': '0 1 0 0', 'G': '0 0 1 0', 'T': '0 0 0 1'}
            for seq in seqs:
                if self.targets['dtype'].iat[0] == "hamming":
                    charlist = [charmap[letter] for letter in seq]
                    binseq.append(" ".join(charlist))
                else: # leven
                    binseq.append(seq)

            rand_seqs = self.nmslib_index.knnQueryBatch(binseq, k=2, num_threads=num_threads)
            distlist = []
            for i in rand_seqs:
                distlist.append(i[1][0])
            zipped = list(zip(seqs, distlist))
            dist_seqs = sorted(zipped, reverse=True, key=lambda x: x[1])
            sort_seq = [item[0] for item in dist_seqs][0:n]

            #sort_dist
            if self.targets['dtype'].iat[0] == "hamming":
                sort_dist = [item[1] / 2 for item in dist_seqs][0:n]  ### ? does divide by 2 holds for leven???
            else:
                sort_dist = [item[1] for item in dist_seqs][0:n]  ### ? does divide by 2 holds for leven???

            minimum_hmdist = int(min(sort_dist))
            sm_count += 1
    except IndexError as e:
        raise e

    total_ncontrolsearched = search_mult * n
    self.ncontrolsearched = total_ncontrolsearched
    randomdf = pd.DataFrame(data={"Sequences": sort_seq, "Hamming distance": sort_dist})

    def create_name(seq):
        return "Cont-" + hashlib.md5(seq.encode()).hexdigest()
    randomdf['name'] = randomdf["Sequences"].apply(create_name)
    randomdf = randomdf[["name", "Sequences", "Hamming distance"]]
    randomdf.head()
    return (min(sort_dist),
            statistics.median(sort_dist),
            randomdf)

get_neighbors(configpath, num_threads=2)

Get nearest neighbors for sequences removing sequences that have neighbors less than the Hamming distance threshold. For the list of all targets calculate the (knum) nearest neighbors. filter out targets with close neighbors and Writes a dictionary to self.neighbors: self.neighbors[seq]{target: seq_obj, neighbors: {seqs:[s1, s1, ...], dist:[d1, d1,...]}}

Parameters:

Name	Type	Description	Default
configpath	str	Path to a parameter config file	required
num_threads	int	Number of threads	2

Returns:

Type	Description
None	None

Source code in guidemaker/core.py

def get_neighbors(self, configpath, num_threads=2) -> None:
    """
    Get nearest neighbors for sequences removing sequences that
    have neighbors less than the Hamming distance threshold.
    For the list of all targets calculate the (knum) nearest neighbors.
    filter out targets with close neighbors and
    Writes a dictionary to self.neighbors:
    self.neighbors[seq]{target: seq_obj, neighbors: {seqs:[s1, s1, ...], dist:[d1, d1,...]}}

    Args:
        configpath (str): Path to a parameter config file
        num_threads (int): Number of threads

    Returns:
        None
    """
    with open(configpath) as cf:
        config = yaml.safe_load(cf)

    ef = config['NMSLIB']['ef']

    # unique_targets = self.targets.loc[self.targets['isseedduplicated']
    #     == False]['target'].tolist()
    # For indexing we need to use all targets -- for checking off-targets. For searching neighbours remove seed duplicated and one wiht restriction site.
    unique_targets = self.targets.loc[(self.targets['isseedduplicated']==False) & (self.targets['hasrestrictionsite']==False)]['target'].tolist()
    if self.targets['dtype'].iat[0] == "hamming":
        unique_bintargets = self._one_hot_encode(unique_targets)  # search unique seed one
    else:
        unique_bintargets = unique_targets

    self.nmslib_index.setQueryTimeParams({'efSearch': ef})
    results_list = self.nmslib_index.knnQueryBatch(unique_bintargets,
                                           k=self.knum, num_threads=num_threads)
    neighbor_dict = {}
    target_lookup = getattr(self, 'notduplicated_targets', self.targets['target'].values)
    for i, entry in enumerate(results_list):
        queryseq = unique_targets[i]
        hitseqidx = entry[0].tolist()
        editdist = entry[1].tolist()
        if self.targets['dtype'].iat[0] == "hamming":
            # check that the closest sequence meets the min. dist. requirment. We multiply by 2 b/c each 
            # base is one hot encoded. e.g. 1000 vs 0100 = 2 differences
            if editdist[1] >= 2 * self.editdist:
                neighbors = {"seqs": [target_lookup[x] for x in hitseqidx],
                            "dist": [int(x / 2) for x in editdist]} 
                neighbor_dict[queryseq] = {"target": unique_targets[i],
                                        "neighbors": neighbors}
        else:
           if editdist[1] >= self.editdist: 
                neighbors = {"seqs": [target_lookup[x] for x in hitseqidx],
                            "dist": [int(x) for x in editdist]}
                neighbor_dict[queryseq] = {"target": unique_targets[i],
                                        "neighbors": neighbors}
    self.neighbors = neighbor_dict

extend_ambiguous_dna(seq, max_expansion=256)

Return list of all possible sequences given an ambiguous DNA input

Parameters:

Name	Type	Description	Default
seq	str	A DNA string	required
max_expansion	int	Maximum permitted combinations (default 256)	256

Return

List[str]: A list of DNA string with expanded ambiguous DNA values

Source code in guidemaker/core.py

def extend_ambiguous_dna(seq: str, max_expansion: int = 256) -> List[str]:
    """
    Return list of all possible sequences given an ambiguous DNA input

    Args:
        seq(str): A DNA string
        max_expansion(int): Maximum permitted combinations (default 256)

    Return:
        List[str]: A list of DNA string with expanded ambiguous DNA values
    """
    ambiguous_dna_values = {
    "A": "A",
    "C": "C",
    "G": "G",
    "T": "T",
    "M": "AC",
    "R": "AG",
    "W": "AT",
    "S": "CG",
    "Y": "CT",
    "K": "GT",
    "V": "ACG",
    "H": "ACT",
    "D": "AGT",
    "B": "CGT",
    "X": "GATC",
    "N": "GATC",
    }
    prod_count = 1
    for char in seq.upper():
        prod_count *= len(ambiguous_dna_values.get(char, char))
    if prod_count > max_expansion:
        raise ValueError(f"Restriction enzyme sequence '{seq}' produces {prod_count} combinations, exceeding maximum expansion limit of {max_expansion}.")

    extend_list = []
    for i in product(*[ambiguous_dna_values[j] for j in seq.upper()]):
        extend_list.append("".join(i))
    return extend_list

get_fastas(filelist, input_format='genbank', tempdir=None, max_decompressed_bytes=500 * 1024 * 1024)

Saves a Fasta and from 1 or more Genbank files (may be gzipped)

Parameters:

Name	Type	Description	Default
filelist	str	Genbank file to process	required
max_decompressed_bytes	int	Maximum decompressed size limit in bytes (default 500MB)	500 * 1024 * 1024

Returns:

Type	Description
	None

Source code in guidemaker/core.py

def get_fastas(filelist, input_format="genbank", tempdir=None, max_decompressed_bytes=500 * 1024 * 1024):
    """
    Saves a Fasta and from 1 or more Genbank files (may be gzipped)

    Args:
        filelist (str): Genbank file to process
        max_decompressed_bytes (int): Maximum decompressed size limit in bytes (default 500MB)

    Returns:
        None
    """
    try:
        fastpath = os.path.join(tempdir, "forward.fasta")
        bytes_written = 0
        with open(fastpath, "w") as f1:
            for file in filelist:
                if is_gzip(file):
                    f_in = gzip.open(file, 'rt')
                else:
                    f_in = open(file, 'r')
                with f_in as f:
                    records = SeqIO.parse(f, input_format)
                    for record in records:
                        record_str = f">{record.id}\n{str(record.seq).upper()}\n"
                        bytes_written += len(record_str.encode('utf-8'))
                        if bytes_written > max_decompressed_bytes:
                            raise ValueError(f"Decompressed genome file size exceeds maximum limit of {max_decompressed_bytes / (1024*1024):.0f}MB")
                        f1.write(record_str)
        return fastpath
    except Exception as e:
        logger.exception("An error occurred in the input file %s" % file)
        raise e

GuideMaker: The command line interface A command line Software to design gRNAs pools in non-model genomes and CRISPR-Cas systems

main(arglist=None)

Run The complete GuideMaker workflow.

Source code in guidemaker/cli.py

def main(arglist: list = None):
    """Run The complete GuideMaker workflow."""
    # Set up logging
    parser = myparser()
    args = parser.parse_args(arglist)
    logger = _logger_setup(args.log)
    parserval(args)



    try:
        with open(args.config) as cf:
            config = yaml.safe_load(cf)
    except:
        print("Could not parse the configuration file.")
        raise SystemExit(1)

    try:
        logger.info("Configuration data loaded from {}:".format(args.config))
        logger.info(config)
    except:
        print("Could not find config file, exiting.")
        raise SystemExit(1)

    try:

        if args.tempdir:
            if not os.path.exists(args.tempdir):
                logger.warning("Specified location for tempfile (%s) does not \
                                 exist, using default location." % args.tempdir)
                os.mkdir(args.tempdir)
                tempdir = args.tempdir
            else:
                tempdir = tempfile.mkdtemp()
        else:
            tempdir = tempfile.mkdtemp(prefix='guidemaker_', dir=args.tempdir)
        logger.info("Temp directory is: %s" % (tempdir))
        if args.genbank:
            logger.info("Writing fasta file from genbank file(s)")
            fastapath = guidemaker.get_fastas(args.genbank, input_format="genbank", tempdir=tempdir)
        elif args.fasta:
            fastapath = guidemaker.get_fastas(args.fasta, input_format="fasta", tempdir=tempdir)
        logger.info("Identifying PAM sites in the genome")
        pamobj = guidemaker.core.PamTarget(args.pamseq, args.pam_orientation, args.dtype)
        seq_record_iter = SeqIO.parse(fastapath, "fasta")
        pamtargets = pamobj.find_targets(
            seq_record_iter=seq_record_iter, target_len=args.guidelength)
        tl = guidemaker.core.TargetProcessor(
            targets=pamtargets, lsr=args.lsr, editdist=args.dist, knum=args.knum)
        lengthoftl = len(tl.targets)
        logger.info("Checking guides for restriction enzymes")
        tl.check_restriction_enzymes(restriction_enzyme_list=args.restriction_enzyme_list)
        logger.info("Number of guides removed after checking for restriction enzymes: %d",
                     (lengthoftl - len(tl.targets)))
        logger.info("Identifying guides that are unique near the PAM site")
        tl.find_unique_near_pam()
        logger.info("Number of guides with non unique seed sequence: %d",
                     (tl.targets.isseedduplicated.sum()))
        logger.info("Indexing all potential guide sites: %s. This is the longest step." %
                     len(list(set(tl.targets['target'].tolist()))))
        tl.create_index(num_threads=args.threads, configpath=args.config)
        logger.info(
            "Identifying guides that have a hamming distance <= %s to all other potential guides", str(args.dist))
        tl.get_neighbors(num_threads=args.threads, configpath=args.config)
        logger.info("Formatting guide target data")
        tf_df = tl.export_bed()
        if args.raw_output_only:
            tf_df.to_csv(os.path.join(args.outdir, "rawguides.csv.gz"), index=False, header=["Chromosome", "Start", "Stop","gRNA", "Strand"])
            logger.info("Raw guides options was selected Guidemaker, so has completed opperations")
            raise SystemExit(0)

        logger.info("Create GuideMaker Annotation object")
        if args.genbank:
            anno = guidemaker.core.Annotation(annotation_list=args.genbank, annotation_type="genbank",
                                              target_bed_df=tf_df)
        elif args.gff:
            anno = guidemaker.core.Annotation(annotation_list=args.gff, annotation_type="gff",
                                              target_bed_df=tf_df)
        logger.info("Identify genomic features")
        anno.get_annotation_features()
        logger.info("Total number of %s in the input genome: %d" % anno.locuslen())
        logger.info("Find genomic features closest the guides")
        anno._get_nearby_features()
        logger.info("Select guides that start between +%s and -%s of a feature start" %
                     (args.before, args.into))
        anno._filter_features(before_feat=args.before, after_feat=args.into)
        logger.info("Select description columns")
        anno._get_qualifiers(configpath=args.config)
        logger.info("Format the output")
        anno._format_guide_table(tl)
        prettydf = anno._filterlocus(args.attribute_key, args.filter_by_attribute)
        # prettydf = anno.filter_pretty_df
        if args.doench_efficiency_score:
            logger.info("Creating Efficiency Score based on Doench et al. 2016 - only for NGG PAM...")
            prettydf = guidemaker.core.get_doench_efficiency_score(df=prettydf, pam_orientation=args.pam_orientation, num_threads=args.threads)

        if args.cfd_score:
            logger.info("Calculating CFD score for assessing off-target activity of gRNAs")
            prettydf = guidemaker.core.cfd_score(df=prettydf)

        fd_zero = prettydf['Feature distance'].isin([0]).sum()
        logger.info("Number of Guides within a gene coordinates i.e. zero Feature distance: %d", fd_zero)
        if not os.path.exists(args.outdir):
            os.makedirs(args.outdir)
        csvpath = os.path.join(args.outdir, "targets.csv.gz")
        prettydf.to_csv(csvpath, index=False)
        if args.controls > 0:
            logger.info("Creating random control guides")
            contpath = os.path.join(args.outdir, "controls.csv.gz")
            seq_record_iter = SeqIO.parse(fastapath, "fasta")
            cmin, cmed, randomdf = tl.get_control_seqs(seq_record_iter,
                                                       configpath=args.config,
                                                       length=args.guidelength,
                                                       n=args.controls,
                                                       num_threads=args.threads)
            randomdf.to_csv(contpath)
            logger.info("Number of random control searched: %d" % tl.ncontrolsearched)
            logger.info("Created %i control guides with a minimum distance of %d and a median distance of %d" % (
                args.controls, cmin, cmed))
            logger.info("Percentage of GC content in the input genome: " +
                        "{:.2f}".format(tl.gc_percent))
            logger.info("Total length of the genome: " + "{:.1f} MB".format(tl.genomesize))


        logger.info("GuideMaker completed, results are at %s" % args.outdir)
        logger.info("PAM sequence: %s" % args.pamseq)
        logger.info("PAM orientation: %s" % args.pam_orientation)
        logger.info("Genome strand(s) searched: %s" % "both")
        logger.info("Total PAM sites considered: %d" % lengthoftl)
        logger.info("Guide RNA candidates found: %d" % len(prettydf))
    except Exception as e:
        logger.exception("GuideMaker terminated with errors. See the log file for details.")
        raise SystemExit(1)
    try:
        if args.plot:
            logger.info("Creating Plots...")
            guidemaker.core.GuideMakerPlot(prettydf=prettydf, outdir=args.outdir)
            logger.info("Plots saved at: %s" % args.outdir)
    except Exception as e:
        logger.exception(e)
        raise SystemExit(1)
    try:
        if not args.keeptemp:
            shutil.rmtree(tempdir)
    except UnboundLocalError as e:
        logger.exception(e)
        raise SystemExit(1)
    except AttributeError as e:
        logger.exception(e)
        raise SystemExit(1)

doench_predict.py. Simplified module to run the model 'V3_model_nopos' from Doench et al. 2016 for on-target scoring.

For use in Guidemaker https://guidemaker.org. Adam Rivers, Unites States Department of Agriculture, Agricultural Research Service

The core code, https://github.com/MicrosoftResearch/Azimuth, is in Python2 and does not run well given changes to packages. Miles Smith worked on porting to Python3 in this repo: https://github.com/milescsmith/Azimuth, including a new branch that used Poetry to build. The work is not complete.

This work is derivative of that BSD 3-clause, Modified licensed work. The key changes are: 1. Much of the code needed for tasks other thant prediction of the V3_model_nopos was removed. 2. The Calculation of NGGX features was re-written. A bug that prevented scaling to thousands guides efficiently. 3. the Pickle model and scikit-learn were replaced with an Onnx model ('V3_model_nopos.onnx"), and medadata file ("V3_model_nopos_options.json") and onnxruntime for better persistence, security, and performance.

Reference:

John G. Doench, Nicolo Fusi, Meagan Sullender, Mudra Hegde, Emma W. Vaimberg, Katherine F. Donovan, Ian Smith, Zuzana Tothova, Craig Wilen , Robert Orchard , Herbert W. Virgin, Jennifer Listgarten, David E. Root. Optimized sgRNA design to maximize activity and minimize off-target effects for genetic screens with CRISPR-Cas9. Nature Biotechnology Jan 2016, doi:10.1038/nbt.3437.

concatenate_feature_sets(feature_sets, keys=None)

Combine features

Given a dictionary of sets of features, each in a pd.DataFrame, concatenate them together to form one big np.array, and get the dimension of each set

Parameters:

Name	Type	Description	Default
feature_sets	dict	a Dict of feature sets as pandas DataFrames	required

Returns:

Type	Description
tuple	(tuple: inputs(numpy.ndarray), dim (tuple)

Source code in guidemaker/doench_predict.py

def concatenate_feature_sets(feature_sets: dict, keys: List[str] = None) -> tuple:
    """ Combine features

    Given a dictionary of sets of features, each in a pd.DataFrame,
    concatenate them together to form one big np.array, and get the dimension
    of each set

    Args:
        feature_sets (dict): a Dict of feature sets as pandas DataFrames

    Returns:
        (tuple: inputs(numpy.ndarray), dim (tuple)
    """
    if feature_sets == {}:
        raise AssertionError("no feature sets present")
    if keys is None:
        keys = list(feature_sets.keys())
    feature_1 = feature_sets[keys[0]].shape[0]
    for assemblage in feature_sets:
        feature_2 = feature_sets[assemblage].shape[0]
        if feature_1 != feature_2:
            raise AssertionError(
                f"not same # items for features {keys[0]} and {assemblage}"
            )
    num_sets = feature_sets[keys[0]].shape[0]
    inputs = np.zeros((num_sets, 0))
    feature_names = []
    dim = {}
    dimsum = 0
    for assemblage in keys:
        inputs_set = feature_sets[assemblage].values
        dim[assemblage] = inputs_set.shape[1]
        dimsum = dimsum + dim[assemblage]
        inputs = np.hstack((inputs, inputs_set))
        feature_names.extend(feature_sets[assemblage].columns.tolist())
    return inputs, dim, dimsum, feature_names

predict(seq, model_file=MODEL, model_metadata=MODEL_META, pam_audit=True, length_audit=False, num_threads=1)

Predicts regression scores from sequences.

Parameters:

Name	Type	Description	Default
seq (numpy.ndarray) numpy array of 30 nt sequences with 25 nt of guide, NGG pam in 25		27 and the following 2 nts.	required
model_file	str	file path of the onnx Boosted Gradient Regressor model file without position data	MODEL
model_metadata	str	file path of the json model parameters metadata file.	MODEL_META
pam_audit	bool	check PAM of each sequence.	True
length_audit(bool)		check length of each sequence.	required

Returns:

Type	Description
array	An array with regression values.

Source code in guidemaker/doench_predict.py

def predict(
    seq: np.ndarray,
    model_file: Optional[Path] = MODEL,
    model_metadata: Optional[Path] = MODEL_META,
    pam_audit: bool = True,
    length_audit: bool = False,
    num_threads: int = 1
) -> np.array:
    """Predicts regression scores from sequences.

    Args:
        seq (numpy.ndarray) numpy array of 30 nt sequences with 25 nt of guide, NGG pam in 25:27 and the following 2 nts.
        model_file (str): file path of the onnx Boosted Gradient Regressor model file without position data
        model_metadata (str): file path of the json model parameters metadata file.
        pam_audit (bool): check PAM of each sequence.
        length_audit(bool) : check length of each sequence.

    Returns:
        (numpy.array): An array with regression values.

     """
    if not isinstance(seq, np.ndarray):
        raise AssertionError("Please ensure seq is a numpy array")
    if len(seq[0]) <= 0:
        raise AssertionError("Make sure that seq is not empty")
    if not isinstance(seq[0], str):
        raise AssertionError(
            f"Please ensure input sequences are in string format, i.e. 'AGAG' "
            f"rather than ['A' 'G' 'A' 'G'] or alternate representations"
        )
    # Open onnx runtime session
    sess = rt.InferenceSession(model_file)
    with open(model_metadata, "r") as f:
        learn_options = json.load(f)
    x_df = pd.DataFrame(
        columns=["30mer", "Strand"],
        data=list(zip(seq, ["NA" for x in range(len(seq))])),
    )

    feature_sets = parallel_featurize_data(
        x_df,
        learn_options,
        pam_audit=pam_audit,
        length_audit=length_audit,
        num_threads=num_threads
    )
    inputs, *_ = concatenate_feature_sets(feature_sets)
    preds = sess.run(None, {'float_input': inputs.astype(np.float32)})[0]
    return preds

doench_featurization.py. Simplified feature extraction to run the model 'V3_model_nopos' from Doench et al. 2016 for on-target scoring.

For use in Guidemaker https://guidemaker.org. Adam Rivers, Unites States Department of Agriculture, Agricultural Research Service

Core code https://github.com/MicrosoftResearch/Azimuth is in Python2 and does not run well given changes to packages. Miles Smith worked on porting to Python3 in this repo: https://github.com/milescsmith/Azimuth. including a new branch that used Poetry to build. The work is not complete.

This work is derivitive of that BSD 3-clause licensed work. The key changes are: 1. Much of the code needed for tasks other than prediction of the V3_model_nopos was removed. 2. The Calculation of NGGX features was re-written to fix a bug that prevented scaling to thousands guides efficiently. 3. the Pickle model and scikit-learn were replaced with an Onnx model and onnxruntime for better persistance, security, and performance.

Reference:

John G. Doench, Nicolo Fusi, Meagan Sullender, Mudra Hegde, Emma W. Vaimberg, Katherine F. Donovan, Ian Smith, Zuzana Tothova, Craig Wilen , Robert Orchard , Herbert W. Virgin, Jennifer Listgarten, David E. Root. Optimized sgRNA design to maximize activity and minimize off-target effects for genetic screens with CRISPR-Cas9. Nature Biotechnology Jan 2016, doi:10.1038/nbt.3437.

Tm_feature(data, pam_audit=True, learn_options=None)

assuming '30-mer'is a key get melting temperature features from: 0-the 30-mer ("global Tm") 1-the Tm (melting temperature) of the DNA:RNA hybrid from positions 16 - 20 of the sgRNA, i.e. the 5nts immediately proximal of the NGG PAM 2-the Tm of the DNA:RNA hybrid from position 8 - 15 (i.e. 8 nt) 3-the Tm of the DNA:RNA hybrid from position 3 - 7 (i.e. 5 nt)

Source code in guidemaker/doench_featurization.py

def Tm_feature(data, pam_audit=True, learn_options=None):
    """
    assuming '30-mer'is a key
    get melting temperature features from:
        0-the 30-mer ("global Tm")
        1-the Tm (melting temperature) of the DNA:RNA hybrid from positions 16 - 20 of the sgRNA,
        i.e. the 5nts immediately proximal of the NGG PAM
        2-the Tm of the DNA:RNA hybrid from position 8 - 15 (i.e. 8 nt)
        3-the Tm of the DNA:RNA hybrid from position 3 - 7  (i.e. 5 nt)
    """

    if learn_options is None or "Tm segments" not in learn_options:
        segments = [(19, 24), (11, 19), (6, 11)]
    else:
        segments = learn_options["Tm segments"]

    sequence = data["30mer"].values
    featarray = np.ones((sequence.shape[0], 4))

    for i, seq in enumerate(sequence):
        if pam_audit and seq[25:27] != "GG":
            raise Exception(f"expected GG but found {seq[25:27]}")
        rna = False
        featarray[i, 0] = Tm.Tm_NN(seq, nn_table=Tm.RNA_NN2)  # 30mer Tm
        featarray[i, 1] = Tm.Tm_NN(
            seq[segments[0][0]: segments[0][1]], nn_table=Tm.RNA_NN2
        )  # 5nts immediately proximal of the NGG PAM
        featarray[i, 2] = Tm.Tm_NN(
            seq[segments[1][0]: segments[1][1]], nn_table=Tm.RNA_NN2
        )  # 8-mer
        featarray[i, 3] = Tm.Tm_NN(
            seq[segments[2][0]: segments[2][1]], nn_table=Tm.RNA_NN2
        )  # 5-mer

    feat = pd.DataFrame(
        featarray,
        index=data.index,
        columns=[
            f"Tm global_{rna}",
            f"5mer_end_{rna}",
            f"8mer_middle_{rna}",
            f"5mer_start_{rna}",
        ],
    )

    return feat

check_feature_set(feature_sets)

Ensure the number of features is the same in each feature set

Parameters:

Name	Type	Description	Default
feature_sets	dict	the feature set dictionary	required

Returns: None

Source code in guidemaker/doench_featurization.py

def check_feature_set(feature_sets: dict) -> None:
    """Ensure the number of features is the same in each feature set

    Args:
            feature_sets (dict): the feature set dictionary
    Returns:
        None
    """
    if feature_sets == {}:
        raise AssertionError("no feature sets present")

    num = None
    for ft in feature_sets:
        num2 = feature_sets[ft].shape[0]
        if num is None:
            num = num2
        else:
            if num < 1:
                raise AssertionError("should be at least one individual")
            if num != num2:
                raise AssertionError(
                    "Number of individuals do not match up across feature sets"
                )

    for item in feature_sets:
        if np.any(np.isnan(feature_sets[item])):
            raise Exception(f"found Nan in set {item}")

countGC(s, length_audit=True)

GC content for only the 20mer, as per the Doench paper/code

Source code in guidemaker/doench_featurization.py

def countGC(s: str, length_audit: bool=True) -> int:
    """
    GC content for only the 20mer, as per the Doench paper/code
    """
    if length_audit:
        if len(s) != 30:
            raise AssertionError("seems to assume 30mer")
    #return len(s[4:24].replace("A", "").replace("T", ""))
    return s[4:24].count("G") + s[4:24].count("C")

featurize_data(data, learn_options, pam_audit=True, length_audit=True)

Creates a dictionary of feature data

Parameters:

Name	Type	Description	Default
data pd.DataFrame		of 30-mer sequences in column 1 and strand in column 2	required
learn_options	dict	dict of model training parameters	required
pam_audit	bool	should a check of GG at position 25:27 be performed?	True
length_audit	bool	should sequence length be checked?	True

Returns:

Type	Description
dict	Returns a dict containing pandas dataframs of features

Source code in guidemaker/doench_featurization.py

def featurize_data(data: pd.DataFrame, learn_options: dict, pam_audit: bool=True, length_audit :bool=True) -> dict:
    """Creates a dictionary of feature data

    Args:
        data pd.DataFrame: of 30-mer sequences in column 1 and strand in column 2
        learn_options (dict): dict of model training parameters
        pam_audit (bool): should a check of GG  at position 25:27 be performed?
        length_audit (bool): should sequence length be checked?

    Returns:
        (dict): Returns a dict containing pandas dataframs of features
    """

    logger.info("Creating features for Doench et al. 2016 score prediction")

    # Coerce and explicitly restore the column name expected by the rest of the script
    data = pd.DataFrame(data)
    data.columns = ["30mer"] + list(data.columns[1:])

    if np.any(data["30mer"].str.len() != 30):
        raise AssertionError(f"Sequences should be 30 nt long")



    feature_sets = {}

    if learn_options["nuc_features"]:
        # spectrum kernels (position-independent) and weighted degree kernels (position-dependent)
        logger.info("Creating nucleotide features")
        feature_sets["_nuc_pd_Order1"], feature_sets["_nuc_pi_Order1"], feature_sets["_nuc_pd_Order2"], feature_sets["_nuc_pi_Order2"] = get_nuc_features(data)

    logger.info("Verifying nucleotide features")
    check_feature_set(feature_sets)

    if learn_options["gc_features"]:
        logger.info("Creating GC features")
        gc_above_10, gc_below_10, gc_count = gc_features(data, length_audit)
        feature_sets["gc_above_10"] = pd.DataFrame(gc_above_10)
        feature_sets["gc_below_10"] = pd.DataFrame(gc_below_10)
        feature_sets["gc_count"] = pd.DataFrame(gc_count)
        logger.info("gc features complete")

    if learn_options["include_NGGX_interaction"]:
        logger.info("Creating ggx features")
        feature_sets["NGGX"] = nggx_interaction_feature(data, pam_audit)

    if learn_options["include_Tm"]:
        logger.info("Creating Tm features")
        feature_sets["Tm"] = Tm_feature(data, pam_audit, learn_options=None)


    check_feature_set(feature_sets)
    logger.info("final feature check complete")

    return feature_sets

get_nuc_features(data)

Create first and second order nucleotide features

Parameters:

Name	Type	Description	Default
data pd.DataFrame		of 30-mer sequences in column 1 and strand in column 2	required

Returns:

Type	Description
tuple	Returns a tuple pwith 4 Pandas dataframes

Source code in guidemaker/doench_featurization.py

def get_nuc_features(data: pd.DataFrame) -> tuple[pd.DataFrame, pd.DataFrame, pd.DataFrame, pd.DataFrame]:
    """ Create first and second order nucleotide features

        Args:
            data pd.DataFrame:      of 30-mer sequences in column 1 and strand in column 2

        Returns:
            (tuple): Returns a tuple pwith 4 Pandas dataframes

    """
    seqlen = 30
    # create first header
    nuc_pi_Order1_header = [x[0] for x  in product('ATCG', repeat=1)]

    # create second header
    nuc_pd_Order1_header = []
    for i in range(seqlen):
        nuc_pd_Order1_header.extend([x + "_" + str(i) for x in  nuc_pi_Order1_header ])

    # create third header
    nuc_pi_Order2_header = [ "".join(x) for x in product('ATCG', repeat=2)]

    # create forth header
    nuc_pd_Order2_header = []
    for i in range(seqlen - 1):
        nuc_pd_Order2_header.extend([x + "_" + str(i) for x in  nuc_pi_Order2_header ])

    # Create lists for holding features
    nuc_pd_Order1_list = []
    nuc_pi_Order1_list = []
    nuc_pd_Order2_list = []
    nuc_pi_Order2_list = []

    # Create identity matricies and lookups for one hot encoding
    o1_id = np.eye(4)
    o2_id = np.eye(16)
    o1_lookup = {x : i for i, x in enumerate(nuc_pi_Order1_header)}
    o2_lookup = {x : i for i, x in enumerate(nuc_pi_Order2_header)}

    # not feasable ot use pandas.get_dummies for this
    def one_hot(seq, idmat, lookup):
        """One hot endode  an iterable matching the items in a lookup

        Args:
            seq (str): a 30mer sequence string
            idmat (np.array): a Numpy identity matrix
            lookup (dict): a dictionary matchin the substring to the row in idmat

        Returns:
            pd.Series: one hot encoding
        """
        featurevect = []
        num_cols = idmat.shape[1]
        for let in seq:
            if let in lookup:
                pos = lookup[let]
                featurevect.extend(list(idmat[pos, :]))
            else:
                featurevect.extend([0.0] * num_cols)
        return pd.Series(featurevect)


    def sliding_window(iterable, n):
        """Create a generator of substrings

            Args:
                iterable (iterable): an itterable object like a string or list
                n (int): the size of the window or kmer

            Returns:
                a genrator of substrings
        """
        # sliding_window('ABCDEFG', 4) -> ABCD BCDE CDEF DEFG
        it = iter(iterable)
        window = deque(islice(it, n), maxlen=n)
        if len(window) == n:
            yield tuple(window)
        for x in it:
            window.append(x)
            yield tuple(window)

    for seq in data["30mer"]:
        # add order 1 frequency features
        pi1dict = dict.fromkeys(nuc_pi_Order1_header, 0)
        for let in seq:
            if let in pi1dict:
                pi1dict[let] += 1
        nuc_pi_Order1_list.append(pi1dict)
        # add order 1 positon features
        nuc_pd_Order1_list.append(one_hot(seq, o1_id, o1_lookup))
        # create list of 2mers
        seq_2mers =  [ "".join(x) for x in list(sliding_window(seq, 2))]
        # add order two frequency features
        pi2dict = dict.fromkeys(nuc_pi_Order2_header, 0)
        for let in seq_2mers:
            if let in pi2dict:
                pi2dict[let] += 1
        nuc_pi_Order2_list.append(pi2dict)
        # add order 2 positon features
        nuc_pd_Order2_list.append(one_hot(seq_2mers, o2_id, o2_lookup))

    # Create DataFrames
    nuc_pd_Order1 = pd.DataFrame(data=nuc_pd_Order1_list)
    nuc_pd_Order1.columns = nuc_pd_Order1_header
    nuc_pi_Order1 = pd.DataFrame(data=nuc_pi_Order1_list, columns=nuc_pi_Order1_header)
    nuc_pd_Order2 = pd.DataFrame(data=nuc_pd_Order2_list)
    nuc_pd_Order2.columns = nuc_pd_Order2_header
    nuc_pi_Order2 = pd.DataFrame(data=nuc_pi_Order2_list, columns=nuc_pi_Order2_header)

    return nuc_pd_Order1, nuc_pi_Order1, nuc_pd_Order2, nuc_pi_Order2

nggx_interaction_feature(data, pam_audit=True)

One hot encode the sequence of NX aroung pam site NGGX

Parameters:

Name	Type	Description	Default
data	DataFrame	A dataframe of 30-mer and strand (filled with NA)	required
pam_audit	bool	should check of GG at position 25:27 be performed?	True

Returns:

Type	Description
DataFrame	A dataframe with 16 columns containing one hoe encoding of NX data

Source code in guidemaker/doench_featurization.py

def nggx_interaction_feature(data: pd.DataFrame, pam_audit: bool=True) -> pd.DataFrame:
    """ One hot encode the sequence of NX aroung pam site NGGX

    Args:
        data (pandas.DataFrame): A dataframe of 30-mer and strand (filled with NA)
        pam_audit (bool): should check of GG  at position 25:27 be performed?

    Returns:
        (pandas.DataFrame):  A dataframe with 16 columns containing  one hoe encoding of NX data
    """
    # function completly replaced from Doench et al. (2016). The old function had time complexity of O^2
    sequence = data["30mer"].values
    nxcombos = []
    for i in product('ACGT', repeat=2):
        nxcombos.append("".join(list(i)))
    nxlist = []
    # check that GG is where we think
    for seq in sequence:
        if pam_audit and seq[25:27] != "GG":
            raise Exception(f"expected GG but found {seq[25 :27]}")
        nxlist.append(seq[24] + seq[27])
    nxs = pd.Series(nxlist)
    feat_nx = pd.get_dummies(nxs, columns=nxcombos)
    # Add any missing columns
    for i, header in enumerate(nxcombos):
        if header not in feat_nx.columns:
            feat_nx.insert(i, header, 0)
    feat_nx.columns = ["NGGX" + i + "_0" for i in nxcombos]
    feat_nx = feat_nx.astype('float64')
    return feat_nx

normalize_features(data, axis)

input: pd.DataFrame of dtype=np.float64 array, of dimensions mean-center, and unit variance each feature

Source code in guidemaker/doench_featurization.py

def normalize_features(data, axis):
    """
    input: pd.DataFrame of dtype=np.float64 array, of dimensions
    mean-center, and unit variance each feature
    """
    data -= data.mean(axis)
    data /= data.std(axis)
    # remove rows with NaNs
    data = data.dropna(1)
    if np.any(np.isnan(data.values)):
        raise Exception("found NaN in normalized features")
    return data

organism_feature(data)

Human vs. mouse

Source code in guidemaker/doench_featurization.py

def organism_feature(data):
    """
    Human vs. mouse
    """
    organism = np.array(data["Organism"].values)
    feat = pd.DataFrame(pd.DataFrame(organism))
    import pdb

    pdb.set_trace()
    return feat

parallel_featurize_data(data, learn_options, pam_audit=True, length_audit=True, num_threads=1)

Use multprocessing to divide up the creation of ML features for Doench scoring Creates a dictionary of feature data

Parameters:

Name	Type	Description	Default
data pd.DataFrame		of 30-mer sequences in column 1 and strand in column 2	required
learn_options	dict	dict of model training parameters	required
pam_audit	bool	should a check of GG at position 25:27 be performed?	True
length_audit	bool	should sequence length be checked?	True

Returns:

Type	Description
dict	Returns a dict containing pandas dataframs of features

Source code in guidemaker/doench_featurization.py

def parallel_featurize_data(data: pd.DataFrame, learn_options: dict, pam_audit: bool=True, length_audit: bool=True, num_threads: int=1) -> dict:
    """ Use multprocessing to divide up the creation of ML features for Doench scoring
        Creates a dictionary of feature data

        Args:
            data pd.DataFrame:      of 30-mer sequences in column 1 and strand in column 2
            learn_options (dict):   dict of model training parameters
            pam_audit (bool):       should a check of GG  at position 25:27 be performed?
            length_audit (bool):    should sequence length be checked?

        Returns:
            (dict): Returns a dict containing pandas dataframs of features

    """
    if num_threads > 1:
        dflist = np.array_split(data, num_threads)
        partial_fd = partial(featurize_data,learn_options=learn_options, pam_audit=pam_audit, length_audit=length_audit )
        with Pool(processes=num_threads) as pool:
            result = pool.map(partial_fd, dflist)
        featdict = dict.fromkeys(result[0].keys())
        for featkey in featdict.keys():
            tempdflist = []
            for d1 in result:
                tempdflist.append(d1[featkey])
            featdict[featkey] = pd.concat(tempdflist)
        return featdict
    else:
        return featurize_data(data=data, learn_options=learn_options, pam_audit=pam_audit, length_audit=length_audit)

cfd_score_calculator.py This is a modified version of the CDF score calculator in Doench et al. (2016) for use in Guidemaker (https://guidemaker.org) Adam Rivers, USDA Agricultural Research Service

We score only the CFD for off targets with a NGG site, we do not collect these non-matching PAM off targets Guidemaker. For this reason we omit the PAM scoring portion of CFD. For that reason we omit the pam scoring part of the doench et al. (2016) script. Results are identical for all off-targets that are scored.

Very few off trargets with non-pam matching sites would interact with targets in a small geneome (The highest scoring non-Pam,NGT, has a score of 0.3). Additionally we require all our guides have a distance of at least 2 by default so any off targets would have a score below the 0.2 threshold most people use.

We also modified the script to score pam sites longer than 20 by ignoring the 5' end past 20 and for shorter pam's by only scoring the sites present.

calc_cfd(wt, off, mm_scores=None)

Calculate the CFD score using precalculated weights

Parameters:

Name	Type	Description	Default
wt	str	wild-type gRNA sequence, excluding the PAM Cas9 site	required
off	str	off target sequence, excluding the PAM Cas9 site	required

Returns:

Type	Description
float	CDF score of the pair

Source code in guidemaker/cfd_score_calculator.py

def calc_cfd(wt: str, off: str, mm_scores=None) -> float:
    """Calculate the CFD score using precalculated weights

    Args:
        wt: wild-type gRNA sequence, excluding the PAM Cas9 site
        off: off target sequence, excluding the PAM Cas9 site

    Returns:
        (float): CDF score of the pair

    """
    guidelen = check_len(wt, off)
    if mm_scores is None:
        mm_scores, _ = get_mm_pam_scores()
    score = 1.
    off = off.upper().replace('T', 'U')
    wt = wt.upper().replace('T', 'U')
    s_list = list(off)
    wt_list = list(wt)
    basecomp = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A', 'U': 'A'}
    for i, sl in enumerate(s_list):
        if (guidelen - 20 - i) <= 0:
            if wt_list[i] != sl:
                key = 'r' + wt_list[i] + ':d' + basecomp[sl] + ',' + str(20 + i + 1 - guidelen)
                score *= mm_scores[key]
    return score

check_len(wt, off)

Verify the lengths of guide and off target match returning the length

Parameters:

Name	Type	Description	Default
wt	str	the guide type guide sequence	required
off	str	the off target sequence	required

Returns:

Type	Description
int	the length of the data

Source code in guidemaker/cfd_score_calculator.py

def check_len(wt: str, off: str) -> int:
    """Verify the lengths of guide and off target match returning the length

    Args:
        wt: the guide type guide sequence
        off: the off target sequence

    Returns:
        (int): the length of the data

    """
    wtl = len(wt)
    offl = len(off)
    assert (wtl == offl), "The lengths wt and off differ: wt = {}, off = {}".format(str(wtl), str(offl))
    return wtl

get_mm_pam_scores()

load json file of mismatch scores and PAM scores

Returns:

Type	Description
tuple	dict of mismatch scores, dict of pam scores

Source code in guidemaker/cfd_score_calculator.py

def get_mm_pam_scores() -> Tuple[Dict, Dict]:
    """load json file of mismatch scores and PAM scores

    Returns:
        (tuple):dict of mismatch scores, dict of pam scores

    """
    try:
        with open(MODEL_META) as dat:
            scores = json.load(dat)
        mm_s = scores['mm']
        pam_s = scores['pam']
        return mm_s, pam_s
    except (FileNotFoundError, IOError):
        raise Exception("Could not find file with reference mismatch scores and PAM scores")

Web Application

Run web App.

Source code in guidemaker/app.py

def main(arglist: list = None):
    """Run web App."""
    header = "GuideMaker"
    subheader = "Software to design CRISPR-Cas guide RNA pools in non-model genomes 🦠 🧬"
    st.markdown(f'<strong style="font-family:Hoefler Text;font-size: 36px;color: #0021A5">{header}</strong>',
                unsafe_allow_html=True)
    st.markdown(
        f'<strong style="font-family:Hoefler Text;font-size: 18px;color: #FA4616">{subheader}</strong>', unsafe_allow_html=True)

    st.markdown("---")

    sessionID = str(uuid.uuid1())
    # st.write(sessionID)
    logfilename = sessionID + "_log.txt"

    # Create a downloads directory within the streamlit static asset directory
    # and write output files to it. Binary downloader as html file  has limited
    # file size for conversion. So, we need to write it to local folder.
    STREAMLIT_STATIC_PATH = pathlib.Path(st.__path__[0]) / 'static'
    DOWNLOADS_PATH = (STREAMLIT_STATIC_PATH / "downloads")
    if not DOWNLOADS_PATH.is_dir():
        DOWNLOADS_PATH.mkdir()

    # Define input parameters and widgets

    multiple_files_gbk = st.sidebar.file_uploader("Upload one or more Genome file [ .gbk, .gbk.gz]", type=[".gbk", ".gz",".gbff"], accept_multiple_files=True)
    genome = list( map(lambda x: x.getvalue(), multiple_files_gbk))

    multiple_files_fasta = st.sidebar.file_uploader("Upload one or more fasta file [ .fasta, .fasta.gz]", type=[".fasta", ".gz",".fna"], accept_multiple_files=True)
    fasta = list( map(lambda x: x.getvalue(), multiple_files_fasta))


    multiple_files_gff= st.sidebar.file_uploader("Upload gff/gtf file if you are using fasta [ .gff, .gtf]", type=[".gff", ".gtf"], accept_multiple_files=True)
    gff = list( map(lambda x: x.getvalue(), multiple_files_gff))

    DemoGenome = st.sidebar.selectbox("OR Use Demo GBK",['Carsonella_ruddii.gbk.gz','Pseudomonas_aeruginosa.gbk.gz'])
    DemoGenomepath = os.path.join(DATA_DIR,DemoGenome)
    demo = open(DemoGenomepath,"rb")
    #st.write("You selected this option ",xx)



    pam = st.sidebar.text_input("Input PAM Motif [ E.g. NGG ] ", "NGG")
    restriction_enzyme_list = st_tags_sidebar(label = 'Restriction Enzymes[e.g. NGRT]:',
                                                                  text  = 'Enter to add more', 
                                                                  value = ['NGRT'])
    #restriction_enzyme_list= st.sidebar.text_input("Restriction Enzymes list [ E.g. NGRT ] ", "NGRT")
    pam_orientation = st.sidebar.selectbox(
        "PAM Orientation [ Options: 3prime, 5prime ]", ("3prime", "5prime"))
    guidelength = st.sidebar.number_input('Guidelength [ Options: 10 - 27 ]', 10, 27, value=20)
    lsr = st.sidebar.number_input('Length of seed region[ Options: 0 - 27 ]', 0, 27, value=10)
    #dtype = st.sidebar.selectbox(
    #  "Type of Edit Distance [ Options: hamming, leven ]", ("hamming", "leven"))
    dist = st.sidebar.number_input('Edit Distance [Options: 0 - 5 ]', 0, 5, value=2)
    before = st.sidebar.number_input('Before [Options: 1 - 500 ]', 1, 500, value=100, step=50)
    into = st.sidebar.number_input('Into [Options: 1 - 500 ]', 1, 500, value=200, step=50)
    knum = st.sidebar.number_input('Similar Guides[Options: 2 - 20 ]', 2, 20, value=3)
    controls = st.sidebar.number_input('Control RNAs', 1, 1000, value=10, step=100)
    #threads = st.sidebar.number_input('Threads [ Options: 2, 4, 6, 8]', 2, 8, step=2)


    # Validate restriction enzyme tags to prevent CLI option injection (F-01)
    clean_enzyme_list = []
    for tag in restriction_enzyme_list:
        tag_str = str(tag).strip()
        if tag_str.startswith('-') or any(c not in 'ACGTMRWSYKVHDBXNacgtmrwsykvhdbxn' for c in tag_str):
            st.error(f"Invalid restriction enzyme tag '{tag_str}'. Tags cannot start with '-' or contain non-IUPAC DNA characters.")
            return
        clean_enzyme_list.append(tag_str.upper())

    gbk_filename = f"input_{sessionID}.gbk"
    fasta_filename = f"input_{sessionID}.fasta"
    gff_filename = f"input_{sessionID}.gff"

    scriptorun = None
    input_context = None

    if genome:
        input_context = genome_connect(genome, gbk_filename)
        input_args = ["-i", gbk_filename]
    elif fasta and gff:
        @contextmanager
        def combined_fasta_gff():
            with fasta_connect(fasta, fasta_filename) as f_conn:
                with gff_connect(gff, gff_filename) as g_conn:
                    yield (f_conn, g_conn)
        input_context = combined_fasta_gff()
        input_args = ["-f", fasta_filename, "-g", gff_filename]
    elif demo:
        input_context = genome_connect(demo, gbk_filename)
        input_args = ["-i", gbk_filename]

    if input_context is not None:
        args = ["guidemaker"] + input_args + [
            "-p", pam,
            "--guidelength", str(guidelength),
            "--pam_orientation", pam_orientation,
            "--lsr", str(lsr),
            "--dtype", "hamming",
            "--dist", str(dist),
            "--outdir", sessionID,
            "--log", logfilename,
            "--into", str(into),
            "--before", str(before),
            "--knum", str(knum),
            "--controls", str(controls),
            "--threads", "2",
            "--cfd_score",
            "--doench_efficiency_score",
            "--restriction_enzyme_list"
        ] + clean_enzyme_list
        scriptorun = args

    if st.sidebar.button("SUBMIT"):
        if scriptorun and input_context:
            with input_context:
                run_command(scriptorun)

    if os.path.exists(sessionID):
        source = pd.read_csv(os.path.join("./", sessionID, 'targets.csv.gz'), low_memory=False)

        accession_list = list(set(source['Accession']))
        for accession in accession_list:
            accession_df = source[source["Accession"] == accession]
            accession_info = f"**Accession:** {accession}"
            st.markdown(accession_info)
            st.write(guidemakerplot(accession_df))

        # F-02: Isolated downloads path per session
        session_downloads_path = DOWNLOADS_PATH / sessionID
        session_downloads_path.mkdir(parents=True, exist_ok=True)

        targets_out_path = session_downloads_path / "targets.csv.gz"
        controls_out_path = session_downloads_path / "controls.csv.gz"

        # Targets
        target_tab = f"✅ [Target Data](downloads/{sessionID}/targets.csv.gz)"
        targets = pd.read_csv(os.path.join("./", sessionID, 'targets.csv.gz'), low_memory=False)
        targets.to_csv(str(targets_out_path), index=False)

        # Controls
        control_tab = f"✅ [Control Data](downloads/{sessionID}/controls.csv.gz)"
        controls = pd.read_csv(os.path.join("./", sessionID, 'controls.csv.gz'), low_memory=False)
        controls.to_csv(str(controls_out_path), index=False)

        # logs
        with st.expander("Results"):
            st.write(target_tab)
            st.write(control_tab)
            if os.path.exists(logfilename):
                st.write(get_binary_file_downloader_html(
                    logfilename, '✅ Log File'), unsafe_allow_html=True)

    # Parameters Dictionary
    image = Image.open(guidemaker.APP_PARAMETER_IMG)
    optionals = st.expander("Parameter Dictionary", False)
    optionals.image(image, caption='GuideMaker Parameters', width='stretch')

    with st.expander("Designing Experiments with GuideMaker Results"):
        intro_markdown = read_markdown_file(guidemaker.APP_EXPERIMENT_FILE)
        st.markdown(intro_markdown, unsafe_allow_html=True)

    st.markdown("""
    ##### API documentation 📖

    API documentation for the module can be found [here](https://github.com/USDA-ARS-GBRU/GuideMaker)


    ##### License information ©️

    *Guidemaker was created by the United States Department of Agriculture - Agricultural Research Service (USDA-ARS). As a work of the United States Government this software is available under the CC0 1.0 Universal Public Domain Dedication (CC0 1.0)*

    """)


    # Cleanup session-specific files
    try:
        shutil.rmtree(sessionID, ignore_errors=True)
        if os.path.exists(logfilename):
            os.remove(logfilename)
        for fn in [gbk_filename, fasta_filename, gff_filename]:
            if os.path.exists(fn):
                os.remove(fn)
    except Exception:
        pass