This seems like a very simple task to do and I'd assume that samtools view would be able to do this. Basically all I want to do is to return all reads that a entirely contained within a given region.
samtools view bamfile "chr1:100-1000" returns any read the overlaps that region. For example a read may start before 100 or extend beyond 1000. I only what those that map entirely inside that region.
•
written
by
george.wiggins •
1 answer
Use bedops and convert2bed (bam2bed) with a process substitution that specifies your region of interest:
$ bam2bed < reads.bam | bedops -e 100% - <(echo -e "chr1\t100\t1000") > answer.bed
Using bedops -e 100% asks for reads from read.bam that are contained entirely within chr1:100-1000.
Any reads that are found are written to answer.bed.
For RNA-seq data, you might wish to use the --split operand with bam2bed. See the documentation for more information.
References:
Log in to answer this question.
-
How to extract and analyze variant sequences from gnomAD v4 exome VCF
written by Rohan •
-
Extract aligned reads from a bam file for only selected coordinates
written by star715 •
-
samtools view region and include duplicate reads
written by mgene •
-
extracting mapped reads from the bam file
written by KVC_bioinfo
-
samtools return seek position
written by peiqi1122 •
-
Extract parts of reads from BAM file that overlap a specific region of genome
written by Chadi Saad
-
How can I extract sequence from a SAM file, with coordinates
written by SaltedPork
-
Reading CRAM files in pysam
written by QVINTVS_FABIVS_MAXIMVS
-
Getting read counts per loci
written by JameB11 •
-
Annotating RNA-Seq exon junctions
written by george.wiggins •
