I am trying to select some reference genome region of a bam file, but got an error
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.
• 3,112 views
•
link
updated
by
lieven.sterck
•
written
by
schlogl
•
written
by
schlogl
0 answers
No answers yet.
Log in to answer this question.
no need to shout (== caps in the title, I changed it for you this time)
sorry. Thank you. Paulo
did you double check the correct syntax for querying a specific region?
Yep! Look here http://www.htslib.org/doc/samtools-view.html[1]
can you post the first line of the fasta file you're using? (= the header line)
AF086833.2 Ebola virus - Mayinga, Zaire, 1976, complete genome
as you figured out by now, you were not using the correct sequence name ;)
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
by
The result is:
AF086833.2=\=AF086833You got it bro. Thank you so much! 8)
see my post above, the "correct" name of your reference is AF086833.2 (note the .2)
still not worked even with the correct name! I don't get thet error, but also get a empty file!
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)
First thing to check if the chromosome name matches in all locations. If that checks out try
Hi bro...well! I tried that but didn't work too Genomax! 8(