mgkit.workflow.add_gff_info module¶
Add more information to GFF annotations: gene mappings, coverage, taxonomy, etc..
Uniprot Command¶
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.
Note
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
Coverage Command¶
Adds coverage information from BAM alignment files to a GFF file, using the
function 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
arguments.
Hint
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 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.
Note
From version 0.4.2 the behaviour is different in the following ways:
- the GFF file is loaded in memory to get information that helps reduce
the memory size
- the depth file is scanned and each time a sequence in the GFF file is
found, coverage data is added, this way memory is freed more often
- the GFF file is written at the very end, so chaining multiple commands
(one for each sample) doesn’t give any advantage now
- there’s no need to use the -aa option of samtools depth, since some
assumptions are made when scanning the fils (see next point)
- when the depth file is scanned completely, the sequences in the GFF file
with coverage information are set coverage 0.0, instead of raising an error
average coverage is only calculated if requested, with the -m option
The GFF file is first loaded in memory to get information about the maximum array size for each sequence, then the script proceeds for each depth file, to scan it until a) all sequences in the GFF are found or b) until the end of the file. If b) is the case, all sequences with not found in the depth file are set to a coverage of 0.0. A warning with the number of sequences lacking coverage is printed, for diagnostic purposes. If the the mean coverage is requested with the -m option, coverage is then calculated
To create samtools depth files, you can use this command:
$ samtools depth [bam_file] > [depth_file]
where [bam_file] is the alignment to be used and [depth_file] where the
file will be stored (by redirection). You could create a gzipped file, as
the file can be huge (and in my experience compressing reduces by ~10x the
file size)::
$ samtools depth [bam_file] | gzip > [depth_file.gz]
Warning
if version < 0.4.2:
The -aa options must be used to pass information about all base pairs and sequences coverage in the BAM/SAM file.
$ 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.
Kegg Information¶
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
done using 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
dictionary
HDF5
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.
Note
the dictionary content is loaded after the table files and its keys and corresponding values takes precedence over the text files.
Warning
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
mgkit.io.blast.parse_accession_taxa_table()
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.
Adding information about Bins¶
Warning
The taxonomy used must have been created using the taxon-utils import command on the same DB the results is based upon. At the moment, only PhyloPhlan is supported
The script needs to know to which Bin the annotation belongs to, which can be defined at the moment in 3 different ways (the priority given is in reverse order):
passing the Bin ID defined in the first column of PhyloPhlan results
pass the attribute in the GFF file that contains the Bin ID
the seq_id (FASTA headers) contain the Bin ID
The first way to operate relies on -i to indicate the exact bin_id to use. This can only be passed once, so the whole input GFF will set to that Bin ID.
The secondo option -a relies on an attribute being set in the GFF, this can be used if annotations in the GFF file comes from different Bins
The third option -h relies on the fact that the sequences were renames in a way that can help find the Bin ID. It relies on the sequences being renamed by fasta-utils rename in a form similar to NCBI. By default expects something like this as seq_id:
BinID.fa|contig_1221|xXAS
The script will split the header -s ‘|’, then picking the first element -x 0 and finally replacing -e .fa from it:
BinID.fa|contig_1221|xXAS -> BinID.fa -> BinID
Changes¶
Changed in version 0.5.7: added bins command and –manual option to script
Changed in version 0.4.2: added -m option for cov_samtools command, to calculate the average coverage for an annotation (cov attribute). Fixed loading compressed files
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.
-
mgkit.workflow.add_gff_info.add_bin_information(annotation, bin_info, taxon_id)[source]¶ New in version 0.5.7.
Adds bin information to the GFF annotation
-
mgkit.workflow.add_gff_info.add_uniprot_info(annotations, email, force_taxon_id, taxon_id, lineage, eggnog, enzymes, kegg_orthologs, protein_names, mapping, info_cache)[source]¶