Different alignment from same minimap2 version

[airichli]

Hi,

Thanks for developing clair3-RNA. I am trying to repeat the performance evaluation in your paper using the raw reads. Although I use the same version of minimap2, same ref fasta, same raw reads as you use (https://s3.amazonaws.com/gtl-public-data/giab/bams/dRNA/03_30_23_R941_DRS_NA24143_dRNA_Guppy_6.4.6_rna_hac_prom.pass.NoU.fastq.gz.hg38.bam).

I found in my alignment, reads tend to have more mismatches, although every parameters are the same with you, it is interesting. I wonder if you could give me some advice on this matter?

test_data.zip

The attached is the my alignment and the alignment extract from the bam file above, and the IGV is as below:

As you can see in the figure, although with same parameters (as can be seen from bam header), mine has 3 mismatches and 1 deletion, yours have 2 deletions and 1 mismatch. More mismatches will lead to lower precision. Actually, this is not a random event, most of my alignments are not consistent with yours (alignment position is correct, just minor difference in the assignment of some mismatch and deletion), do you have any idea why there is such a difference? Thank you so much.

Best,
Aifu

[xianyu0623]

Hi Aifu,

Thanks for using Clair3-RNA!

I noticed that you're using a different thread setting, which might be the cause of the issue you're seeing. I recommend trying the same thread configuration to see if that resolves it.

If the problem persists, you may want to consider opening an issue in the minimap2 repository. My assumption is that since combinations like "3 mismatches + 1 deletion" and "2 deletions + 1 mismatch" could have the same alignment score, different machines might take slightly different actions when selecting the final result.

Let me know how it goes!

Best,
Xian.

[airichli]

Hi Xian,

I used 32 thread and did not solve the problem. To test your second hypothesis (different machines has different actions), could you please run my test data in your machine? It should be finished in about one second.

test_data_2.zip

Maybe I can add a bit more details about my problem, as shown in the IGV figure I posted before, for the read f665b36c-aada-4138-a846-4fb2abc72ec7, I used the same command as yours yet got the cigar (Cigar-1):
605S15M1D7M3D28M2I3M1I30M1I14M1D12M2D12M1I58M1D44M2I8M2D12M2D48M1D22M1D1M2D8M1I34M3I91M1I8M1I8M1I46M1D73M2I15M2D37M3D81M1D19M2I47M2D15M2D20M1D19M2D10M1I13M

The cigar (Cigar-2) in your alignment for this same read is:
605S15M1D7M3D28M2I3M1I30M1I14M1D12M2D12M1I50M1I3M1D4M1D44M2I8M2D12M2D48M1D22M1D1M2D8M1I34M3I91M1I8M1I8M1I46M1D73M2I15M2D37M3D81M1D19M2I47M2D15M2D20M1D19M2D10M1I13M

If you run my test data and still get the Cigar-2, then it should be the problem of the machine. Thanks for your assistance.

Anyone else could also try this if you have time, thanks.

This scenario does not just occur in this read, many other reads have similar issues, resulting in ~27% precision for dRNA002 HG004 data in my result, which is far below the 66% precision mentioned in your manuscript.

[xianyu0623]

Hi Aifu,

Yes, I got the same alignment results as you.

However, I did notice something that might explain the discrepancies in alignment results—where did you get the FASTQ files? Did you convert them from the aligned BAM files using samtools fastq? I ask because your sequences seem to match those in the BAM file, which made me curious.

In our case, we didn’t run the alignment ourselves. We're working directly with the aligned BAM files, specifically the merged set for HG004. (https://s3.amazonaws.com/gtl-public-data/giab/bams/dRNA/03_30_23_R941_DRS_NA24143_dRNA_Guppy_6.4.6_rna_hac_prom.pass.NoU.fastq.gz.hg38.bam)

Cheers,
Xian.

[airichli]

Hi Xian,

Thanks for your timely reply. The fastq could be extracted from the bam file.
I think the alignment step may affect the performance of clair3-RNA. I tried minimap2 versions from 2.17 to 2.29, all output the same alignment cigar, which is different with the cigar in the bam file fro GIAB (e.g., for the read mentioned here ). It is interesting that GIAB generated bam files tend to have fewer mismatch and more indels, which may result in higher precision in your evaluation. Actually with my alignment, the cDNA R9.4 data have only 27% precision.
I would like to run clair3-rna on the cDNA R10.4 version data and see its performance, could you please tell me where do you download it? thanks.

[airichli]

I have found R10 cDNA data. Thanks.