This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Annotating gene ID and genomic features from bed files

Hello, I have a bed file from a CLIP-seq experiment that has the columns: chr, start, stop, -log10(p-value), log2fold enrichment and strand (The bed file was generated with bedtools intersect of three bed files). I wish to annotate the intervals with gene name and genomic features (3'UTR, exon, intron, etc). I have found several ways to just annotate the gene name but I would really like to include the genomic features if possible. Any recommendations? Thanks!

genome clip
• 1,993 views • link
updated by Alex Reynolds Alex Reynolds • written by nornorbinks • nornorbinks

1 answer

Using BEDOPS convert2bed, you can make a file of named features:

$ wget -qO- ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/gencode.v28.annotation.gff3.gz \
    | gunzip --stdout - \
    | convert2bed -i gff - \
    | awk -vFS="\t" -vOFS="\t" '{ print $1, $2, $3, $4"|"$8 }' - \
    > named-features.bed

Adjust the GENCODE download URL as needed for your reference genome and desired release version.

Then map your sorted BED file to the named features (assuming your BED file is called clip-seq.bed):

$ bedmap --echo --echo-map-id --delim '\t' <(sort-bed clip-seq.bed) named-features.bed > answer.bed

The name and feature category will be in the last column of answer.bed.

If you need to translate Ensembl IDs to HGNC names, one option is to modify the awk statement above to parse out the HGNC name from the optional attributes in the GFF file. As they may not always be present, another option is the use of the mygene service. An example Python script is available here that demonstrates retrieval for mm10 records: A: Get gene names from ensembl ID or gene region

• 0 views • link
written by Alex Reynolds Alex Reynolds

Log in to answer this question.

More posts like this