Assigning reads to genes in the absence of genomic annotation

[moldach]

I'm looking for a method to calculate read counts from SAM/BAM files that does not require an annotation file (GTF or BED files).

I've done a de novo transcriptome assembly of an organism which does not have an annotated genome. While there is an available genome for a closely related species, it has not been annotated.

Which program(s) can be used to output count data to a text file suitable for differential expression analysis with DESeq?

If I blastx my transcriptome against uniref 50 (or uniprot_sprot, I'm not sure which database is ideal - but that is another post all in itself) how can this output be integrated with the count data set (either in DESeq or prior to that)?

[SylvainL]

Hi,

did you map your reads versus the closely realted genome of directly on your transcriptome assembly? If you did on the transcriptome, simply use bedtools to get the number of reads on each transcripts...

[colindaven]

Suggestions

#1 Remap reads to closely related genome. Satisfied with mapping rate ?

#2 Use gmap (easy) or Maker to map your de novo assembled transcripts to the related genome. Again, satisfied ? View both sets in a genome browser.

#3 if unsatisfied with #2 perhaps use Trinity genome guided or cufflinks to recreate transcripts.

#4 Quantify - ie using featureCounts - transcripts from #2 or #3.

Forget blast for this kind of approach.

[moldach]

Hi,

did you map your reads versus the closely realted genome of directly on your transcriptome assembly? If you did on the transcriptome, simply use bedtools to get the number of reads on each transcripts...

I mapped directly on my transcription assembly.

I couldn't find any reference to getting the number of reads on each transcript (maybe it's just worded differently?) from the documentation of bedtools. However, I found a Biostars link that suggested using the multicov sub-command in the bedtools suite.

However, according to the documentation the multicov from BEDtools requires genome annotation. For example:

>bedtools multicov –bams run.bam -bed genes.bed

Are you talking about another sub-command or can multicov be run without the bed file?

[moldach]

Suggestions

#1 Remap reads to closely related genome. Satisfied with mapping rate ?

I had tried mapping at one point some-time-ago to the closely related un-annotated genome. Unfortunately, I used Bowtie2. I now know better; you need to use a splice-junction aware aligner.

Suggestions
#2 Use gmap (easy) or Maker to map your de novo assembled transcripts to the related genome. Again, satisfied ? View both sets in a genome browser.

So GMAP maps and aligns with this command:

>gmap -d <genome> -A <cdna_file>

And it would output SAM files.

What I don't understand is how (or if) GMAP annotates this genome?
The documentation for maker on the other hand clearly states it annotates but I can't find anything in the GMAP documentation.

Will gmap and Maker output an annotation file including chromosomal coordinates of features (GTF)? It says that this is a required file to use featureCounts

[moldach]

Can anyone help?

[SylvainL]

Hi,

since you aligned directly on your transcriptome, I guess your reference contains all the transcripts, so you can get the counts for each by using samtools idxstats

[colindaven]

A GMAP command which produces GFF3 output might look like this:

~/gmap-2015-07-23/bin/gmap -f gff3_gene -D gmap/ -d mygenome.fasta.gmap -B 5 -t 12 --intronlength=50000 --totallength=1000000 -p 3 --npaths=1 transcripts.fa > transcripts.gff3

This is a nice GFF3 which can be used directly by "bedtools multicov"

If you want to use featureCounts for read counting try using ngsutils to convert from gff3 to gtf.

http://ngsutils.org/modules/gtfutils/

[shi]

featureCounts works with both GTF and GFF formats. I think it should be fine if you directly provide your GFF3 annotation to featureCounts program.