Is there a way to fetch a read efficiently from a BAM file, using Pysam or a similar module (from Python), by its read ID?
For example, if I have a list of read IDs, "read_ids", I want to do something like:
bam_file = pysam.Samfile(bam_filename, "rb")
for read_id in read_ids:
# fetch the read id?
my_aligned_read = bam_file.fetch(read_id)
is there a way to do this? The indexed/sorted BAM format should have all this information I am just wondering how to retrieve it.
thanks.
1 answer
the short answer is "no", not without some programming on your own. The BAM index is by location, not by read-id. You will have to create your own index to access by name.
I had written some code to do this using tokyo cabinet to save an index of sam header to file position (from which you can then read the SAM info). That code is here. (if you dont like the tokyo cabinet dependency, Istvan Albert pointed out that it's just as well to use the bsddb module that comes with python up to at least 2.6.)
I believe you can write your own index using screed as well--by default it supports fasta and fastq formats--but I have not tried that. It uses an sqlite backend.
What about the -n option of samtools sort?
-n Sort by read names rather than by chromosomal coordinates"
Log in to answer this question.
-
Expect Library - Bam-Readcount
written by tayebwajb
-
Sam To Bam - Loss Of Data Or Just Great Compression?
written by Clare
-
Samtools View Of A Sam File Entire Chromosome
written by dfernan
-
Adding Custom Bed And Bam Tracks To UCSC Browser By URL
written by user
-
B/Sam To Fasta In R
written by laemtao •
-
Emule "Samtools View" Using Pysam
written by Geparada
-
Filtering Stream Bam Output Idealy From Pysam/Python
written by User 9996
-
Getting Unmapped Reads: Comparing Fastq To Bam
written by User 9996
-
How To Read Bam.Ebd, Bam.Ebi Format Dataset
written by Lerong •
-
How To Get Number Of Reads In Bam File Efficiently In Python?
written by User 9996
I am hitting the same issue ..just wondering what solution worked out for you ??