Read coverage for specific region of the genome

[shi]

I have RNA-seq data from 14 clones producing my protein of interest. I used Tophat2 to align my reads to the genome. In order to check for ER stress I would like to see whether a specific 26bp intron in XBP-1 is removed or not.

To solve this my plan is to calculate the read coverage of bp position 35303-35329 in scaffold NW_003614142.1 and check whether this is lower than the surrounding regions = bp region 35250-35302 and 35329 -35359.

Is there an easy way to do this? The information should be pretty straight forward to get from the accepted_hits.bam file. As far as i can tell samtools mpileup (http://samtools.sourceforge.net/samtools.shtml) should be able to do this but i can't get it to work (my programming stills are very limited). Any input to how you would extract the information or how i should correct the mpileup idea are more than welcome

samtools mpileup -b accepted_hits.bam -r NW_003614142.1 > example

You can do this using the featureCounts program included in Subread package (http://subread.sourceforge.net).

# create a tab-delimited annotation file ("annot.txt")
GeneID	Chr	Start	End	Strand
1	NW_003614142.1	35303	35329	+
2	NW_003614142.1	35250	35302	+
3	NW_003614142.1	35329	35359	+

# count reads overlapping with each region
featureCounts -a annot.txt -i your_reads.bam -o counts.txt -F 'SAF' -b -f -O

[mikep]

Try samtools view with the -c parameter, which will generate a count summary of your regions.

[Kaas]

Hi

Thank you for the answer. That worked fine... but the results did not match what i found in CLC genomic workbench. Seems that it can only be used for DNA and not RNA data with introns as it probably only use the information in the bam file of start of annealing and length of annealing. Results can be seen here, but the intron has the same coverage as the flanking regions using bedtools.

[dariober]

I have RNA-seq data from 14 clones producing my protein of interest. I used Tophat2 to align my reads to the genome. In order to check for ER stress I would like to see whether a specific 26bp intron in XBP-1 is removed or not.

To solve this my plan is to calculate the read coverage of bp position 35303-35329 in scaffold NW_003614142.1 and check whether this is lower than the surrounding regions = bp region 35250-35302 and 35329 -35359.

Is there an easy way to do this? The information should be pretty straight forward to get from the accepted_hits.bam file. As far as i can tell samtools mpileup (http://samtools.sourceforge.net/samtools.shtml) should be able to do this but i can't get it to work (my programming stills are very limited). Any input to how you would extract the information or how i should correct the mpileup idea are more than welcome

samtools mpileup -b accepted_hits.bam -r NW_003614142.1 > example

Hi- You can use coverageBed in bedtools http://bedtools.readthedocs.org/en/latest/:

## Bed file of intervals where to count reads (i.e. your intron):
cat b.bed
NW_003614142.1 35303 35329
NW_003614142.1 35250-35302
NW_003614142.1 35329 -35359
...

## Count reads in bed intervals
coverageBed -abam accepted_hits.bam -b b.bed

[Kaas]

I have RNA-seq data from 14 clones producing my protein of interest. I used Tophat2 to align my reads to the genome. In order to check for ER stress I would like to see whether a specific 26bp intron in XBP-1 is removed or not.

To solve this my plan is to calculate the read coverage of bp position 35303-35329 in scaffold NW_003614142.1 and check whether this is lower than the surrounding regions = bp region 35250-35302 and 35329 -35359.

Is there an easy way to do this? The information should be pretty straight forward to get from the accepted_hits.bam file. As far as i can tell samtools mpileup (http://samtools.sourceforge.net/samtools.shtml) should be able to do this but i can't get it to work (my programming stills are very limited). Any input to how you would extract the information or how i should correct the mpileup idea are more than welcome

samtools mpileup -b accepted_hits.bam -r NW_003614142.1 > example