I have a bam file and I wanted to extract the sequence of a specific region. for example ch1:1-337, the sequence starts EXACTLY at position 1 and ends EXACTLY at position 337 and I don't want anything else that starts and ends between them.
I tried:
samtools view -b input.bam ch1:1-337 > output.bam
and:
samtools view -b -L regions.bed input.bam > output.bam
but the result was every sequence that starts in the region or ends in it.
I can extract the sequence from igv but that will take too much time to do for multiple regions (> 150 regions).
Any ideas?
•
written
by
malkomalktho •
samtools view -b --expr 'pos==1 && endpos==337' input.bam "chr1:1:337" > output.bam
•
written
by
Pierre Lindenbaum
thank you for the answer, but i got an empty bam file. i also replaced "ch1:1-337" with "chr1:1-337" it didn't work either.
•
written
by
malkomalktho •
i also replaced "ch1:1-337" with "chr1:1-337"
ah yes, fixed.
but i got an empty bam file
are you sure such read exists ? can you please show us one alignment ?
the read did exist. i extracted all reads using igv manually. thank you though. for example:
HWI-EAS746_0007:7:46:8030:13231#0|GAT 16 chr1 29304 2 27M *0 0 CCCGGACGCCTCCTCACCTGCGAGCCC ca_da^aaaa\aaad^ecc^effafff NM:i:1 MD:Z:26G0 XP:Z:~88F~~~~~~~~~~~~~~~~~~~~~~~
link to bam file (example): http://hgdownload-test.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeRikenCage/releaseLatest/wgEncodeRikenCageA549CellPapAlnRep1.bam
this read starts at position 29304 not at 1 as you asked.
the 1 and 337 were just examples
Log in to answer this question.
Perhaps this is what you are looking for: Extract from a .bam file the reads that map to specific area in the reference genome