Hi,
I was wondering if there is an easy way to find gene density (by counting the number of genes for example in a 10-Mbp window) in R across the whole genome?
Is this normally done by downloading the transcript data bed file?
Thank you
Yes it's very easy in R.
First of all, install and load the Homo.sapiens package:
> biocLite('Homo.sapiens')
> library(Homo.sapiens)
> human.genes = genes(TxDb.Hsapiens.UCSC.hg19.knownGene)
> human.genes
GRanges object with 23056 ranges and 1 metadata column:
seqnames ranges strand | gene_id
<Rle> <IRanges> <Rle> | <character>
1 chr19 [ 58858172, 58874214] - | 1
10 chr8 [ 18248755, 18258723] + | 10
100 chr20 [ 43248163, 43280376] - | 100
1000 chr18 [ 25530930, 25757445] - | 1000
10000 chr1 [243651535, 244006886] - | 10000
... ... ... ... ... ...
9991 chr9 [114979995, 115095944] - | 9991
9992 chr21 [ 35736323, 35743440] + | 9992
9993 chr22 [ 19023795, 19109967] - | 9993
9994 chr6 [ 90539619, 90584155] + | 9994
9997 chr22 [ 50961997, 50964905] - | 9997
The human.genes object contains the coordinates of all human genes.
Then, create another GRanges object with the coordinates of all the sliding windows:
human.windows = tileGenome(seqinfo(human.genes), tilewidth=100000, cut.last.tile.in.chrom=T)
You can add a column with the number of overlapping genes by using countOverlaps:
> human.windows$totgenes = countOverlaps(human.windows, human.genes)
> human.windows
GRanges object with 31419 ranges and 1 metadata column:
seqnames ranges strand | totgenes
<Rle> <IRanges> <Rle> | <integer>
[1] chr1 [ 1, 100000] * | 3
[2] chr1 [100001, 200000] * | 1
[3] chr1 [200001, 300000] * | 0
[4] chr1 [300001, 400000] * | 0
[5] chr1 [400001, 500000] * | 0
... ... ... ... ... ...
[31415] chrUn_gl000245 [1, 36651] * | 0
[31416] chrUn_gl000246 [1, 38154] * | 0
[31417] chrUn_gl000247 [1, 36422] * | 0
[31418] chrUn_gl000248 [1, 39786] * | 0
[31419] chrUn_gl000249 [1, 38502] * | 0
Let's verify if this is correct. For example, the first window describes that the first 100,000 bases in chromosome 1 should contain three genes:
> subset(human.genes, seqnames=='chr1' & end < 100000)
GRanges object with 3 ranges and 1 metadata column:
seqnames ranges strand | gene_id
<Rle> <IRanges> <Rle> | <character>
100287102 chr1 [11874, 14409] + | 100287102
653635 chr1 [14362, 29961] - | 653635
79501 chr1 [69091, 70008] + | 79501
Let's take a window with more than 10 genes, and verify if this is correct:
In summary, you now have a human.windows object, where the totgenes column contains the count of genes per window. You can access it with human.windows$totgenes.
Note: you may want to merge all the overlapping genes into a single one. To do so, just use reduce(human.genes) at the beginning.
BEDOPS is pretty useful for solving this problem.
If you're not stuck using R, you could use something like this to get chromosome sizes, depending on your genome build:
$ mysql --user=genome --host=genome-mysql.cse.ucsc.edu -N -A -e "SELECT chrom,0,size FROM chromInfo;" hg19 | sort-bed - > hg19.chromSizes.bed
To run the window-to-gene mapping step, you will first need a sorted BED file called genes.bed. Assuming you are working with human hg19, you could get that via the following:
$ wget -O - ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_19/gencode.v19.annotation.gtf.gz \
| gunzip -c \
| grep -w "gene" \
| gtf2bed \
> gencode.v19.genes.bed
Once you have both chromosome sizes and genes in BED format, you could generate 10 Mb windows across each chromosome with bedops --chop and pipe the result to bedmap --count to get the number of genes over each window within each chromosome:
$ bedops --chop 10000000 hg19.chromSizes.bed | bedmap --echo --count - gencode.v19.genes.bed > hg19.10MbWindowsWithCountsOfOverlappingGencodeV19Genes.bed
All of this would depend on your reference genome or build, and what annotations you want to use to get genes, but the general commands and process are all identical.
This sounds like a job for bedtools.
Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
version 2.3.6
hello, sorry to bring this question up after so long, but do you know if there ways to do the same thing with my own genomic sequences and annotation gff file rather than using what is incorporated in bioconductor? Thanks a lot!
sure, you can use makeTxDbFromGFF from GenomicRanges to create the same TxDB object as above, then follow the same instructions.
hey, sorry to bother you again but i'm having problems when looking for genes that fall into each bin. to take a step further from finding genes in each bin, I was also trying to find genes with copy number variations(CNV) in each bin. I have already have a list of genes (with chromosome names, starting and ending positions for each gene) and converted into a GRange object using makeGRangefromdataframe function (plus i also added length of each chromosome at seqlength column), and then used countOverlap function for it and the window vector. it gave me number of CNV-affected genes in each bin, but when i tried to add # of cnv genes in all bins, the sum was larger than the total number of CNV-affected genes from my original files. Are you able to tell if I did those steps right, or why the total number of CNV-genes did not match? Thank you so much!