Getting Reads from Bam file

[empyrean]

Hello every one.

I have a bam file which is approx 50gb from hiseq paired end data. I aligned it using BWA and now i am intersted in getting reads at particular location. i.e. for examle in chromosome 1 from (100000 - 101000). I would like to make a subset out from bam file with the above range. How can i do that ?

Thank you for your help!!

[iansealy]

Personally, I'd use samtools view to extract the region of interest into another BAM file. If you then want the reads then you could use Picard's SamToFastq.

Cheers,
Ian

[empyrean]

How can i specify the parameters for samtools view for extracting reads in that region?

I see the below options in samtools.

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
-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

[alpesh]

A region can be presented, for example, in the following format: ‘chr2’ (the whole chr2), ‘chr2:1000000’ (region starting from 1,000,000bp) or ‘chr2:1,000,000-2,000,000’ (region between 1,000,000 and 2,000,000bp including the end points). The coordinate is 1-based.

example

samtools view aln.sorted.bam chr2:20,100,000-20,200,000

[empyrean]

thank you.. it worked