Mapping transcript coordinates to genomic positions

[kylepalos]

Hi there,

I hope this is the appropriate forum this sort of post. There seems to be general interest in mapping the predicted m6A sites back to the genome. The script from xpore that gets referenced in various issues may not be the most approachable for many users. A potentially easier solution for users more comfortable in R:

# generate minimal m6Anet output
m6Anet_output <- data.frame(transcript_id = "AT1G01040.1",
                            transcript_position = "100",
                            n_reads = "100",
                            probability_modified = "1",
                            kmer = "AGACA",
                            mod_ratio = "1")

# convert to more general bed format however the user prefers
# this will likely need incorporating other data types (e.g., transcript 2 gene mapping file)
m6Anet_bed <- data.frame(chr = "1", # this is the chromosome the transcript belongs to
                         start = 100, # this is the position of the predicted mod
                         end = 100, # same position
                         gene = "AT1G01040", # this is the geneID the transcript belongs to)
                         strand = "+") # strand of the gene

# convert to GRanges                         
library(GenomicFeatures)     
library(Repitools)
gr <- GRanges(seqnames = m6Anet_bed$gene, ranges = IRanges(start=m6Anet_bed$start,
                                                                   end=m6Anet_bed$end), 
              strand = m6Anet_bed$strand)

# load in annotation file
txdb <- makeTxDbFromGFF("annotation.gff3")

# isolate transcripts and genes from txdb
transcripts <- transcriptsBy(txdb, by = c("gene"))

# use mapFromTranscripts in the GenomicFeatures package
map2genome <- mapFromTranscripts(gr, transcripts)

# convert to dataframe
map2genomedf <- annoGR2DF(map2genome)

I hope this is useful and appropriate here. Please feel free to suggest users try this code if they have issues converting transcript to genomic coordinates.

[chrishendra93]

hi @kylepalos , thanks a lot for providing this. I will refer to your comment here in case someone needs help with the transcript mapping!

Best regards,

Christopher Hendra

[vetmohit89]

Hello,

I am wondering if there is anyway to functionally annotate ENST_id and position data to their cDNA functional annotation?

I have data in such format:

ENST00000337304.2 | 1367
ENST00000404241.6 | 616
ENST00000396680.3 | 36
ENST00000404241.6 | 667
ENST00000404241.6 | 617
ENST00000396680.3 | 52
ENST00000674920.3 | 260
ENST00000404241.6 | 618
ENST00000396680.3 | 889
ENST00000337304.2 | 1417

Thank you
Mohit

[hengjwj]

Hi @kylepalos,

Thank you for providing this piece of code for the community. Did you try reproducing the authors' findings or some other way to verify that the code works?

I ran m6Anet inference on SGNex_Hek293T_directRNA_replicate1_run1 using the provided HCT116 model and tried to convert the transcriptomic coordinates in data.site_proba.csv to genomic using the code here but found that I couldn't get all the same positions as in supplementary table 4. I got 118,438 rows in data.site_proba.csv but only managed to find 168 rows that matched the genomic coordinates in supplementary table 4. This number goes up to 1,728 if I offset the transcript_position by +1 but it's still much less than the 72,709 in supplementary table 4. Did you have better luck? Hope to hear from you!

Joel

Hi there,

I hope this is the appropriate forum this sort of post. There seems to be general interest in mapping the predicted m6A sites back to the genome. The script from xpore that gets referenced in various issues may not be the most approachable for many users. A potentially easier solution for users more comfortable in R:

# generate minimal m6Anet output
m6Anet_output <- data.frame(transcript_id = "AT1G01040.1",
                            transcript_position = "100",
                            n_reads = "100",
                            probability_modified = "1",
                            kmer = "AGACA",
                            mod_ratio = "1")

# convert to more general bed format however the user prefers
# this will likely need incorporating other data types (e.g., transcript 2 gene mapping file)
m6Anet_bed <- data.frame(chr = "1", # this is the chromosome the transcript belongs to
                         start = 100, # this is the position of the predicted mod
                         end = 100, # same position
                         gene = "AT1G01040", # this is the geneID the transcript belongs to)
                         strand = "+") # strand of the gene

# convert to GRanges                         
library(GenomicFeatures)     
library(Repitools)
gr <- GRanges(seqnames = m6Anet_bed$gene, ranges = IRanges(start=m6Anet_bed$start,
                                                                   end=m6Anet_bed$end), 
              strand = m6Anet_bed$strand)

# load in annotation file
txdb <- makeTxDbFromGFF("annotation.gff3")

# isolate transcripts and genes from txdb
transcripts <- transcriptsBy(txdb, by = c("gene"))

# use mapFromTranscripts in the GenomicFeatures package
map2genome <- mapFromTranscripts(gr, transcripts)

# convert to dataframe
map2genomedf <- annoGR2DF(map2genome)

I hope this is useful and appropriate here. Please feel free to suggest users try this code if they have issues converting transcript to genomic coordinates.