Question: Failure in extracting reads from bam file due to unknown reference
0
gravatar for marongiu.luigi
9 months ago by
Germany, Mannheim, UMM
marongiu.luigi520 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 9 months ago by marongiu.luigi520
1

How can I check the name of the reference?

samtools view -H AlnSrt.bam | less
ADD REPLYlink modified 9 months ago • written 9 months ago by h.mon30k

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 9 months ago by marongiu.luigi520
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 9 months ago by jkbonfield400

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 9 months ago • written 9 months ago by marongiu.luigi520

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 9 months ago by Carlo Yague5.0k

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

ADD REPLYlink modified 9 months ago • written 9 months ago by marongiu.luigi520
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: 1525 users visited in the last hour