Annotate each peak to gene starts within 10 kb of the peak with python ?

Question

[Kevin_Smith]

How can I do this with peaks in a bed file and genes in a txt file ? Also I will like to see the nearest gene to each peak.

Thank you very much

Answer 1

[Alex Reynolds]

Via command-line toolkit, to map a set of names of overlapping genes to each peak within 10kb:

$ bedmap --echo --echo-map-id-uniq --range 10000 peaks.bed genes.bed > annotation-answer.bed

To get the nearest gene to each peak:

$ closest-features --closest peaks.bed genes.bed > nearestgene-answer.bed

To use with Python, call with subprocess or similar library.

[Kevin_Smith]

I'm using pybedtools, do you know how I can use your code with this module?

Thanks!

[Alex Reynolds]

To use with Python, call with subprocess or similar library.

[Kevin_Smith]

I'm new with python , I'm using pybedtool. I don't know how to do it. Do you can give me more specific instructions?

[Alex Reynolds]

For example:

import subprocess
subprocess.call('bedmap --echo --echo-map-id-uniq --range 10000 peaks.bed genes.bed > annotation-answer.bed', shell=True)

More details are available via the documentation. Search engines will also point you to a good review of the subprocess library, among others.

[Fidel]

Using pybedtools you can do something like the following (I simply adapted the code found in pybedtools documentation):

from pybedtools import BedTool

peaks = BedTool('peaks.bed.gz')
genes = BedTool('hg19.gff')

nearby = genes.closest(peaks, d=True, stream=True)

for gene in nearby:
    if int(gene[-1]) < 10000:
        print gene.name