Missing gene - blastn evalue rounded to 0.0

[crcardenas]

I am running mitofinder on assembled contigs to search for mtDNA/Genomes. However, I am missing CO1 from the output. I thought maybe I just didn't capture that mtGene. However, to confirm my suspicions, I ran blastn using the same 5 reference genes and found a contig or two that should contain that gene (three hits actually, but the longest blast alignment fits the expected length of 1.5k bp). For some reason the tmp files directory generated by mitofinder has the CO1 blast output empty (actually all of them are) so I cannot check the blast scores used by mitofinder.

Here is an example of my sanity check where all CO1 evalues are 0.0

contig_2018	NC_036507.1@COX1	91.063	884	77	2	1	883	364	1246	0.0	1194
contig_2018	NC_046469.1@COX1	90.611	884	81	2	1	883	364	1246	0.0	1171
contig_2018	NC_016469.1@COX1	90.498	884	82	2	1	883	364	1246	0.0	1166
contig_2018	NC_018339.1@COX1	89.581	883	92	0	1	883	355	1237	0.0	1122
contig_2018	NC_044759.1@COX1	89.052	886	91	3	1	883	364	1246	0.0	1094
...
contig_1535	NC_036507.1@COX1	90.537	1173	102	4	15804	16968	1	1172	0.0	1543
contig_1535	NC_044759.1@COX1	89.608	1174	115	6	15804	16970	1	1174	0.0	1485
contig_1535	NC_016469.1@COX1	89.632	1167	112	4	15804	16962	1	1166	0.0	1476
contig_1535	NC_046469.1@COX1	89.632	1167	112	4	15804	16962	1	1166	0.0	1476
contig_1535	NC_018339.1@COX1	89.261	1164	116	5	15813	16968	1	1163	0.0	1448
...
contig_2876	NC_016469.1@COX1	90.944	1546	138	2	3549	5093	1	1545	0.0	2078
contig_2876	NC_036507.1@COX1	90.980	1541	137	2	3549	5088	1	1540	0.0	2074
contig_2876	NC_046469.1@COX1	90.962	1538	137	2	3549	5085	1	1537	0.0	2069
contig_2876	NC_044759.1@COX1	90.201	1541	145	3	3549	5086	1	1538	0.0	2004
contig_2876	NC_018339.1@COX1	90.137	1531	151	0	3558	5088	1	1531	0.0	1991

My understanding is that blast will round a very small e-value to 0 if its a good match.

However, I noticed that this could be an issue with the way the max value is added in the python script:
line 84 mitofinder

parser.add_argument('-e', '--blast-eval', help='e-value of blast program used for contig identification and annotation. Default = 0.00001', type=float,
						default=0.00001, dest='blasteVal')

I was also missing some other genes before adding mtDNA reference from a sister genus, which I think might explain why some of the genes were missing. Here is a good example, ND4: ~ 1.3k bp reference genes. Which now appear with near zero evalues for the longest alignment length contig.

contig_2359	NC_016469.1@ND4	90.280	751	73	0	259262	260012	1290	540	0.0	983
contig_2359	NC_036507.1@ND4	90.280	751	73	0	259262	260012	1290	540	0.0	983
contig_2359	NC_046469.1@ND4	89.614	751	78	0	259262	260012	1290	540	0.0	955
contig_2359	NC_044759.1@ND4	89.418	756	76	4	259262	260017	1290	539	0.0	950
contig_2359	NC_018339.1@ND4	89.390	754	78	2	259262	260015	1290	539	0.0	948
...
contig_3886	NC_016469.1@ND4	89.389	311	32	1	837341	837651	839	1148	1.64e-107	390
contig_3886	NC_018339.1@ND4	89.389	311	32	1	837341	837651	839	1148	1.64e-107	390
contig_3886	NC_036507.1@ND4	89.644	309	29	3	837341	837648	839	1145	1.64e-107	390
contig_3886	NC_046469.1@ND4	88.746	311	34	1	837341	837651	839	1148	3.56e-104	379
contig_3886	NC_044759.1@ND4	87.539	321	36	4	837335	837653	832	1150	7.70e-101	368

I would normally try tweaking the code myself, but I am using the singularity container and don't know how to edit anything there.

Any advice would be greatly appreciated!

---

To clarify, I am using oxford nanopore output, so this might be part of the issue. However, I have kept the max length for contigs as default (25k bp); so I do not think it would be in conflict with the framework of using mitofinder for UCE assemblies. It may still not be appropriate, but I thought you would want to be made aware of this issue.

system information:
ubuntu server, mitofinder version: mitofinder_v1.4.2.sif, singularity version: 3.7.1
for my own blastn: BLAST 2.12.0+

[RemiAllio]

Hi!

Thank you for your feedback.
I am a little surprised by what you found.

Is it possible for you to share the listed contigs, your command line, and your reference so I can check all of this?

Best,
Rémi

[crcardenas]

Hello Rémi,

Interestingly, when I subset it like this, I end up with COX1. However, I still see disagreement between my blastn output and mitofinder:
I see the following genes in my blast: COX1, COX2, COX3, ND2, ND3, ND4, ND5, ND6, and rrnS. However, mitofinder recovers these genes in addition to ND1,ND4L, and CYTB; but not rrnS.

It looks like blastx is used at some point, so that might explain part of this but; so I think my initial prediction about the 0.0 value was inaccurate.

However, I cannot be certain given that the mitofinder blast files are inconsistently empty ( those found in CBX1139_subset/CBX1139_subset_tmp/*_blast_out.txt) and there is inconsistencies between some of those missing genes.

Please see the attached files: CBX1139_subset_files.zip

reference.fasta
-this is the file found in the CBX1139_subset/CBX1139_subset_tmp/contig_id_database.fasta its what I used to create the blast database.
CBX1139_subset.fasta contans any contigs below the default 25k threshold.

Hopefully everything I have provided is clear and I didn't miss anything.

Best,
Cody

Example mitofinder script

#!/bin/bash
source /local/anaconda3/bin/activate
conda activate singularity
cd /home/cody/CALOSOMA_Genomes/NANOPORE/blobtools/CBX1139/mitofinder
for i in 1; do
# can be setup to iterate over a directory of files; cp assemblies over (singularity does not like symlinks) and adjust -B,-j,-a, etc.
/Calosoma_WORKING/MitoFinder/ mitofinder_v1.4.1.sif \
singularity run  \
-B /home/cody/mito/:/home/cody/mito/,/data/raw/Calosoma/NANOPORE/blobtools/CBX1139/mitofinder/:/data/raw/Calosoma/NANOPORE/blobtools/CBX1139/mitofinder/ mitofinder_v1.4.2.sif \
-j CBX1139_subset \
-a /data/raw/Calosoma/NANOPORE/blobtools/CBX1139/mitofinder/CBX1139_subset.fasta \
-r /data/raw/Calosoma/NANOPORE/blobtools/CBX1139/mitofinder/reference.gb \
-m 64 -p 8 -o 5 --rename-contig no \
--min-contig-size 400 --max-contig-size 30000  
done |& tee -a ./mitofinder_`date +%Y%m%d-%H%M%S`.log

example blastn script

for i in 1; do 
#see previous comment about iteration, less necessary here.
DB="../CBX1139/mitofinder/find_COI"; 
blastn -db ${DB}/test.fasta -query CBX1139_subset.fasta \
-outfmt "6" -max_target_seqs 10  -evalue 1e-45 -num_threads 3 \
-out blastn_mito_CBX1139_blast.out; 
done