2
$\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?

For example get a sub BAM file which contains read line numbers [1, 3, 8, 70] of the original BAM.

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.

Edit:

  1. Size of the BAM file is ~7GB. Its human exome.

  2. I'm using Julia (https://github.com/BioJulia/XAM.jl/) to iterate over the BAM and it takes me ~300 seconds to iterate over all the records. I feel this is too slow.

  3. I'm running a function on each read of the BAM file and it decides which bucket the read has to go into. So this is specific subsampling.

I think samtools view - and awk route seems viable.

Improve this question
$\endgroup$
4
  • $\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__ ♦
    Jul 24 at 13:24
  • 1
    $\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
    Jul 24 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
    Jul 24 at 17:25
  • $\begingroup$ Julia is supposed to be fast (never used it). samtools will be faster. $\endgroup$
    – M__ ♦
    Jul 24 at 17:31

1 Answer 1

Reset to default
2
$\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.

Improve this answer
$\endgroup$

Your Answer

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

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