function/abundance prediction from reads

[thedude]

I've got a more general question. Many of the recent metagenomics studies deal with read lengths that complicate assembly or render it impossible. After some preprocessing people hence blastx (or translate and blastp, HMMer etc.) these reads against nucleotide or protein databases (or COGs, KEGG etc.) to annotate them. So far so good, but none of these papers ever answer my simple question: as I understand it 454, Solexa, actually any kind of sequencing method, will produce many reads that map to the same gene of the same genome, just 'shifted' slightly. That means for one say dna polymerase gene i'd get 10 reads with similarly good blast scores or hmm hits. What I don't get is how can these papers say 'we got 10 hits in this and that COG category' - how can they be sure it's not reads of the same gene, the same physical piece of the DNA sample (I'm not talking about very similar genes from closely related species here, I'm talking about the exact same piece of DNA, so it's more a technical question for that bit). If my idea of the multiple reads per gene is true and you functionally annotate from the reads directly (as outlined above) you can only talk about relative abundances (of hits) but not about actual counts, no? Help would be much appreciated.

Cheers.
Rob

[risc303]

There are many reads matching one transcript. For each transcript (mRNA isoform) you'd get 100s reads. Count them all for each transcript, then for each gene then for each GO term. "Digital gene expression".

[thedude]

Mmm...not sure if I get that. I actually meant the raw DNA reads as produced in metagenomic studies. Not metatranscriptomics or so, so no mRNA isoforms. What I mean is that the same part of the same gene (from exactly the same genome) will be represented by a varying number of reads - even if that gene would appear only once in the whole sample. Then BLASTing these reads against families (like e.g. COG) would make this family appear to get several hits (several genes with that function). So you can never say 'I found 10 DNA polymerase genes' cause the hits could all come from the same single sequence having been 'read' several times, right? You'd have to bin all reads that match with a given score over a given length first and assume they all represent the same gene from a single genome. But that's not what people do in the studies I've read. That's why I'm still slightly confused...

Thanks for your reply!

[risc303]

Now I see what you mean. I agree, I'd also assemble a gene from multiple reads first and then count "DNA Pol" genes... Which papers are you referring to?
It could be that in that metagenomics studies they don't get complete coverage of a genome with there reads, so there is space between (one ore less reads per gene), but the reads itself are long enough to identify a family.