Curation#

MATLAB functions in RAVEN/src/curation of the RAVEN toolbox.

Functions#

Function	Summary
`curateModelFromTables`	Curate or add mets, rxns and genes from tables.
`downloadGenomeData`	Download genome annotation files from NCBI.
`getGeneData`	Build a gene-mapping table from NCBI genome annotation files.
`processProteinFastaFile`	Rename protein FASTA headers using a gene mapping table.
`renameModelGenes`	Replace gene identifiers in a RAVEN model.

Reference#

curateModelFromTables#

function newModel=curateModelFromTables(model,metsInfo,varargin)

Curate or add mets, rxns and genes from tables.

Curates existing metabolites, reactions and genes and adds new ones from tabular data files.

If the *.tsv files contain metabolites, reactions or genes that are already present in the model, the information in the model is overwritten. This includes empty annotations in the *.tsv files, unless 'overwriteWithEmpty' is false. Metabolites are matched by metaboliteName[comp]; reactions by the stoichiometry of their reactants and products; genes by their gene name. Entities that are not present in the model are added.

Columns are identified by their header (case-insensitive, and spaces around a header are ignored, also inside its quotes), so their order does not matter and optional columns can be left out. Every field that a table supplies is written to new and existing entities alike; a field the model does not have yet is created, with its default value for the other entities, unless all supplied values for it are empty. An empty grRules cell leaves the rule of an existing reaction unchanged, and so does an empty lb, ub or rev cell: an empty bound means that no bound is given. A new reaction without a lb gets a lower bound of -Inf, and one without an ub gets model.annotation.defaultUB (Inf if the model has none).

Parameters

model : struct: RAVEN model structure to be curated.
metsInfo : char: Path to a *.tsv file with metabolite information, or 'none' to skip metabolite curation. Columns: metNames, comps (both required), formula, charge, inchi, metNotes, and any number of MIRIAM-namespace columns. Singular forms and the RAVEN field names (metFormulas, metCharges, inchis) are accepted as headers too.

Name-Value Arguments

genesInfo : char: Path to a *.tsv file with gene information, or 'none'. Columns: genes (required), geneShortNames (or geneShortName), and MIRIAM.
rxnsCoeffs : char: Path to a *.tsv file with reaction stoichiometric coefficients, or 'none'. Columns (all required): rxnIdx (or index), rxnNames, metNames, comps, coefficient. One row per (reaction, metabolite) pair.
rxnsInfo : char: Path to a *.tsv file with reaction information, or 'none'. Columns: rxnIdx (or index), rxnNames (both required), grRules, lb, ub, rev, subSystems, eccodes, rxnNotes, rxnReferences, rxnConfidenceScores, and MIRIAM. An index appears in at most one row of this table, otherwise RAVEN:badInput is raised (the rxnsCoeffs file has one row per metabolite of the reaction and repeats it). The two reaction tables are matched on index and name, and reactions or indices that appear in only one of them raise RAVEN:badInput.
metPrefix : char: Prefix used to create new metabolite ids (for example 's_' for yeast-GEM, 'M_' for the BiGG default) (default 'M_').
rxnPrefix : char: Prefix used to create new reaction ids (default 'R_').
overwriteWithEmpty : logical: if true, an empty cell in a row of an existing metabolite, gene or reaction empties the corresponding field of the model. If false, an empty cell leaves the existing value unchanged and only the cells with a value are written (default true; [] also means the default). Must be a scalar logical, 0 or 1. A cell with only blanks counts as empty in both modes. Empty grRules cells never overwrite. An empty or NaN lb, ub or rev cell always keeps the existing value. For the charge and the confidence score, an empty or NaN cell is treated alike: when overwriteWithEmpty is true and at least one other row of the column has a value, it becomes 0; a column in which no row has a value is not written, and the existing values stay. The same holds for a MIRIAM column: a column in which no row of an existing entity has a value does not clear the annotations.

Returns

newModel : struct: Curated RAVEN model structure.

Examples

newModel = curateModelFromTables(model, metsInfo, genesInfo, ...
                rxnsCoeffs, rxnsInfo, metPrefix, rxnPrefix);

Notes

In all three info tables, any column whose header is not one of the listed core fields is treated as a MIRIAM annotation namespace and stored on the matching entity.

