Transferring genomic features on new coordinates

Question

I have a eukaryotic genome for which an updated sequence for a chromosome was recently obtained. I want to map RNAseq reads on the genome (and perform other downstream analyses) and would like to use the most up-to-date information possible (so the 'new' sequence of that chromosome).

However, before performing the read mapping, I'd like to update the annotation of my genome (GFF3) to match with the 'new' coordinates of this chromosome.

Basically, I would like to transfer information from the 'old-coordinates' GFF3 to obtain a 'new-coordinates' GFF3, and if possible conserve all the information/hierarchy of the file (gene, mRNA, exon, etc.).

Example:

chrXX   source  gene    222 5942    .   -   .   ID=gene_1;Name=gene_1;length=5720
chrXX   source  mRNA    222 5942    .   -   .   ID=gene_1.1;Parent=gene_1;Name=gene_1.1;length=5720
chrXX   source  exon    222 5794    .   -   .   ID=gene_1.1.2;Parent=gene_1.1
chrXX   source  exon    5889    5942    .   -   .   ID=gene_1.1.1;Parent=gene_1.1
chrXX   source  CDS 222 5794    .   -   1   ID=CDS:gene_1.1.2;Parent=gene_1.1;Name=gene_1.1
chrXX   source  CDS 5889    5942    .   -   0   ID=CDS:gene_1.1.1;Parent=gene_1.1;Name=gene_1.1

... should be updated to ...

chrXX   source  gene    333 6053    .   -   .   ID=gene_1;Name=gene_1;length=5720
chrXX   source  mRNA    333 6053    .   -   .   ID=gene_1.1;Parent=gene_1;Name=gene_1.1;length=5720
chrXX   source  exon    333 5905    .   -   .   ID=gene_1.1.2;Parent=gene_1.1
chrXX   source  exon    6000    6053    .   -   .   ID=gene_1.1.1;Parent=gene_1.1
chrXX   source  CDS 333 5905    .   -   1   ID=CDS:gene_1.1.2;Parent=gene_1.1;Name=gene_1.1
chrXX   source  CDS 6000    6053    .   -   0   ID=CDS:gene_1.1.1;Parent=gene_1.1;Name=gene_1.1

The approach I tried (but I am not sure it is the way to do it):

1. Extract sequence of all the features that are going to be updated to fasta (bedtools getfasta)
2. Map these sequences to the genome with the updated chromosome (gmap with --nosplicing since sequences we are mapping correspond to genomic regions).
3. Create the updated GFF3 file. When doing that, I apply one supplementary rule: if a feature mapped on another chromosome, but there was also another alignment on the updated chromosome, prioritize the one that correspond to the same (updated) chromosome.

What would be the proper way to do such a thing? Any suggestion of methods/tools is welcome!

Answer 1

I think the standard way of doing this is to make a chain file, then use it to do a liftOver of the annotations:

mkdir psl
for i in ../ci3/rm/masked/*.masked; do blat ../ci2.2bit $i -tileSize=12 -fastMap -minIdentity=98 psl/`basename $i .fa.masked`.psl -noHead -minScore=100; done

Translate psl files to chains in the directory chain:

mkdir chain
for i in psl/*.psl; do axtChain -linearGap=medium -psl $i ../ci2.2bit ../ci3/ci3.2bit chain/`basename $i .psl`.chain; done

Merge short chains into longer ones into the directory chainMerge:

mkdir chainMerge
chainMergeSort chain/*.chain | chainSplit chainMerge stdin -lump=50

concat and sort the chains:

cat chainMerge/*.chain > all.chain
chainSort all.chain all.sorted.chain

Need info about chromosome sizes for netting:

twoBitInfo ../ci3/ci3.2bit ci3.chromInfo
twoBitInfo ../ci2.2bit ci2.chromInfo

Netting: identify alignable regions from chains:

mkdir net
chainNet all.sorted.chain ci2.chromInfo ci3.chromInfo net/all.net /dev/null

Finally, select the right alignable regions using the nets, creating a "liftOver" file:

netChainSubset net/all.net all.chain ci2ToCi3.liftOver

Run liftOver:

CrossMapy.py bed ci2ToCi3.liftOver test.hg18.bed

Answer 2

In addition to the answer by heathobrien you might take a look at RATT: Rapid Annotation Transfer Tool. The approach is essentially the same, only using more modern approaches to for instance whole genome alignments. It uses nucmer, which when using the latest version (4.0) allows you to do whole genome alignments using parallel computing.