incorrect hits with bad alignment still gives "1.0" read assignment probability

[Robvh-git]

Dear,

We are analyzing nematode full 18S data (Nanopore reads) and are stumbling into some peculiar results. We use a curated database which only contains 18S nematode sequences with a few 'outgroups' (i.e. non-nematode).

In our soil nematode data, we find a lot of reads of a certain marine worm, which is not possible. The reference of this marine worm is correct in our database (checked).

When I run emu with outputting the read assignment distributions, I find the reads that have a score of "1.0" for this marine worm.

When I filter this read from the emu input data, and use BLAST to identify it, it shows it is actually some insect which is not in our curated emu reference database. We found out later that there indeed could be some insect DNA in this set of samples.

Then I aligned (in Geneious Prime) this read (which had a 1.0 prob. of the marine worm) with the reference read of that marine worm, the alignment is actually quite poor with a similarity of 77%.

I'm a bit worried by this result, as I don't understand why Emu (minimap2) would still not put this read in the "unclassified" as it is clearly not a taxon from our database. I think this happens because we have a limited database and Emu is trying to force it into a taxon, eventhough it is not correct. Is there a way to stop this forcing?

Do you know why this read is not classified as "unclassified/unassigned" and if there is a way to add some threshold parameter to filter poor aligned reads?

[kdc10]

Can you look at the minimap alignment and see type of score minimap gives these reads when aligned to your database?

Emu will classify reads as the database sequences that minimap generates an alignment. Unassigned in emu = unmapped in minimap2. So, if you don't want the classification to show up via emu, then you will need to the minimap alignment restriction such that it does not return an alignment. I suggest you look at the minimap alignment and the minimap2 parameters to try to find which parameter adjustment would meet your needs.

[Robvh-git]

@kdc10 thank you for the quick reply. I dived a bit deeper into the minimap alignment and this is what I did and found:

After running emu, when I open the read-assignment-distribution.tsv and I sort (descending) the column of the marine worm (taxid: 121012), I see that "sequence3679" has a prob. of 1.0 for this taxon:

I then grep (grep "sequence3679" barcode33_emu_alignments.sam > minimap_alignment_read3679.sam) the alignments of this sequence from the alignments.sam file. You can find that result here:

minimap_alignment_read3679.txt

and only the alignment of the read with the marine worm (taxid 121012) (= grep "121012" minimap_alignment_read3679.sam > only_marineworm_read3679.sam).

only_marineworm_read3679.txt

I am not expert in reading .sam files, but do I understand correctly that:

1. The alignment with the marine worm - taxid 121012 - is actually very poor based on the "256" Flag, and the "0" in the 5th column? If yes, why is this then receives a prop. score of 1.0 in the emu-assignment-distributions.tsv?
2. The first alignment (taxid48693:nem_db:4805) receives a mapping quality (5th col) of "1" which is very poor.

What conclusions do you draw from this .sam output? I'm a bit puzzled why such poor assignments are included in the final emu output.

So, if you don't want the classification to show up via emu, then you will need to the minimap alignment restriction such that it does not return an alignment. I suggest you look at the minimap alignment and the minimap2 parameters to try to find which parameter adjustment would meet your needs.

is there a way to do this within emu? Or do you mean edit minimap code in Emu?
If not, then would this be a feature request because I cannot imagine I'm the only one who would want some alignment quality control. Else it seems impossible to me to actually know if your taxonomic assignments make any sense? Or do I miss something?

thanks for the help.

(Sorry for the open en close spam, something went wrong)

[kdcurry]

Minimap2 and emu use different metrics for defining which alignment each thinks is "better". Emu relies solely on the CIGAR strings (number of matches/mismatches/insertions, etc), and does not take read quality into account. It seems in this case, minimap2 favored the alignment to 48693, while emu favored 121012.

If you are concerned about low read quality reads being included in results, I suggest you run some sort of read quality filtering before running emu. Perhaps this would filter out this read.

[Robvh-git]

If you are concerned about low read quality reads being included in results, I suggest you run some sort of read quality filtering before running emu. Perhaps this would filter out this read.

@kdc10 sorry for the late reply, I was on holidays. We do filter on read quality (>Q20). I did not mean the quality of the reads (phred scores) but the quality of the alignment itself. In the sam file, the alignment with the marine worm is obviously very poor, yet emu still classifies it as being that taxon. Why is this taxon not labeled as "unassigned", as there is no real proper alignment?

In that sense, my question is, how can you know in the final emu output if the aligned taxa actually came from good reliable alignments or if it is mainly these kind of poor alignments?

You mention this in your previous comment:

So, if you don't want the classification to show up via emu, then you will need to the minimap alignment restriction such that it does not return an alignment. I suggest you look at the minimap alignment and the minimap2 parameters to try to find which parameter adjustment would meet your needs.

You mean I would have to edit the minimap command in the emu code?

[kdcurry]

Yes, currently you will need to edit the minimap command in the emu code. If there is a specific minimap parameter that you have success with, let us know we can update the emu code such that you pass the minimap parameter in directly from an emu command.

[Robvh-git]

Hi @kdcurry / @kdc10

I did some quick tests/checks and I think that there is no need to edit the minimap commands, but there should be a function in emu that would allow to filter the .sam files resulting from Minimap2.

Does emu now just selects the primary alignment of both strands of each read (so flag 0x0/0 and 16/0x10)?

How I currently see it, it is essential to include (at least) the following filtering steps:

1. filter om MAPQ (= 5th field in .sam file) to determine how "unique/probable" the alignment is. This could be done using samtools, e.g.: samtools view -q 20 minimap_out.sam > minimap_output_filtered.sam

2. filter on percentage identity: minimap2 provides the NM tag (e.g.: NM:i:545) that should show the edit distance (mismatches+indels) in number of bases between the input read and the reference read (if I'm not mistaken). By dividing that by the length of the query (SEQ column, 10th field), you should obtain the percentage ID.

Basically this would result in 2 additional emu arguments that allow you to set the MAPQ fitler value, and the perc ID you want to filter on.

I didn't validate this and these are thus suggestions/options. Even though I have experience in bio-informatics, I'm no "hardcore" bio-informatician/computer scientists, so I'm not comfortable enough to implement this in the emu code myself.

However, I think it is essential to include these steps inemu as now it is difficult to trust the emu results as you have literally no idea if the taxonomic identification is credible or not. You just get an output, regardless of its quality/confidence, and this thus seems dangerous.

[kdcurry]

Emu incorporates all primary and secondary alignments given by minimap2 to generate classification likelihoods, as described in the manuscript: https://www.nature.com/articles/s41592-022-01520-4.

I agree it would be ideal to have a more robust metric for marking reads and unclassified/unknown in emu. We do not have the bandwidth to develop and test this right now, but I welcome anyone from the community to do so.

[Robvh-git]

I agree it would be ideal to have a more robust metric for marking reads and unclassified/unknown in emu.

If I'm honest, I think that it is not ideal, but that is essential to produce reliable results.

We do not have the bandwidth to develop and test this right now, but I welcome anyone from the community to do so.

Unfortunate, but understandable. I hope that someone has the time to implement this soon.

[treangen]

Hi @Robvh-git, to follow up: I am reopening this issue as we are now looking into it and exploring potential solutions etc.