Cell values: the charge, confidence score, lb, ub, rev, rxnIdx and coefficient cells are read as plain decimal numbers (digits with an optional sign, decimal point and exponent), Inf or NaN in any case. Only the lb, ub and rev cells also accept true and false in any case (read as 1 and 0). Any other text, including text with a comma, hexadecimal numbers and expressions, raises RAVEN:badInput; no expression is evaluated. A charge must be a finite whole number (or empty or NaN). A rev must be 0, 1, true or false (or empty or NaN). The rev column sets the reversibility of existing reactions only, independently of lb and ub, which are written as given; for new reactions the reversibility is derived from lb and ub.

Gene rules: a grRules cell must be a well-formed gene rule (genes joined by the words and or or in any case, or by &, &&, | or ||, with balanced parentheses); a rule with an operator that lacks an operand, empty parentheses or two genes without an operator between them raises RAVEN:badInput before the model is changed. A valid rule is stored in RAVEN's standard format. Whether the genes exist is checked when the rule is applied.

Genes: a gene id that the genesInfo table lists more than once raises RAVEN:badInput, whether the model has the gene or not.

Bounds: a lower bound above the upper bound is written as given. solveLP rejects such a model when it is run ("Invalid set of bounds for reaction(s): ...").

Reaction tables: a row of the rxnsCoeffs table without a rxnIdx, rxnNames, metNames and comps value is ignored; a row with any of them needs a rxnIdx and a coefficient that is a finite number other than 0. A row of the rxnsInfo table without any value is ignored; a row with a value needs a rxnIdx. A rxnIdx is a whole number of at least 1; the indices (and names) of the two reaction tables must match, with indices compared by their numeric value (1, 1.0 and 1e0 are the same index). Otherwise RAVEN:badInput is raised.

downloadGenomeData#

function [gffFile, faaFile] = downloadGenomeData(accession, varargin)

Download genome annotation files from NCBI.

Downloads the GFF3 annotation and protein FASTA (.faa) files for an NCBI genome assembly accession using the NCBI Datasets v2 API, which returns both files in a single archive. The files are saved as _genomic.gff and _protein.faa in outputDir. They are input files for getGeneData.

Parameters

accession : char: NCBI genome assembly accession, for example 'GCF_000002595.2'. The prefix is 'GCF_' (RefSeq) or 'GCA_' (GenBank).

Name-Value Arguments

outputDir : char: Directory where downloaded files are saved (default: current working directory).
verbose : logical: Print download progress to the command window (default: true).

Returns

gffFile : char: Full path to the downloaded GFF3 annotation file.
faaFile : char: Full path to the downloaded protein FASTA file.

Examples

% Download to current directory
[gff, faa] = downloadGenomeData('GCF_000002595.2');

% Save to a specific directory
[gff, faa] = downloadGenomeData('GCF_000002595.2', 'data/');

% Suppress progress messages
[gff, faa] = downloadGenomeData('GCF_000002595.2', 'data/', false);

Notes

An internet connection is required. When both files are already present in outputDir, no download takes place; deleting them forces a new download. An error is raised when NCBI Datasets does not return both a GFF3 file and a protein FASTA file for the accession.

getGeneData#

function geneTable = getGeneData(accession, varargin)

Build a gene-mapping table from NCBI genome annotation files.

Parses the GFF3 annotation file of an NCBI genome assembly and produces a table mapping locus tags to gene symbols, protein IDs and other identifiers. The table can be passed to renameModelGenes to change gene identifiers in a RAVEN model. If the input is an assembly accession, downloadGenomeData downloads the GFF3 and protein FASTA (.faa) files to downloadDir, unless both are already present there.

Parameters

accession : char: NCBI genome assembly accession, for example 'GCF_000002595.2', or the path to a local GFF3 annotation file. An error is raised if the input is neither an existing file nor an accession of the form 'GCF_' or 'GCA_' followed by digits and an optional version.

Name-Value Arguments

outputFile : char: Path of the tab-delimited .tsv file to which the table is written. No file is written if omitted or empty (default '').
downloadDir : char: Directory where genome files are downloaded if not already present. Only used when the input is an accession (default: current working directory).

Returns

geneTable : table

MATLAB table with one row per distinct combination of gene and protein accession (a CDS feature with a protein accession, repeated CDS rows are merged), with the columns below. A column entry is empty when the annotation lacks the value.

locus_tag : locus identifier (for example 'Cre01.g000001')
old_locus_tag : earlier locus identifier
GeneID : NCBI Gene ID (for example '5723799')
gene_name : gene symbol (for example 'rbcL')
GenBank_protein : protein accession (for example 'XP_001698190.2'), as in the protein FASTA headers
UniProt : UniProt accession from the Dbxref attribute of the gene

Examples

% Fetch and parse gene data for Chlamydomonas reinhardtii
geneTable = getGeneData('GCF_000002595.2');

