1
$\begingroup$

I've got some data on genomic regions of interest in a bedgraph-like format from the Mus musculus [mm10] genome, and would like to find out the nearest gene to these regions.

Some of the regions will reside within a gene (in which case I'd like to know that they're within a gene), some regions may reside within multiple genes (in which case I'd like to know all gene names), and some regions are not within a gene (in which case I'd like to know the closest gene: the shortest sequence gap between either end of the gene and either end of the tagged genomic region. In this case I want to do an unstranded search, because it's possible that I might have regions describing suppressing transcripts. Here's some example data:

Ref         Start             End      fwd     rev
chrX     73716152        73716152        0      -2
chrX    167207094       167209162        0      -3
chrY     30770844        30772724        3       0
chr1     24613189        24613641        0       3
chr1     24613971        24614812        0       4
chr1     35788729        35789387       -5       0
chr1     93151586        93160225        0      -2
chr1    107511454       107525597       -2       0
chr1    130882067       130887413        0      -3
chr1    134182466       134189942        2       0
chr1    135818954       135832956       -2       0
chr1    149520498       149521157        0      -3
chr1    166166500       166166500       -2       0
chr2      5379059         5379374        0      -2
chr2     32379110        32381837        0      -2
chr2     32381854        32381858        0      -2
chr2     32384633        32384633        0      -2
chr2     32384636        32387773        0      -3
chr2    120539390       120563799        0      -2
chr2    125946785       125975527        2       0
chr2    127633311       127633684        0      -2
chr2    127634841       127656210        0      -3
chr2    164562717       164562717       -2       0
chr2    164562719       164562720       -2       0
chr2    164568498       164568499       -2       0
chr2    181188007       181188388       -2       0

I'd like to add another 'Gene' column to this table that gives the official names of the nearest gene(s) for each region defined by <Ref>:<Start>-<End>.

$\endgroup$
1
  • $\begingroup$ What would be the desired output format? Would you like the name of the gene in that new column? Are you interested in the nearest gene in the same strand or it doesn't matter? What happen if multiple genes are in those regions? $\endgroup$
    – llrs
    Jan 16 '18 at 21:09
4
$\begingroup$

You're looking for bedtools closest, which thankfully defaults both considering overlapping genes and outputting all hits in case of ties.

  bedtools closest -a foo.txt -b genes.bed | awk '{foo = sprintf("%s:%i-%i", $6, $7, $8); OFS="\t"; print $1, $2, $3, $4, $5, foo}'

You'll need to strip the header and ensure the file is sorted. I've piped the output to awk to reformat things. The output then is something like:

chr1    24613189    24613641    0   3   chr1:24613188-24613971
chr1    24613971    24614812    0   4   chr1:24613973-24614651
chr1    35788729    35789387    -5  0   chr1:35806973-35810462
chr1    93151586    93160225    0   -2  chr1:93151348-93160948
chr1    93151586    93160225    0   -2  chr1:93151582-93160870
chr1    107511454   107525597   -2  0   chr1:107511422-107525598

Those results are for the Gencode M15 release (the BED file I used has transcripts rather than genes, so you'll get a bit different coordinates). As an aside, the gene name would normally be output as well, just add column 9 to the output.

$\endgroup$
0
1
$\begingroup$

With the prompting from @Devon_Ryan's answer, I came up with a solution that worked well for me. I downloaded the mm10 RefSeq gene locations from UCSC [Fields: chrom, txStart, txEnd, name2], and sorted both input files using bedtools sort. From there the bedtools closest command got me almost to where I wanted to go:

(echo -e "Chr\tStart\tEnd\tfwd\trev\tm_Chr\tm_Start\tm_End\tGene\tDistance"; 
  bedtools closest -d -a regions.bed -b NCBI_RefSeq_mm10.bed) > annotated.tsv

Output example:

Chr     Start   End     fwd     rev     m_Chr   m_Start m_End   Gene    Distance
chrX    73716152        73716152        0       -2      chrX    73686177        73716204        Bcap31  0
chrX    73716152        73716152        0       -2      chrX    73686177        73717815        Bcap31  0
chrX    73716152        73716152        0       -2      chrX    73686177        73716401        Bcap31  0
chrX    73716152        73716152        0       -2      chrX    73686177        73716502        Bcap31  0
chrX    167207094       167209162       0       -3      chrX    167207092       167213479       Tmsb4x  0
chrX    167207094       167209162       0       -3      chrX    167207093       167209218       Tmsb4x  0
chrY    30770844        30772724        3       0       chrY    30725674        30738546        Gm29582 32299
chr1    24613189        24613641        0       3       chr1    24257682        24587535        Col19a1 25655
chr1    24613189        24613641        0       3       chr1    24257682        24587535        Col19a1 25655
chr1    24613189        24613641        0       3       chr1    24257682        24587535        Col19a1 25655
chr1    24613189        24613641        0       3       chr1    24257682        24587535        Col19a1 25655
chr1    24613971        24614812        0       4       chr1    24257682        24587535        Col19a1 26437
chr1    24613971        24614812        0       4       chr1    24257682        24587535        Col19a1 26437
chr1    24613971        24614812        0       4       chr1    24257682        24587535        Col19a1 26437
chr1    24613971        24614812        0       4       chr1    24257682        24587535        Col19a1 26437
chr1    35788729        35789387        -5      0       chr1    35869594        35881164        1110002O04Rik   80208
chr1    93151586        93160225        0       -2      chr1    93151353        93160935        2310007B03Rik   0
chr1    93151586        93160225        0       -2      chr1    93151353        93160957        2310007B03Rik   0
chr1    93151586        93160225        0       -2      chr1    93151354        93160948        2310007B03Rik   0
chr1    93151586        93160225        0       -2      chr1    93151354        93160870        2310007B03Rik   0
chr1    107511454       107525597       -2      0       chr1    107500891       107525600       Serpinb2        0
chr1    107511454       107525597       -2      0       chr1    107511422       107525600       Serpinb2        0
chr1    107511454       107525597       -2      0       chr1    107511509       107525600       Serpinb2        0
chr1    130882067       130887413       0       -3      chr1    130882073       130887411       Il24    0
chr1    134182466       134189942       2       0       chr1    134182169       134190031       Chil1   0
chr1    134182466       134189942       2       0       chr1    134182403       134190031       Chil1   0
chr1    134182466       134189942       2       0       chr1    134182408       134190031       Chil1   0
chr1    135818954       135832956       -2      0       chr1    135818597       135833341       Lad1    0
chr1    149520498       149521157       0       -3      chr1    149410728       149449931       Gm29398 70568
chr1    166166500       166166500       -2      0       chr1    166130459       166166510       Gpa33   0
chr2    5379059         5379374         0       -2      chr2    5293456         5676046         Camk1d  0
chr2    5379059         5379374         0       -2      chr2    5293456         5714762         Camk1d  0
chr2    5379059         5379374         0       -2      chr2    5293456         5714762         Camk1d  0
chr2    5379059         5379374         0       -2      chr2    5293456         5714762         Camk1d  0
chr2    5379059         5379374         0       -2      chr2    5293458         5383533         Camk1d  0
chr2    5379059         5379374         0       -2      chr2    5293458         5569754         Camk1d  0
chr2    5379059         5379374         0       -2      chr2    5293458         5601302         Camk1d  0
chr2    5379059         5379374         0       -2      chr2    5293458         5715350         Camk1d  0
chr2    5379059         5379374         0       -2      chr2    5298992         5715357         Camk1d  0
chr2    32379110        32381837        0       -2      chr2    32379100        32381915        1110008P14Rik   0
chr2    32381854        32381858        0       -2      chr2    32379100        32381915        1110008P14Rik   0
chr2    32384633        32384633        0       -2      chr2    32384636        32387739        Lcn2    3
chr2    32384636        32387773        0       -3      chr2    32384636        32387739        Lcn2    0
chr2    120539390       120563799       0       -2      chr2    120506829       120563831       Zfp106  0
chr2    120539390       120563799       0       -2      chr2    120507433       120563851       Zfp106  0
chr2    120539390       120563799       0       -2      chr2    120507433       120563853       Zfp106  0
chr2    125946785       125975527       2       0       chr2    125859108       125984298       Galk2   0
chr2    125946785       125975527       2       0       chr2    125859228       125983587       Galk2   0
chr2    125946785       125975527       2       0       chr2    125866111       125983587       Galk2   0
chr2    125946785       125975527       2       0       chr2    125866218       125983587       Galk2   0
chr2    125946785       125975527       2       0       chr2    125866219       125983587       Galk2   0
chr2    125946785       125975527       2       0       chr2    125866222       125984298       Galk2   0
chr2    125946785       125975527       2       0       chr2    125920565       125983587       Galk2   0
chr2    127633311       127633684       0       -2      chr2    127633225       127656695       Mal     0
chr2    127633311       127633684       0       -2      chr2    127633225       127656695       Mal     0
chr2    127634841       127656210       0       -3      chr2    127633225       127656695       Mal     0
chr2    127634841       127656210       0       -3      chr2    127633225       127656695       Mal     0
chr2    164562717       164562717       -2      0       chr2    164562560       164568510       Wfdc2   0
chr2    164562717       164562717       -2      0       chr2    164562715       164568506       Wfdc2   0
chr2    164562719       164562720       -2      0       chr2    164562560       164568510       Wfdc2   0
chr2    164562719       164562720       -2      0       chr2    164562715       164568506       Wfdc2   0
chr2    164568498       164568499       -2      0       chr2    164562560       164568510       Wfdc2   0
chr2    164568498       164568499       -2      0       chr2    164562715       164568506       Wfdc2   0
chr2    181188007       181188388       -2      0       chr2    181187221       181188506       Ppdpf   0
chr2    181188007       181188388       -2      0       chr2    181187342       181188504       Ppdpf   0

After that it was just a matter of aggregating the gene names for each unique region with a quick R script:

library(dplyr);
data.tbl <-
    read.delim("annotated.tsv") %>%
    group_by(Chr, Start, End, fwd, rev) %>%
    summarise(Genes = paste(unique(Gene), collapse=","), Distance=min(Distance));

write.table(data.tbl, "collapsed_annotated.tsv",
            sep="\t", row.names=FALSE, quote=FALSE);

Collapsed output example:

Chr     Start   End     fwd     rev     Genes   Distance
chr1    24613189        24613641        0       3       Col19a1 25655
chr1    24613971        24614812        0       4       Col19a1 26437
chr1    35788729        35789387        -5      0       1110002O04Rik   80208
chr1    93151586        93160225        0       -2      2310007B03Rik   0
chr1    107511454       107525597       -2      0       Serpinb2        0
chr1    130882067       130887413       0       -3      Il24    0
chr1    134182466       134189942       2       0       Chil1   0
chr1    135818954       135832956       -2      0       Lad1    0
chr1    149520498       149521157       0       -3      Gm29398 70568
chr1    166166500       166166500       -2      0       Gpa33   0
chr2    5379059         5379374         0       -2      Camk1d  0
chr2    32379110        32381837        0       -2      1110008P14Rik   0
chr2    32381854        32381858        0       -2      1110008P14Rik   0
chr2    32384633        32384633        0       -2      Lcn2    3
chr2    32384636        32387773        0       -3      Lcn2    0
chr2    120539390       120563799       0       -2      Zfp106  0
chr2    125946785       125975527       2       0       Galk2   0
chr2    127633311       127633684       0       -2      Mal     0
chr2    127634841       127656210       0       -3      Mal     0
chr2    164562717       164562717       -2      0       Wfdc2   0
chr2    164562719       164562720       -2      0       Wfdc2   0
chr2    164568498       164568499       -2      0       Wfdc2   0
chr2    181188007       181188388       -2      0       Ppdpf   0
chrX    73716152        73716152        0       -2      Bcap31  0
chrX    167207094       167209162       0       -3      Tmsb4x  0
chrY    30770844        30772724        3       0       Gm29582 32299
$\endgroup$
1
$\begingroup$

Another option is BEDOPS closest-features, which efficiently locates nearest annotations, overlapping or otherwise.

$ closest-features --closest --delim '\t' intervals.bed genes.bed > answer.bed

It's easy to reformat output by piping to awk, say.

If files are very large, you can parallelize work into "map-reduce" form by specifying the chromosome name as an option:

$ closest-features --chrom chr1 --closest --delim '\t' intervals.bed genes.bed > answer.chr1.bed
$ closest-features --chrom chr2 --closest --delim '\t' intervals.bed genes.bed > answer.chr2.bed
...

To reduce, just bedops -u all the answer.*.bed files.

Running these per-chromosome jobs in parallel makes the work take about as long as searching chr1, which tends to be the largest chromosome.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.