Hi,
I am using pysam in python to extract mate reads of a list of target alignments, but get error "ValueError: mate MG01HX01:920:H3YVTCCX2:8:1223:25530:16850: is unmapped". I am just wondering whether pysam can get mate reads for both mapped and unmapped reads? Below is my code. I really appreaciated it if you can give adice for this issue. Thank you!
No, pysam cannot get the mates when the mate is unmapped read.
The reason for this is that in order to fetch the mate, it must know where the mate is in the BAM file. It does this using the index file. The a bam index tells the parser the position of the reads mapping to each part of each chromosome start in the BAM file as a number of bytes from the start of the file.
As an unmapped read doesn't have a position, the parser cannot find offset to look for it at. In theory the BAM specification could specify that unmapped in a contig, and give its position in the index. Unfortunately it doesn't.
interesting point, never thought about it that way,
what the OP could do is build a dictionary with unmapped reads (hopefully these are few and far between) store the unmapped alignments there and when the mate is unmapped select it from that storage.
No, not directly.
There are two ways to do this.
until_eof=True because otherwise fetch doesn't return unmapped reads. Then test every read for read.is_unmapped.Neither is particularly satisfactory.
Not also that even where you do have properly mapped pairs, AlignmentFile.mate is slow. It is specifically noted in the documentation that it is not suitable for high-throughput work.
You can get all unmapped read pairs efficiently by running AlignmentFile.fetch('*').
see: https://github.com/pysam-developers/pysam/issues/424#issuecomment-2192755497
AFAIK fetch('*') only returns reads where both mates are unmapped.
When a read pair has 1 mate that's mapped and one that's unmapped, the unmapped mate will be located next to their mapped mate in the file, so should be returned by the regular interval query (ie. fetch('chr1:12345-54321')). In contrast, the unmapped read pairs are located at the end of the BAM/CRAM file, and so require a special query.
Log in to answer this question.
please advice.
I think you should check your data to ensure that you do have the mate in the BAM file.
And only then troubleshoot the code. Make sure that
MG01HX01:920:H3YVTCCX2:8:1223:25530:16850can be found twice in the datanow check the flags etc.