% Save the mapping table to a specific file
geneTable = getGeneData('GCF_000002595.2', 'chlamy.tsv');

% Save to a file and download to a specific directory
geneTable = getGeneData('GCF_000002595.2', 'chlamy.tsv', 'data/');

Notes

When outputFile is provided the table is written as UTF-8 encoded tab-delimited text with a header row, which readtable and renameModelGenes can read. A CDS whose gene has no locus_tag, or whose gene feature is missing, gives a row with an empty locus_tag.

processProteinFastaFile#

function processedFaaFile = processProteinFastaFile(faaFile, geneTable, headerCol, varargin)

Rename protein FASTA headers using a gene mapping table.

Reads a protein FASTA file (such as one downloaded by downloadGenomeData) and replaces each sequence header with the value from the specified geneTable column. A sequence is matched to a row through the GenBank_protein accession, which is the first token of the original FASTA header. Sequences whose accession is not found in geneTable keep their original header. The numbers of renamed and kept sequences are printed.

Parameters

faaFile : char: Path to the protein FASTA file (.faa) to process.
geneTable : char | table: Either a path to a tab-delimited .tsv file produced by getGeneData, or a MATLAB table variable. It must contain the columns 'GenBank_protein' and the column named by headerCol.
headerCol : char: Name of the geneTable column whose values will replace each FASTA header (for example 'locus_tag', 'gene_name', 'GenBank_protein').

Name-Value Arguments

outputDir : char: Directory where the processed FASTA file is saved. It is created if missing. The output file name is the original base name with '_processed' appended before the extension (default: current working directory).

Returns

processedFaaFile : char: Full path to the written processed FASTA file.

Examples

% Rename headers using locus_tag
[~, faa]  = downloadGenomeData('GCF_000002595.2');
geneTable = getGeneData('GCF_000002595.2', 'chlamy.tsv');
outFile   = processProteinFastaFile(faa, geneTable, 'locus_tag');

% Load both inputs from saved files, use gene_name as header
outFile = processProteinFastaFile('GCF_000002595.2_protein.faa', 'chlamy.tsv', 'gene_name');

% Save the processed file to a specific directory
outFile = processProteinFastaFile(faa, 'chlamy.tsv', 'locus_tag', 'results/');

Notes

GenBank_protein must contain accessions equal to the first whitespace-delimited token of each FASTA header (after '>'), for example WP_012345678.1. getGeneData takes the accession from the protein_id attribute of the CDS features of the GFF3 annotation, or from a GenBank entry of the Dbxref attribute when protein_id is empty.

renameModelGenes#

function model = renameModelGenes(model, geneTable, fromCol, toCol)

Replace gene identifiers in a RAVEN model.

Replaces the identifiers in model.genes and model.grRules that are listed in column fromCol of a mapping table with the identifiers in column toCol. The table is produced by getGeneData or supplied manually. After renaming, model.rxnGeneMat is rebuilt.

Parameters

model : struct: RAVEN model struct. It must contain model.genes and model.grRules.
geneTable : table | char: Gene-mapping table, either (a) a MATLAB table variable, or (b) the path to a tab-delimited .tsv file, which is read with readtable.
fromCol : char: Name of the column in geneTable whose values match the identifiers in model.genes (for example 'locus_tag').
toCol : char: Name of the column in geneTable with the new identifiers (for example 'gene_name').

Returns

model : struct

Model struct with these fields changed:

model.genes : gene list with the new identifiers.
model.grRules : gene-reaction rules with the new identifiers.
model.rxnGeneMat : sparse matrix rebuilt for the new gene list.

Examples

model     = readYAMLmodel('iCre1355.yml');
geneTable = getGeneData('GCF_000002595.2', 'chlamy.tsv');
model     = renameModelGenes(model, geneTable, 'locus_tag', 'gene_name');

% Supply a TSV file path directly
model = renameModelGenes(model, 'chlamy.tsv', 'locus_tag', 'gene_name');

Notes

Genes with no entry in fromCol, or whose toCol value is empty, are left unchanged, and a warning lists them.
If a fromCol value occurs in several rows, the first row with a non-empty toCol value is used. Values are trimmed of surrounding whitespace.
An identifier is replaced only where it is delimited by whitespace, parentheses or the start or end of the rule, so 'gene1' is not replaced inside 'gene10'.
Unnecessary outer parentheses are removed from the grRules and the grRules are standardised with standardizeGrRules before and after renaming. A warning is issued when standardizeGrRules reports rules with potential issues.