Hi,
I'm trying to extract reads based on their start coordinate in a bam file. I've tried using samtools view but that seems to give all reads covering that region, not originating there.
Apologies if this has been asked elsewhere but I couldn't find an answer via Google.
Many thanks.
•
written
by
fizoca •
Something like:
samtools view -h -q 10 in.bam chr1:1000-1500 | awk 'BEGIN{OFS="\t"}{if($1 ~ @) {print} else {if($4 >= 1000 || $4 <=1500) {print} else {}}}' | samtools view -Sbo output.bam -
One could restructure this to get around the default print action in awk, but this will also get the job done.
Thanks very much for the reply. I'm getting a couple of syntax errors when I try this if you wouldn't mind taking a look? I'm not great with awk.
awk: cmd. line:1: BEGIN{OFS="\t"}{if($1 ~ @) {print} else {if($4 >= 99795689 || $4 <=99795700) {print} else {}}}
awk: cmd. line:1: ^ syntax error
awk: cmd. line:1: BEGIN{OFS="\t"}{if($1 ~ @) {print} else {if($4 >= 99795689 || $4 <=99795700) {print} else {}}}
awk: cmd. line:1: ^ syntax error
$1 ~ @ should probably be something like $1 ~ /^"@"/.
I'm replying to an old thread, sorry.
While the awk solution is legitimate, samtools FILTER EXPRESSIONS are capable of filtering the starting position, so, if you want to get the reads that originate within the region chr1:1000-1500 you can use the following command:
# optional: append chr1:1000-1500 to the command to speed up the processing (the BAM need to be sorted and indexed)
samtools view -h -q 10 -e '1000 <= pos && pos <= 1500' in.bam
Log in to answer this question.
For me the question is not clear , I will assume that you have list of coordinate that you want to extract in such case it will be like this
samtools view -b -L ROI.bed file.bam > ROI_file.bamwhere region of interest ROI.bed will contains the start position you are interested in
otherwise give example to what you want to extract
Sorry I'll try to explain myself better. From my understanding if I run the command:
output.bam will contain all reads that fall between 1000-1500. I want to extract reads that start at any coordinate within this region. So for example a read that started at 1500 would also be extracted, even if most of the read does not lie within this region.
Apologies if I'm not correctly understanding how samtools view works!
clear now thanks :)
If your ultimate goal is to get the number of 5' position of reads at each genomic position, then you can try bedtools genomecov -5.
EDIT : it is not a very clear question indeed.