Graph module

Provides main class and helpers classes and functions to handle pangenome graph and related data structures

Imports

Path conversion utility functions

These functions are utility, but because they are directly related to graph structures, they are places in this module and not in Util module

calcNodeLengths

 calcNodeLengths (graph)

Simple function that calculates node lengths (in visual columns).

If it is nucleotide graph, it will actually calculate a number of nucleotides in each node, but if it is non-nucleotide graph, then it will return 1 for each node.

/home/pigrenok/.pyenv/versions/3.10.9/envs/pygengraph/lib/python3.10/site-packages/fastcore/docscrape.py:225: UserWarning: Unknown section Return
  else: warn(msg)
/home/pigrenok/.pyenv/versions/3.10.9/envs/pygengraph/lib/python3.10/site-packages/fastcore/docscrape.py:225: UserWarning: Unknown section Example
  else: warn(msg)

initialPathAnalysis

 initialPathAnalysis (graph, nodeLengths)

This function creates auxiliary data structures to make it easier to work with paths and their relationships with nodes.

getNodesStructurePathNodeInversionRate

 getNodesStructurePathNodeInversionRate (pathNodeArray, pathDirArray,
                                         pathLengths,
                                         inversionThreshold=0.5)

Generate a dict of dicts which stores information about inversion rate of each node for each path.

convertPathsToGraph

 convertPathsToGraph (fullPath, doSorting=False, v2=False)

getPathNodeInversionRate

 getPathNodeInversionRate (pathNodeArray, pathDirArray, pathLengths,
                           inversionThreshold=0.5)

deprecated Generate a dict of dicts which stores information about inversion rate of each node for each path.

/home/pigrenok/.pyenv/versions/3.10.9/envs/pygengraph/lib/python3.10/site-packages/fastcore/docscrape.py:225: UserWarning: Unknown section Notes
  else: warn(msg)

pathNodeDirToCombinedArray

 pathNodeDirToCombinedArray (pathNodeArray, pathDirArray)

Combines path node and direction arrays (provided by graph.initialPathAnalysis function as pathNodeArray and pathDirArray (output indices 1 and 3)).

getNextNodePath

 getNextNodePath (pathNodeArray, pathLengths)

When following path to find next in graph order, it will not necessarily be π‘˜+1 , it can be π‘˜+𝑝 if π‘˜+1,…,π‘˜+π‘βˆ’1 are not passed by the path.

Create a list of all unique node numbers in each path by - Either set(path) - *** or np.unique(path) Preferable option would be selected according to the selection of the options for next step

And then do one of the following - Just leave the list as it is and in a loop check if π‘˜+𝑝 (with 𝑝=1,… ) is in the path (k+p in pathUnique) Probably slow - Sort the list pathUnique and for every node where we need to find next in order just find its index and take the next one pathUnique[pathUnique.index(k)+1] - Create a dict for each node (except the last one) with key as π‘˜ and value as π‘˜+𝑝 . On break in path at position π‘˜ we get the next as dict[k] - *** Do np.diff(np.sort(pathUnique)) and create a dict for each node after which diff!=1 with key as π‘˜ and value as π‘˜+𝑝 . Then when we get to the break in path at position k. We check k in dict and if True, then the next is dict[k], otherwise it is k+1

This function implements generating the dict as described in *** (triple starred) options, which are currently used options.

Main class holding all data and main methods operating the graph

Graph definition, constructor and some utils for constructor

GenomeGraph

 GenomeGraph (gfaPath=None, doOverlapCleaning=True, paths=None,
              nodes=None, nodesData=None, links=None, pathsDict=None,
              sequenceFilesDict=None, annotationFiles=None,
              pangenomeFiles=None, doBack=False, verbose=True, **kwargs)

This is a constructor of a class GenomeGraph. This class allows to hold not only vanila genome graph, but also a lot of extra information and also manipulate graphs in multiple ways, including sorting graphs, adding and deleting nodes and links, and manipulating metadata etc.

At the moment, there are four ways an instance can be created. It depends on what parameters are passed to the constructor. If parameters for more than 1 option is passed, there is a priority queue which constructor follows. Each option and its priorities are provided below.

  • Priority 1: If you pass gfaPath as actual path to gfa file, then it will be loaded as is. In this case, the following options are available:

    accessionsToRemove: list or None (default). If not None, a list of strings, if any of the string contains in pathname, the path will be ignored while loading.

    isGFASeq: boolean (default: True). Whether the graph should be considered as sequence graph (True) or as gene/block graph (False).

  • Priority 2: If nodes, links and paths are provided (not None), they should be as following:

    nodes: list[str]: a list of strings with node IDs (unique)

    links: dict{int:dict{str:list[tuple(int,str)]}}: a dict with keys integers with 1-basednode numbers (in the order as in self.nodes) from which the link starts. Value is a dict with key of directionality of the from node (β€˜+’ for normal direction or β€˜-’ for inverse). Value is a list of tuples with two values: first integer is 1-based node number of node to which the link is going and second string is β€˜+’ or β€˜-’ representing the directionality in which this node is represented in this link.

    β€˜paths’: list[list[str]]: List, which contains a list for each accession/path, which is represented by a list of strings, each of which has a format β€˜{1-based node number}{directionality}’, where {1-based node number} is an integer 1-based number of node using the order as in nodes, {directionality} is either β€˜+’ for normal direction, or β€˜-’ for inverted.

  • Priority 3: if pathsDict is provided then the graph is created from the paths for multiple accessions. pathsDict is a dict{int:list[str]}; keys are names of accessions, and values are lists of strings of the following format β€˜{node name}{directionality}’, where {node name} is identifiable unique name which identifies the node, {directionality} is either β€˜+’ for normal direction, or β€˜-’ for inverted. Note, that this can create no-sequence graph only (e.g. gene graph). Sequences can be added later on through adding sequences to GenomeGraph.nodesData list.

    An extra optional parameter is:

    nodeNameLengths: list[int] or None, a list of alternative node lengths. By default, each node will be represented as a single cell/column, but if provided, variable length can be provided.

  • Priority 4: If annotationFiles is not None, but is a list of paths to annotation (gff3) files, then the following extra options are available:

    β€˜sequenceFilesDict’: a dict{str:str}, where keys are IDs of accessions used in annotation files and value is a path to FASTA file (relative to gff files). Assumption is that FASTA sequence names are the same as GFF3 sequence names.

    pangenomeFiles: a list[str], a list of GFF3 files for the same intervals as in annotationFiles, ID fields in GFF2 metadata should be the same.

    accOrder: list or None (default), Order of accessions in the graph. If None, accessions will be sorted in alphabetical order.

    chromosome: str or None, if None, create one graph for all chromosomes (not fully implemented, see manual), otherwise, create only one graph for given chromosome.

    doUS: boolean (default: False) Add unrelated sequence blocks between annotated genes/blocks.

    refAnnotationFile: str. If given, it has to be a path to gff3 file with reference annotation with reference notation for gene names. In main annotations reference gene names should be identified either in β€œgene” records under β€œAT” key (prioritised), or in β€œmRNA” records under β€œName” key (fallback). If ATMap is provided, then

    refSequenceFile: str or None (default). If provided with path, then it will be used to obtain sequences of each block/gene.

    transmapFile: a tab delimited file with column β€œOrthogroup”, which contains similarity ID for genes and a column with name given by transmapFileRefCol which contains reference annotation gene names.

    transmapFileRefCol: str or None, a name of column for reference gene names in transmapFile

    refAccession: str or None (default). Accession ID for reference annotation (should be provided if refAnnotationFile if provided).

Import graph from file

From GFA

From Paths

From annotations

Import or create annotation

From annotation files

GenomeGraph.loadAnnotations

 GenomeGraph.loadAnnotations (annotationPath, seqSuffix)

This function should allow adding interval metadata (annotations) to sequence (nucleotide) graph. It has never been properly tested.

Create annotation from nodes (artificial annotation)

GenomeGraph.updateAnnotationFromNodes

 GenomeGraph.updateAnnotationFromNodes (isSeq=True)

This function is used only for primitive block graphs (e.g. gene and chain graphs) if there is no proper annotation available (e.g. graph was created from paths and some extra information about nodes is needed).

It takes β€œname” of each node either from graph.nodes (if isSeq is False) or from graph.nodesData (if isSeq is True).

Parameters ##########

isSeq: Whether it contains names as names or as seq.

Graph sorting

GenomeGraph.generateTremauxTree

 GenomeGraph.generateTremauxTree (byPath=True)

This function generates Tremaux tree for our graph. It is not a simple Tremaux tree and requires an adjustment process, which happens inside the TremauxTree class constructor.

GenomeGraph.treeSort

 GenomeGraph.treeSort (byPath=True, bubblePriorityThreshold=0.5)

This is the main function for sorting graph. It requires some further work, but works relatively well as is.

Export/Save (to GFA)

GenomeGraph.toGFA

 GenomeGraph.toGFA (gfaFile, doSeq=True)

Recording existing graph structures to GFA v1 file + some json and joblib files to preserve extra metadata.

GenomeGraph.addAccessionAnnotation

 GenomeGraph.addAccessionAnnotation (annotationFile, sequenceFile=None)

Ideally, a function should be able to add one accesstion to existing graph. When implemented, _graphFromAnnotation should be using this function.

GenomeGraph.addNode

 GenomeGraph.addNode (nodeID, data=None)

Need testing. Again, there is no point of adding a node to a graph if this node will not be present in any of the paths.

Node Inversion

GenomeGraph.invertNodes

 GenomeGraph.invertNodes ()

This function look at inversion/strand of each node in each path. If more than half of paths passing node in β€œinverted” direction, then inverstion should be switched over (currently inverted passes should become normal and normal passes should become inverted.) It is done every time a graph is loaded. Possibly, it should be possible to not doing it as it will take a lot of time for larger graphs.

Node removal methods

GenomeGraph.removeNodes

 GenomeGraph.removeNodes (nodeIDsToRemove)

This (and related to it) function allows removal of a node from a graph. In normal situation, it should not be done to a graph as it will make it invalid in most cases, it is very important functionality for removal of overlaps (see below).

Node Substitution (not implemented)

Node overlap removal

GenomeGraph.removeOverlaps

 GenomeGraph.removeOverlaps ()

When the graph (nucleotide only) is loaded, overlaps are allowed and can be provided using CIGAR strings (it will not be checked). Such overlaps can appear for instance when a compacted de bruijn graph is used (e.g. generated by CUTTLEFISH). They should be removed in order to make the graph not artificially overcomplicated.

Unfortunately, this current implementation is not working properly and needs to be looked at in details. It is most probably overcomplicated and overthought. It should be relatively easy to do.

Utility methods