Add more information to GFF annotations: gene mappings, coverage, taxonomy, etc..
If the gene_id of an annotation is a Uniprot ID, the script queries Uniprot for the requested information. At the moment the information that can be added is the taxon_id, taxon_name, lineage and mapping to EC, KO, eggNOG IDs.
It’s also possible to add mappings to other databases using the -m option with the correct identifier for the mapping, which can be found at this page; for example if it’s we want to add the mappings of uniprot IDs to BioCyc, in the abbreviation column of the mappings we find that it’s identifier is REACTOME_ID, so we pass -m REACTOME to the script (leaving _ID out). Mapped IDs are separated by commas.
The taxonomy IDs are not overwritten if they are found in the annotations, the -f is provided to force the overwriting of those values.
See also MGKit GFF Specifications for more informations about the GFF specifications used.
As the script needs to query Uniprot a lot, it is recommended to split the GFF in several files, so an error in the connection doesn’t waste time.
However, a cache is kept to reduce the number of connections
Adds coverage information from BAM alignment files to a GFF file, using the
mgkit.align.add_coverage_info(), the user needs to supply for
each sample a BAM file, using the -a option, whose parameter is in the form
sample,samplealg.bam. More samples can be supplied adding more -a
As an example, to add coverage for sample1, sample2 the command line is:
add-gff-info coverage -a sample1,sample1.bam -a sample2,sample2.bam \ inputgff outputgff
A total coverage for the annotation is also calculated and stored in the cov attribute, while each sample coverage is stored into sample_cov as per MGKit GFF Specifications.
Adding Coverage from samtools depth¶
The cov_samtools allows the use of the output of samtools depth command. The -aa options must be used to pass information about all base pairs and sequences coverage in the BAM/SAM file. The command work similarly to coverage, accepting compressed depth files as well. If only one depth file is passed and no sample is passed, the attribute in the GFF will be cov, otherwise the attribute will be sample1_cov, sample2_cov, etc.
To create samtools depth files, this command must be used:
$ samtools depth -aa bam_file
Uniprot Offline Mappings¶
Similar to the uniprot command, it uses the idmapping file provided by Uniprot, which speeds up the process of adding mappings and taxonomy IDs from Uniprot gene IDs. It’s not possible tough to add EC mappings with this command, as those are not included in the file.
The kegg command allows to add information to each annotation. Right now the information that can be added is restricted to the pathway(s) (reference KO) a KO is part of and both the KO and pathway(s) descriptions. This information is stored in keys starting with ko_.
Expected Aminoacidic Changes¶
Some scripts, like snp_parser - SNPs analysis, require information about the expected
number of synonymous and non-synonymous changes of an annotation. This can be
mgkit.io.gff.Annotation.add_exp_syn_count() by the user of the
command exp_syn of this script. The attributes added to each annotation are
explained in the MGKit GFF Specifications
Adding Count Data¶
Count data on a per-sample basis can be added with the counts command. The accepted inputs are from HTSeq-count and featureCounts. The ouput produced by featureCounts, is the one from using its -f option must be used.
This script accept by default a tab separated file, with a uid in the first column and the other columns are the counts for each sample, in the same order as they are passed to the -s option. To use the featureCounts file format, this script -e option must be used.
The sample names must be provided in the same order as the columns in the input files. If the counts are FPKMS the -f option can be used.
Adding Taxonomy from a Table¶
There are cases where it may needed or preferred to add the taxonomy from a gene_id already provided in the GFF file. For such cases the addtaxa command can be used. It works in a similar way to the taxonomy command, only it expect three different type of inputs:
- GI-Taxa table from NCBI (e.g. gi_taxid_nucl.dmp, )
- tab separated table
The first two are tab separated files, where on each line, the first column is the gene_id that is found in the first column, while the second if the taxon_id.
The third option is a serialised Python dict/hash table, whose keys are the gene_id and the value is that gene corresponding taxon_id. The serialised formats accepted are msgpack, json and pickle. The msgpack module must be importable. The option to use json and msgpack allow to integrate this script with other languages without resorting to a text file.
The last option is a HDF5 created using the to_hdf command in taxon-utils - Taxonomy Utilities. This requires pandas installed and pytables and it provides faster lookup of IDs in the table.
While the default is to look for the gene_id attribute in the GFF annotation, another attribute can be specified, using the -gene-attr option.
the dictionary content is loaded after the table files and its keys and corresponding values takes precedence over the text files.
from September 2016 NCBI will retire the GI. In that case the same
kind of table can be built from the nucl_gb.accession2taxid.gz file
The format is different, but some information can be found in
Adding information from Pfam¶
Adds the Pfam description for the annotation, by downloading the list from Pfam.
The options allow to specify in which attribute the ID/ACCESSION is stored (defaults to gene_id) and which one between ID/ACCESSION is the value of that attribute (defaults to ID). if no description is found for the family, a warning message is logged.
Changed in version 0.3.4: removed the taxonomy command, since a similar result can be obtained with taxon-utils lca and add-gff-info addtaxa. Removed eggnog command and added option to verbose the logging in cov_samtools (now is quiet), also changed the interface
Changed in version 0.3.3: changed how addtaxa -a works, to allow the use of seq_id as key to add the taxon_id
Changed in version 0.3.0: added cov_samtools command, –split option to exp_syn, -c option to addtaxa. kegg now does not skip annotations when the attribute is not found.
Changed in version 0.2.6: added skip-no-taxon option to addtaxa
Changed in version 0.2.5: if a dictionary is supplied to addtaxa, the GFF is not preloaded
Changed in version 0.2.3: added pfam command, renamed gitaxa to addtaxa and made it general
Changed in version 0.2.2: added eggnog, gitaxa and counts command
Changed in version 0.2.1.
- added -d to uniprot command
- added cache to uniprot command
- added kegg command (cached)
Changed in version 0.1.16: added exp_syn command
Changed in version 0.1.15: taxonomy command -b option changed
Changed in version 0.1.13.
- added –force-taxon-id option to the uniprot command
- added coverage command
- added taxonomy command
- added unipfile command
New in version 0.1.12.
add_uniprot_info(annotations, email, force_taxon_id, taxon_id, lineage, eggnog, enzymes, kegg_orthologs, protein_names, mapping, info_cache)¶
parse_hdf5_arg(ctx, param, values)¶
split_sample_alg(ctx, param, values)¶
Split sample/alignment option