3
$\begingroup$

Consider a sorted, indexed, only mapped reads containing BAM file.

Is there a way to get a sub BAM based on read line numbers?

This can be done iteratively using a counter but its too slow. Is there a way to do this using samtools view or some tool like it.

$\endgroup$
5
  • $\begingroup$ Welcome to the site . It would be helpful to state the size of the bam file, and particularly how you did this "reiteratively" (code if possible) - is essential because samtools is not the only solution. Like did you use Python and some sort of re, which would explain it? Also just run time, what is slow? Finally any additional information, e.g. it's human always helps, e.g. if it turns out to be bacteria ... thats different. $\endgroup$
    – M__
    Commented Jul 24, 2023 at 13:24
  • 2
    $\begingroup$ Please elaborate the purpose. How would you know that something in line x,y,z, is interesting so you want to subset for it? What is the context? Is it random subsampling or specific subsampling? In the simplest case you print the bam as sam via samtools view -h, skip header lines, and then make sed or awk print only lines overlapping an array of numbers you define upfront. That's almost no bottleneck, it just has to go through the entire file once which is not slow. $\endgroup$
    – ATpoint
    Commented Jul 24, 2023 at 14:32
  • 1
    $\begingroup$ I've added more information in a edit. Can you please give me a sample code on how I can use AWK to print only lines overlapping an array of numbers. I know to print Nth line you could do awk 'NR==N' $\endgroup$
    – papabiceps
    Commented Jul 24, 2023 at 17:25
  • $\begingroup$ Julia is supposed to be fast (never used it). samtools will be faster. $\endgroup$
    – M__
    Commented Jul 24, 2023 at 17:31
  • $\begingroup$ I second @ATpoint comment: selecting arbitrary lines from position sorted BAM is puzzling. If you need a random fraction of the mapped reads, there is samtools view --subsample 0.1234 my.bam. Selecting reads based on read/pair IDs can be done in one step. $\endgroup$
    – darked89
    Commented May 16 at 8:51

1 Answer 1

4
$\begingroup$

Samtools+awk:

samtools view -h aln.bam \
  | awk 'BEGIN{while((getline<"lineno.txt")>0)a[$1]=1}/^@/{++n;print}!/^@/&&a[NR-n]' \
  | samtools view -b --threads=4 > sub.bam

I would recommend mawk which is much faster than other awk.

I believe using Julia should be faster as it would not need to convert BAM to SAM and then convert back. Five min to read through a 7GB BAM file sounds about right. That is not slow.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.