Question: Failure in extracting reads from bam file due to unknown reference
0
gravatar for marongiu.luigi
11 days ago by
Germany, Mannheim, UMM
marongiu.luigi420 wrote:

Hello,

I need to extract reads from an alignment I have done of the whole human genome. I am looking for the gene MACC1 which is located at Chr 7 20134655..20217384, I used this command:

samtools view -b -h <AlnSrt.bam> "Chr7:20134655-20217384" > <subSet.bam>

but I got:

[main_samview] region "Chr7:20134655-20217384" specifies an unknown reference name. Continue anyway.

Using IGV I can see that the reference has a name chr7 so I changed to:

samtools view -b -h <AlnSrt.bam> "chr7:20134655-20217384" > <subSet.bam>

But in that case, nothing happens at all. The unknown reference name comes out even if I use 7:20134655-20217384.

How can I check the name of the reference?

How can I extract these reads?

Thank you

ADD COMMENTlink written 11 days ago by marongiu.luigi420
1

How can I check the name of the reference?

samtools view -H AlnSrt.bam | less
ADD REPLYlink modified 11 days ago • written 11 days ago by h.mon28k

So I got @SQ SN:chr7 LN:159345973 since the length is 159 345 973 and I am looking for 20 134 655, it should work with samtools view -b -h <alnsrt.bam> "chr7:20134655-20217384" > <subset.bam>

ADD REPLYlink written 11 days ago by marongiu.luigi420
1

Is it possible there is simply nothing aligned against that portion of the genome?

Try expanding the range, eg 1Mb either side.

chr7:20034655-20317384

If still nothing, go further. Also do other regions you pick work? If nothing works then maybe the index is broken.

ADD REPLYlink written 11 days ago by jkbonfield200

I am not getting anything even on the other chromosomes, but I can see the reads on IGV, which uses an index: enter image description here I can see that there are reads in this region; how can I check if the index is broken? Anyway, even with the error, a file has been created:

$ samtools view -c <subSet.bam>
9528

Would that be trustful?

Thank you

ADD REPLYlink modified 11 days ago • written 11 days ago by marongiu.luigi420

samtools view -b -h AlnSrt.bam chr7:20134655-20217384 > subSet.bam

But in that case, nothing happens at all.

When you say that nothing happens, do you mean that nothing appears on the terminal ? It is actually what is expected when the command works (there is no error or warning). This is probably when the subSet.bam file has been created. To make sure that it comes from this command and not from something else, delete subSet.bam and run the command again. By the way I don't think you need to quote (" ") the region in samtools view.

ADD REPLYlink written 11 days ago by Carlo Yague4.8k

Good point, I'll try that. Thank you -- Update: it worked.

ADD REPLYlink modified 10 days ago • written 11 days ago by marongiu.luigi420
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1516 users visited in the last hour