Generating multi-species alignments from GTF

[dagarfield]

Hiya folks,

I've got a two part question for you today. I'm attempting to create multi-species alignments for use in a scan for selection. I've one well-sequenced reference and several related species that have been sequenced using 454.

One technique would be to take each gene in the reference genome as a target and use Newbler (or something else) to map the 454 reads to the genes. But I'd rather something more elegant, like an over.chain file (http://genome.ucsc.edu/goldenPath/help/chain.html) that allows me to pull out specific orthologous regions from a whole-genome alignment.

Here are my issues:
1) I can use Newbler to align the 454 reads to the reference genome, but the resulting file formats don't help very much. Does anyone know an elegant way to transform ace (or other Newbler files) into SAM or GTF or PSL? Or is the most "elegant" way through AMOS as suggested in older threads?

2) For a couple of species, I have .gtf files that have already been generated for me. Yay. Does anyone know of an existing program that would allow me to use these .gtf files to get pull down all 454 reads (or at least names) that map to a given region of the reference? Or a way to reverse to .gtf so that the map is REFERENCE-->454 Reads rather than the current 454 Reads --> REFERENCE?

Many thanks, and happy sequencing.

--David

[westerman]

Hi dagerfield. As far as I know there is no more "elegant" way to get a SAM file from 454 output. The Roche reps promise SAM formatted files "in the next release".

As for question #2, it depends on what is in your GTF file. The GTFs that I work with are simply annotation files for a genome. This by itself will not have information about reads. However you may have a different type of GTF file.

When you use gsMapper or runMapping, in the 454ReadStatus.txt file there will be information on where each read mapped to -- chromosome and position. It should be "easy" to convert this information into reference-->reads in some format.

[dagarfield]

Ah, yes. I was a bit unclear about the GTF files. The "features" described by this GTF are 454 reads that have been mapped to the reference genome. The goal is have a pipeline that allows me to specify a given region of the reference genome and pull down all reads from the other species that map to this region. I imagine I can write some sort of hash-table, but I was really hoping there'd be a mechanism to allow me to transform a gtf file into a .chain file such as those used on the UCSC genome browser.

Thanks for the tip about where to find relevant information in the 454 results directory.

-_DG

[jmw86069]

I'm not sure what tools you're ideally wanting to use, but it seems that GTF may not be the most efficient format for the purpose you described as I understand it. (GTF would be way down my list of desired formats for individually mapped reads.) Instead, a format like BAM seems efficient in terms of storing alignments and sequence reads, and allows you to pull all reads within a specific coordinate range fairly quickly. If you had a BAM file per species, you could quickly pull out individual reads and consensus for each species for a given set of gene coordinates on a reference genome. Samtools rocks the house for that.

Also, maybe it's just me, but I think of Newbler as more a assembly tool (and very good de novo assembler) than a mapper. Although you can use their mapping pipeline, you're probably better off either (1) converting their mapping output to BAM or a close relative, or (2) using Newbler to pre-assemble the reads, then map the isotigs/contigs and convert to BAM format for the assembled reads, or (3) using another mapping tool. The advantage of #2 is that it can cut down data volume substantially, though it depends upon how many reads you have per species whether it would make much difference and whether the contigs would be better or worse than individual reads. Also if it's Illumina data (and probably isn't since you mentioned Newbler) you'd probably want to pre-assemble. Longer contigs would be better if you had to compare more divergent species at the peptide level.

Depending upon how similar the species are, something simple like BLAT may suffice? samTools has a psl2sam.pl script, and I have also written a psl2bam.pl script which merges in the sequence reads.

Anyway I could've just said "try BAM format and samTools" but I'm long-winded like that. :-)

[dagarfield]

Hi,

It appears GTF isn't the best format for what I want to do. Unfortunately, a GTF file is what SpBase.org (the urchin model system database) provides. I was trying to avoid remapping. That said, your suggestion about a BAM file is a good one--it sounds like just the right trick. I'll give shifting the GFF file to SAM (and then to BAM) a try. If that doesn't work, then I'll probably just try BLAT.

To throw a tip back your way, if the species are far apart, I've had a lot of luck with lastz from Webb Miller's group.

Thanks,
DG