Calculating Gene Densities From Rnaseq Data

Question

[Ahmetz]

Hi,

Newbie here. I am interested plotting gene densities from my RNAseq data. I have a file with chromosome locations and read numbers for each gene. I want to plot each chromosome and get a bar graph with each bar indicating the number of genes present in 1MB bins. I guess I can count the number of genes in each 1MB region but I feel like there might be an R package out there that I'm not aware of. I'm planing to use qplot for the plotting. Any suggestions?

Thanks a lot!

Ahmet Z.

Answer 1

[Michael]

IRanges can compute the coverage given a set of intervals. That works for aligned RNA seq reads as for the gene bins itself. But what does RNAseq have to do with gene-density? Probably you meant something else?

for gene-density calculation in a 1MB window:

1. use the successiveIRanges function to create bins of the desires size.
2. Load the gene regions from a file, e.g with the rtracklayer package
3. use the function countOverlaps to count the number of genes in each bin

Can't give you more specific help unless you describe exactly what you want to do.

[Ahmetz]

Thanks for the info, I'll check out IRanges. I want to do a similar plot to Fig2 of this paper: http://genome.cshlp.org/content/early/2011/05/02/gr.116335.110.long Each chromosome is plotted with expressed gene densities and RNAseq would give you that information as to whether a gene is present or not. Would this help, in terms of specificity?

[Saima]

Hello, I tried to follow your suggestion for calculating gene density but I could not figure out how the successiveIRanges can be used for creating bins for all chromosomes, can you please elaborate your example?

Answer 2

[Neilfws]

A useful R package for plotting quantitative data in genomic context is GenomeGraphs (publication). It may be useful for your task if you can get the data into an appropriate form.

[Nibua]

nice package! Thank you!

Answer 3

[erdiazval]

I've been working on a similar approach to intersect gene density with cis-variants.

1) get .bam files (GATK's pipeline for SNP calling) and convert them to BED format with bedtools 2) create a genome.bed like-file (chr, length fields) and split it into 1MB chunks with bedtools 3) calculate coverage with the genomecov tool of bedtools 4) estimate genome coverage by 1MB steps for each chromosome with the bedmap tool of bedops (intersect the intervals in the *chunks.bed file against the coverage-per-position field in the coverage.bed file 5) plot this data in a circos plot!

Good luck