Transcript Coordinate Ranges to Genomic Coordinates I have 2 GFF3 files:


Features using transcript IDs as the landmarks. i.e. "CDS" feature types using coordinates from transcript space.
Features using chromosome IDs as the landmarks. i.e. "exon" feature types using coordinates from chromosome space.


I would like to transform the coordinates of the features in file 1 into the coordinate space of file 2. i.e. transcript-based coordinates into genomic-based coodinates.

Here is an example for file1:

cat transcript_orfs.gff3
##gff-version 3
##sequence-region Tx.1 1 4000
Tx.1    ORF_finder  gene    1   4000    .   +   .   ID=1
Tx.1    ORF_finder  mRNA    1   4000    .   +   .   ID=2;Parent=1
Tx.1    ORF_finder  exon    1501    2500    .   +   0   ID=3;Parent=2
Tx.1    ORF_finder  CDS 1501    2500    .   +   0   ID=4;Parent=2

gt sketch -addintrons transcript_orfs.png transcript_orfs.gff3




Here is an example for file1:

cat genome.gff3
##gff-version 3
##sequence-region chr3A_part1 1 454103970
chr3A_part1 genome_assembly gene    1001    6000    .   +   .   ID=Tx
chr3A_part1 genome_assembly mRNA    1001    6000    .   +   .   ID=Tx.1;Parent=Tx
chr3A_part1 genome_assembly exon    1001    3000    .   +   .   ID=Tx.1.exon1;Parent=Tx.1
chr3A_part1 genome_assembly five_prime_UTR  1001    2000    .   +   .   ID=Tx.1.utr5;Parent=Tx.1
chr3A_part1 genome_assembly CDS 2001    3000    .   +   0   ID=Tx.1.cds1;Parent=Tx.1
chr3A_part1 genome_assembly exon    4001    6000    .   +   .   ID=Tx.1.exon2;Parent=Tx.1
chr3A_part1 genome_assembly CDS 4001    5000    .   +   2   ID=Tx.1.cds2;Parent=Tx.1
chr3A_part1 genome_assembly three_prime_UTR 5001    6000    .   +   .   ID=Tx.1.utr3;Parent=Tx.1

gt sketch -addintrons genome.png genome.gff3




I want to convert the coordinates of the features from file1 into chromosome genomic coordinates. I expect to get something like the following output:

cat output.gff3
##gff-version 3
##sequence-region chr3A_part1 1 454103970
chr3A_part1 ORF_finder  gene    1001    6000    .   +   .   ID=1
chr3A_part1 ORF_finder  mRNA    1001    6000    .   +   .   ID=2;Parent=1
chr3A_part1 ORF_finder  exon    2501    3000    .   +   0   ID=3.1;Parent=2
chr3A_part1 ORF_finder  CDS 2501    3000    .   +   0   ID=4.1;Parent=2
chr3A_part1 ORF_finder  exon    4001    4500    .   +   0   ID=3.2;Parent=2
chr3A_part1 ORF_finder  CDS 4001    4500    .   +   1   ID=4.2;Parent=2

gt sketch -addintrons output.png output.gff3




I have looked into using mapFromTranscripts() from Bioconductor's GenomicRanges library but I have made little progress trying to decipher the manual.
This task has the same "flavor" as many I've done before, but each case is so subtly different that it's impossible to write a generalized tool that will work correctly for all circumstances.

R isn't my wheelhouse, but I was able to throw something together pretty quickly with standard Python. This should get you at least 95% of the way to where you're going. :-)

#!/usr/bin/env python3

from argparse import ArgumentParser, FileType
from collections import defaultdict
from re import search, sub

cli = ArgumentParser()
cli.add_argument('txspace', type=FileType('r'))
cli.add_argument('chrspace', type=FileType('r'))
args = cli.parse_args()

exons = defaultdict(list)
for line in args.txspace:
    if line.startswith('#') or line.strip() == '':
        continue
    values = line.split('\t')
    featuretype = values[2]
    if featuretype not in ('CDS', 'exon'):
        continue
    txid = values[0]
    exons[txid].append(values)

for line in args.chrspace:
    values = line.split('\t')
    if len(values) != 9:
        print(line, end='')
        continue
    featuretype = values[2]
    if featuretype not in ('gene', 'mRNA'):
        continue
    print(line, end='')
    if featuretype == 'mRNA':
        txid = search('ID=([^;\n]+)', values[8]).group(1)
        start = int(values[3])
        txexons = exons[txid]
        for exon in txexons:
            exon[0] = values[0]
            exon[3] = str(int(exon[3]) + start)
            exon[4] = str(int(exon[4]) + start)
            exon[8] = sub(r'Parent=[^;\n]+', 'Parent=' + txid, exon[8])
            print('\t'.join(exon), end='')