Reads clipping by refernce coordinates?

[shiningway]

I have a bam file and I’d like to retrieve the part of sequence of each reads that aligns to a certain region in the reference. For example, for reads that cover region chr1:1000-1100, I need to clip them such that only the part aligning to this region remains. The output should be reads that align exactly at this particular region, with all the flanking sequences clipped. Is there any tools available to do this? Thanks.

[TiborNagy]

You can combine two tools:
samtools view input.bam chr1:1000-1100 >slice.bam
bedtools bamtofastq -i slice.bam -fq output.fastq

[dpryan]

It's best to do something like this in python or another language, at least if you have any indels. The general idea is to extract the reads covering that regions (with samtools), extract the start position, parse the CIGAR to determine the position of each nucleotide in the read, output nucleotides in the position of interest.

@TiborNagy, that won't actually do what shiningway wants. If a read is AAACCTGACTGATCGAC and only the last 8 bases overlap the region, then the output for that read should be TGATCGAC.

[shiningway]

I was hoping to find existing tools to do this. But it seems I have to write it myself. Thanks.

It's best to do something like this in python or another language, at least if you have any indels. The general idea is to extract the reads covering that regions (with samtools), extract the start position, parse the CIGAR to determine the position of each nucleotide in the read, output nucleotides in the position of interest.

@TiborNagy, that won't actually do what shiningway wants. If a read is AAACCTGACTGATCGAC and only the last 8 bases overlap the region, then the output for that read should be TGATCGAC.