How can I extract the reads that map in genomic region?

[HeyIamNuria]

Hello,
I have a bam file with the alignment of Illumina reads to a genome, and a bed file with the coordinates of some genomic regions.

I found out the coverage of these regions using Bedtools, and now I want to retrieve the reads that map in these regions, but I do not know how to get the names of the reads.

Is there any available tool to extract the reads that map in a particular genomic region?

I thought Bedtools or Samtools would have it but I cannot find it.
My only other idea is to blast the reads against these regions, but I know that the results would be different.

Thank you for your help

[dpryan]

Assuming your BAM file is sorted and indexed:

samtools view -h -L Regions.bed alignments.bam > alignments_in_regions.sam

[dexterslab]

Assuming your BAM file is sorted and indexed:

samtools view -h -L Regions.bed alignments.bam > alignments_in_regions.sam

Sorry for blatantly hijacking this thread with a follow up question:

Assuming paired-end reads, would this suggested command also extract reads wherein only one of the mates falls into the region? If not, how could I archive that?

[dpryan]

Sorry for blatantly hijacking this thread with a follow up question:

Assuming paired-end reads, would this suggested command also extract reads wherein only one of the mates falls into the region? If not, how could I archive that?

Good question. Looking at the source code for samtools, it appears that it processes reads one at a time (i.e., not as pairs). So, if only one mate in a pair overlapped a region, I expect only that read would be output and not its mate. If you need both mates, it's possible that bedops implements that (I've never looked, but wouldn't be surprised).

[HeyIamNuria]

Devon,
maybe I'm wrong but did you meant to type "-L"?
Because:

"view: invalid option -- 'L'

Usage: samtools view [options] <in.bam>|<in.sam> [region1 [...]]

Options: -b output BAM
-h print header for the SAM output
-H print header only (no alignments)
-S input is SAM
-u uncompressed BAM output (force -b)
-x output FLAG in HEX (samtools-C specific)
-X output FLAG in string (samtools-C specific)
-t FILE list of reference names and lengths (force -S) [null]
-T FILE reference sequence file (force -S) [null]
-o FILE output file name [stdout]
-f INT required flag, 0 for unset [0]
-F INT filtering flag, 0 for unset [0]
-q INT minimum mapping quality [0]
-l STR only output reads in library STR [null]
-r STR only output reads in read group STR [null]
-? longer help"

[dpryan]

Yes, but it seems that you have an earlier version of samtools (0.1.17 maybe?). Try upgrading, as I know the -L option is available in 0.1.18 and 0.1.19 (the most recent).

[dpryan]

BTW, here are the options for the most recent version (0.1.19):

Usage:   samtools view [options] <in.bam>|<in.sam> [region1 [...]]

Options: -b       output BAM
         -h       print header for the SAM output
         -H       print header only (no alignments)
         -S       input is SAM
         -u       uncompressed BAM output (force -b)
         -1       fast compression (force -b)
         -x       output FLAG in HEX (samtools-C specific)
         -X       output FLAG in string (samtools-C specific)
         -c       print only the count of matching records
         -B       collapse the backward CIGAR operation
         -@ INT   number of BAM compression threads [0]
         -L FILE  output alignments overlapping the input BED FILE [null]
         -t FILE  list of reference names and lengths (force -S) [null]
         -T FILE  reference sequence file (force -S) [null]
         -o FILE  output file name [stdout]
         -R FILE  list of read groups to be outputted [null]
         -f INT   required flag, 0 for unset [0]
         -F INT   filtering flag, 0 for unset [0]
         -q INT   minimum mapping quality [0]
         -l STR   only output reads in library STR [null]
         -r STR   only output reads in read group STR [null]
         -s FLOAT fraction of templates to subsample; integer part as seed [-1]
         -?       longer help

You'll also often find the -@ option useful, at least if you ever need to make BAM files.