Suppose I have a list of genomic positions and a reference GTF/ GFF3. Are there any software tools (any language is fine) that maps the position to the cDNA/ CDS position in relate to any transcript(s) that it may overlap?
For example, based on hg38 and GENCODE version 25, I would expect an input of chr1:151406413 would give the following as output (using CDS coordinates): ENST00000271715: 2622, ENST00000368863: 2337, ENST00000392723: 2463, ENST00000409503: 2595, ENST00000491586: 2490, ENST00000529669: 822 and ENST00000531094: 2436.
I know I can just throw it into ANNOVAR or VEP; but I wonder if there's a simpler way without actually predicting the consequences.
Thanks anyone for the answer in advance!
In R, the GenomicFeatures package could be one choice.
Eventually I got to try out this one. The important GenomicFeatures function here is mapToTranscripts:
library(GenomicFeatures)
library(rtracklayer)
gencodeTxDb <- makeTxDbFromGFF (file="gencode.annotation.gtf")
gencodeTx <- transcripts (gencodeTxDb)
names (gencodeTx) <- id2name (gencodeTxDb, "tx")
genomicLocus <- import.bed (con="genomicLocus.BED")
mappedLocus <- mapToTranscripts (x=genomicLocus, transcripts=gencodeTx)
The output is an GRanges object that is easy enough to work with.
On the command line, one option is to download annotations and convert them to BED via convert2bed:
$ wget -qO- ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_25/gencode.v25.annotation.gff3.gz \
| gunzip -c- \
| convert2bed -i gff - \
> annotations.bed
If you just want transcripts, you could filter the result with awk:
$ awk '$8=="transcript"' annotations.bed > transcripts.bed
Convert your input to a sorted BED file with awk and sort-bed:
$ awk FS=":" '{ print $1"\t"$2"\t"($2 + 1); }' input.txt | sort-bed - > input.bed
Then query the input.bed file against annotations.bed or transcripts.bed with bedmap:
$ bedmap --echo --echo-map-id --delim '\t' input.bed transcripts.bed > answer.bed
Once you have made annotations.bed or transcripts.bed, you can repeat queries with any other input.bed file you generate down the road.
Log in to answer this question.