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
•
written
by
Kevin_Smith •
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.
I'm using pybedtools, do you know how I can use your code with this module?
Thanks!
To use with Python, call with subprocess or similar library.
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?
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.
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
Log in to answer this question.
