merenlab
`anvi-export-locus` cannot handle non-sequential gene caller ids
anvi-export-locus in its current form cannot handle gene calls where non-prodigal gene call software does not order numerically. Tests need to be run to see how different external gene call imports can run with the current implementation of anvi-export-locus.
$ anvi-self-test --version
Anvi'o version ...............................: esther (v6.1-master)
Profile DB version ...........................: 31
Contigs DB version ...........................: 14
Pan DB version ...............................: 13
Genome data storage version ..................: 6
Auxiliary data storage version ...............: 2
Structure DB version .........................: 1
$ echo $OSTYPE
darwin18
$ anvi-export-locus -c Staphylococcus_aureus_subsp__aureus_NCTC_8325_GCF_000013425_1-contigs.db \
$ > --num-genes 5,5 \
$ > --search-term Archaeal_23S_rRNA \
$ > --use-hmm \
$ > --hmm-sources Ribosomal_RNAs \
$ > -O s_epidermidis_rrn_locus
Genes to report ..............................: 5 genes before the matching gene, and 5 that follow
Input / Output
===============================================
Contigs DB ...................................: /Users/mschechter/scratch/Staphylococcus_aureus_subsp__aureus_NCTC_8325_GCF_000013425_1-contigs.db
Output directory .............................: /Users/mschechter/scratch
Rev-comp the locus sequence if necessary .....: True
Initialization bleep bloops
===============================================
Mode .........................................: HMM search
Search term ..................................: ['Archaeal_23S_rRNA']
HMM sources being used .......................: Ribosomal_RNAs
Matching genes ...............................: 5 genes matched your search
Exporting locus 1 of 5
===============================================
Contig name ..................................: NC_007795
Contig length ................................: 2,821,361
Num genes in contig ..........................: 2,777
Target gene call .............................: 2,767
Target gene direction ........................: Reverse
First and last gene of the locus (raw) .......: 2762 and 2772
First and last gene of the locus (final) .....: 2762 and 2772
Locus gene call start/stops excess (nts) .....: 2,816,354
Reverse complementing everything .............: True
Output FASTA file ............................: /Users/mschechter/scratch/s_epidermidis_rrn_locus_0001.fa
HERE WE ARE
/Users/mschechter/scratch/s_epidermidis_rrn_locus_0001_sequence.fa
Config Error: You are in big trouble :( The contig name 'NC_007795' in your external gene
callers file does not appear to be in the contigs FASTA file. How did this
happen?
WD
/Users/mschechter/scratch
Original command
anvi-export-locus -c Staphylococcus_aureus_subsp__aureus_NCTC_8325_GCF_000013425_1-contigs.db
--num-genes 5,5
--search-term Bacterial_16S_rRNA
--use-hmm
--hmm-sources Ribosomal_RNAs
-O s_epidermidis_rrn_locus
Genes to report ..............................: 5 genes before the matching gene, and 5 that follow
Input / Output
===============================================
Contigs DB ...................................: /Users/mschechter/scratch/Staphylococcus_aureus_subsp__aureus_NCTC_8325_GCF_000013425_1-contigs.db
Output directory .............................: /Users/mschechter/scratch
Rev-comp the locus sequence if necessary .....: True
Initialization bleep bloops
===============================================
Mode .........................................: HMM search
Search term ..................................: ['Bacterial_16S_rRNA']
HMM sources being used .......................: Ribosomal_RNAs
WARNING
===============================================
There aren't any gene calls that match to the criteria you provided to anvi'o
export locus magic. Is this yet another case of you did everything right yet
anvi'o failed you? If that's the case, let us know :( This class will quietly
kill this process without reporting any error since a lack of hit may be the
expected outcome of some weird processes somewhere.
Anvio isnt seeing an HMM annotation with the search-term: "Bacterial_16S_rRNA"
Which HMM sources does this contigsDB have?
anvi-export-locus -c Staphylococcus_aureus_subsp__aureus_NCTC_8325_GCF_000013425_1-contigs.db --list-hmm-sources -O test
* Archaea_76 [type: singlecopy] [num genes: 76]
* Bacteria_71 [type: singlecopy] [num genes: 71]
* Protista_83 [type: singlecopy] [num genes: 83]
* Ribosomal_RNAs [type: Ribosomal_RNAs] [num genes: 12]
Are there any HMM hits?
anvi-get-sequences-for-hmm-hits -c Staphylococcus_aureus_subsp__aureus_NCTC_8325_GCF_000013425_1-contigs.db -o hmms
Contigs DB ...................................: Initialized: Staphylococcus_aureus_subsp__aureus_NCTC_8325_GCF_000013425_1-contigs.db (v. 14)
Hits .........................................: 119 hits for 4 source(s)
Mode .........................................: DNA sequences
Genes are concatenated .......................: False
Output .......................................: hmms
Are there any HMM hits to Ribosomal RNAs?
$ grep "Ribosomal_RNAs" hmms
>Archaeal_23S_rRNA___Ribosomal_RNAs___9adc507abe9435b965429fd2e696a16f77403b7fa65ef4d31da005b9 bin_id:Staphylococcus_aureus_subsp__aureus_NCTC_8325_GCF_000013425_1-contigs|source:Ribosomal_RNAs|e_value:0|contig:NC_007795|gene_callers_id:2770|start:450754|stop:453650|length:2896
>Archaeal_16S_rRNA___Ribosomal_RNAs___468c00bbcc8c6555b1b15d7939cfc170c0ce171422305e3fcb5f9906 bin_id:Staphylococcus_aureus_subsp__aureus_NCTC_8325_GCF_000013425_1-contigs|source:Ribosomal_RNAs|e_value:1.8e-228|contig:NC_007795|gene_callers_id:2774|start:448831|stop:450370|length:1539
...
Looks like there are only Archaeal_23S_rRNA and Archaeal_16S_rRNA
Try extracting loci around Archaeal_16S_rRNA
$ anvi-export-locus -c Staphylococcus_aureus_subsp__aureus_NCTC_8325_GCF_000013425_1-contigs.db \
--num-genes 5,5 \
--search-term Archaeal_23S_rRNA \
--use-hmm \
--hmm-sources Ribosomal_RNAs \
-O s_epidermidis_rrn_locus \
--overwrite-output-destinations
Genes to report ..............................: 5 genes before the matching gene, and 5 that follow
Input / Output
===============================================
Contigs DB ...................................: /Users/mschechter/scratch/Staphylococcus_aureus_subsp__aureus_NCTC_8325_GCF_000013425_1-contigs.db
Output directory .............................: /Users/mschechter/scratch
Rev-comp the locus sequence if necessary .....: True
Initialization bleep bloops
===============================================
Mode .........................................: HMM search
Search term ..................................: ['Archaeal_23S_rRNA']
HMM sources being used .......................: Ribosomal_RNAs
Matching genes ...............................: 5 genes matched your search
Exporting locus 1 of 5
===============================================
Contig name ..................................: NC_007795
Contig length ................................: 2,821,361
Num genes in contig ..........................: 2,777
Target gene call .............................: 2,767
Target gene direction ........................: Reverse
First and last gene of the locus (raw) .......: 2762 and 2772
First and last gene of the locus (final) .....: 2762 and 2772
Locus gene call start/stops excess (nts) .....: 2,816,354
Reverse complementing everything .............: True
Output FASTA file ............................: /Users/mschechter/scratch/s_epidermidis_rrn_locus_0001.fa
HERE WE ARE
/Users/mschechter/scratch/s_epidermidis_rrn_locus_0001_sequence.fa
Config Error: You are in big trouble :( The contig name 'NC_007795' in your external gene
callers file does not appear to be in the contigs FASTA file. How did this
happen?
Above is the error where I think anvi-export-locus is not able to work properly with gene-calls that are not numbered numerically.
Above is the error where I think anvi-export-locus is not able to work properly with gene-calls that are not numbered numerically.
By this I assume you mean that the gene caller ids are not ordered sequentially in ascending order (with interval of 1).
I can create another instance in which things fail when this is the case. Here is a minimum reproducible example, that uses the infant gut data as a starting point. I swap gene id 5 with gene id 50. Then I try and export a locus that includes genes 4-10. Since gene 5 is now way further up in the contig, things fail due to assumed logic in anvi-export-locus:
#! /usr/bin/env bash
anvi-export-gene-calls --gene-caller prodigal -o external.txt -c CONTIGS.db
anvi-export-contigs -c CONTIGS.db -o fasta.fa
python -c "import pandas as pd;df=pd.read_csv('external.txt',sep='\t');df.iloc[5,0]=50;df.iloc[50,0]=5;df.to_csv('external.txt',sep='\t',index=False)"
anvi-gen-contigs-database -o TEST.db --external-gene-calls external.txt -f fasta.fa
anvi-export-locus -c TEST.db -O LOCUS --gene-caller-ids 4 --num-genes 6 --overwrite-output-destinations
This yields the error:
Config Error: utils.get_most_likely_translation_frame :: sequence has a length less than 3 so
there is nothing to translate.
Should we prevent this from happening at the source, @meren, during contigs database creation? IIRC, you wrote a "external gene calls" sanity check.
Should we prevent this from happening at the source, @meren, during contigs database creation?
I think for each gene caller source, gene calls must be sequential increasing by 1. It is likely that in the future the contigs database will support multiple gene callers to add their stuff. Then we will have a set of gene calls that go from 0 to X, and perhaps another set of gene calls that go from X to Y. So we shouldn't expect them to start from 0, but we should expect them to increment by one without any blank.
If that is the expected behavior, then we have a bug, because the above code successfully makes a contigs db with external genes that do not have sequentially increasing numbers. It goes 0,1,2,3,4,50,6,7,...,49,5,51,52...
A proper solution would also need to check start / stop positions (or start + ((stop - start) / 2) positions?) to make sure it is not just shuffling of the external gene calls file with properly ordered gene caller ids and start/stop positions, but it is a major snafu :(