Finding gene density from reference genome using R
3
0
Entering edit mode
6.0 years ago

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?

Thank you

R biomaRt • 5.8k views
5
Entering edit mode
6.0 years ago

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.

0
Entering edit mode

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!

0
Entering edit mode

sure, you can use makeTxDbFromGFF from GenomicRanges to create the same TxDB object as above, then follow the same instructions.

0
Entering edit mode

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!

4
Entering edit mode
6.0 years ago

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.

0
Entering edit mode
6.0 years ago

This sounds like a job for bedtools.