Subsetting BAM by positions
Hi All,
I have a .cram file and I want to keep only reads that are related to positions I am interested in.
Say the list of positions looks like this (it has about 16k positions):
rs3094315 1 0.020130 752566 G A
rs12124819 1 0.020242 776546 A G
rs28765502 1 0.022137 832918 T C
rs7419119 1 0.022518 842013 T G
rs950122 1 0.022720 846864 G C
rs113171913 1 0.023436 869303 C T
rs13302957 1 0.024116 891021 G A
rs59986066 1 0.024183 893462 C T
rs112905931 1 0.024260 896271 C T
rs6696609 1 0.024457 903426 C T
rs13303368 1 0.024771 914852 G C
I tried using samtools view and the region attributes like this:
samtools view -T reference_genome/hs38DH.chr22.fa.gz NA12878.final.cram chr22:59030570-59030570 chr22:28799103-28799103 chr22:28799108-28799108
and it works great, but when I put the 16k positions it says its too long.
-bash: /usr/bin/samtools: Argument list too long
How do I keep only reads that have a position that appears in my position list?
1 answer
Log in to answer this question.
Assuming that this is not a limitation on
samtoolsside you may be running into limitation on stack space for command line arguments. See if the solution here helps. You could also loop through the positions perhaps?Thank you GenoMax for your super quick reply.
I thought about looping but wanted to avoid it. It seems like a pretty resonable request so I figured there must be a better solution... I'll look into the link you added.
Thanks! Batel.
What GenoMax says here is correct, you're running in the limit of your linux(?) operating system. Looping might be your only option here (or Xargs ?)
You don't have to loop over each of those regions, you can 'reduce' the loop by taking a number of ranges together in the loop