Question: I am trying to select some reference genome region of a bam file, but got an error
0
gravatar for schlogl
5 weeks ago by
schlogl70
Brazil-Florianopolis
schlogl70 wrote:

Hi there, I am trying to use this command:

samtools view -b sorted_SRR1972739.bwa.bam AF086833:3050-3199 > selected.bam

And got this error:

[main_samview] region "AF086833:3050" specifies an unknown reference name. Continue anyway.

Using this command:

samtools view -c selected.bam
0

The file I used as ref genoma was AF086833 (AF086833.2 Ebola virus - Mayinga, Zaire, 1976, complete genome).

All the anterior commands worked well and were:

bwa mem AF086833.fa SRR1972739_1.fastq SRR1972739_2.fastq  > SRR1972739.bwa.sam
samtools view -S -b SRR1972739.bwa.sam > SRR1972739.bwa.bam
samtools sort SRR1972739.bwa.bam -o sorted_SRR1972739.bwa.bam
samtools index sorted_SRR1972739.bwa.bam

Do you have any ideia what I am doing wrong to get this error.

Thank you,

Paulo.

samtools • 158 views
ADD COMMENTlink modified 5 weeks ago by lieven.sterck8.5k • written 5 weeks ago by schlogl70
1

no need to shout (== caps in the title, I changed it for you this time)

ADD REPLYlink written 5 weeks ago by lieven.sterck8.5k

sorry. Thank you. Paulo

ADD REPLYlink written 5 weeks ago by schlogl70
1

did you double check the correct syntax for querying a specific region?

ADD REPLYlink written 5 weeks ago by lieven.sterck8.5k

Yep! Look here http://www.htslib.org/doc/samtools-view.html[1]

ADD REPLYlink modified 5 weeks ago by lieven.sterck8.5k • written 5 weeks ago by schlogl70

can you post the first line of the fasta file you're using? (= the header line)

ADD REPLYlink written 5 weeks ago by lieven.sterck8.5k

AF086833.2 Ebola virus - Mayinga, Zaire, 1976, complete genome

ADD REPLYlink written 5 weeks ago by schlogl70
1

as you figured out by now, you were not using the correct sequence name ;)

ADD REPLYlink written 5 weeks ago by lieven.sterck8.5k
1

What's the output of samtools idxstats sorted_SRR1972739.bwa.bam?

By the way, your command could be a lot shorter, avoiding unnecessary intermediate files. This assumes you are using a fairly recent version of samtools, but you should make sure to work with up-to-date tools anyway.

You can replace the following

bwa mem AF086833.fa SRR1972739_1.fastq SRR1972739_2.fastq  > SRR1972739.bwa.sam
samtools view -S -b SRR1972739.bwa.sam > SRR1972739.bwa.bam
samtools sort SRR1972739.bwa.bam -o sorted_SRR1972739.bwa.bam

by

bwa mem AF086833.fa SRR1972739_1.fastq SRR1972739_2.fastq  | samtools sort -o sorted_SRR1972739.bwa.bam
ADD REPLYlink written 5 weeks ago by WouterDeCoster44k

The result is:

samtools idxstats sorted_SRR1972739.bwa.bam

AF086833.2  18959   15279   11
*   0   0   5450
ADD REPLYlink written 5 weeks ago by schlogl70
2

AF086833.2 =\= AF086833

ADD REPLYlink written 5 weeks ago by genomax89k

You got it bro. Thank you so much! 8)

ADD REPLYlink written 5 weeks ago by schlogl70
1

see my post above, the "correct" name of your reference is AF086833.2 (note the .2)

ADD REPLYlink modified 5 weeks ago • written 5 weeks ago by lieven.sterck8.5k

still not worked even with the correct name! I don't get thet error, but also get a empty file!

ADD REPLYlink written 5 weeks ago by schlogl70
1

it's a fairly small region you're querying, perhaps there are just no reads mapped in that region?

try querying a bigger region (or one that you now for sure has reads mapped) to verify that your command is working (which it should as far as I can see)

ADD REPLYlink written 5 weeks ago by lieven.sterck8.5k
1

First thing to check if the chromosome name matches in all locations. If that checks out try

samtools view sorted_SRR1972739.bwa.bam "AF086833:3050-3199" > selected.bam
ADD REPLYlink written 5 weeks ago by genomax89k

Hi bro...well! I tried that but didn't work too Genomax! 8(

ADD REPLYlink written 5 weeks ago by schlogl70
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: 595 users visited in the last hour