Fastq reads mapped to the desired position on reference genome

Question

[Louis Kok]

Hi all, I would like to ask if there is a way to check if your reads are mapped to certain positions on the reference genome, given the regions to be checked is already defined. For example, if you expected reads to be mapped to genomic coordinates 2000 - 4000 bp, how do you check if certain reads are mapped to it? Thanks a lot.

Answer 1

[Dan D]

This is a possible duplicate. I'll let another mod make the call: How Do I Get The Reads That Correspond To A Particular Gene Of Interest

You can use the view command of samtools. For example:

samtools view [my.bam] chr:start-end

You can then do things like pipe the results into wc -l to get the total number of reads

Answer 2

[Alex Reynolds]

If your regions-of-interest (say, genes) are in BED format and you want to test all of them in one go, you could use bam2bed and bedmap --count to count the number of reads that overlap each region of interest:

$ bam2bed < reads.bam | bedmap --echo --count roi.bed - > answer.bed

[Dan D]

You can also feed the bed file into samtools:

samtools view -L [bedfile] [my.bed]

Answer 3

[xiaoleiusc]

For a customer genome, I use the below command:

samtools view [my.bam] genome_name:start-end | wc -l

To get genome_name, you can use samtools view -H [my.bam], the name is after the SN: