Transcript Coordinate Ranges to Genomic Coordinates

Question

I have 2 GFF3 files:

1. Features using transcript IDs as the landmarks. i.e. "CDS" feature types using coordinates from transcript space.
2. 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.

Answer 1

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='')

Answer 2

The linked gist contains a class I wrote for doing this on GTF files:

https://gist.github.com/IanSudery/d8349c22823a475ceb489c3e8aeb448e

It uses the GTF class from cgat, which can be found here: https://github.com/cgat-developers/cgat-apps

You would use it to do this task like so:

from cgat import gtf
from cgat import iotools
import TranscriptCoordsInterconverter


trans_coord_file = gtf.iterator(iotools.open_file(args[1])
transcripts = dict()

for transcript in gtf.transcript_iterator(trans_coords_file):
    transcript_id = transcript[0].transcript_id
    transcripts[transcript_id] = transcript

trans_coord_file.close()

genome_coord_file = gtf.iterator(iotools.open_file(args[2]))
for transcript in gtf.transcript_iterator(genome_coord_file):

    transcript_id = transcript[0].transcript_id
    contig = transcript[0].contig
    strand = transcript[0].strand

    if transcript_id not in transcripts:
        continue

    converter = TranscriptCoordsInterconverter(transcript)
    for exon in transcripts[transcript_id]:
        genome_intervals = converter.transcript_interval2genome_intervals((exon.start,exon.end))
        for interval in genome_intervals:
            new_entry = GTF.Entry().fromGTF(exon)
            new_entry.contig = contig
            new_entry.strand = strand
            new_entry.start = interval[0]
            new_entry.end = interval[1]
            print (str(new_entry)+"\n")