SEQanswers

Go Back   SEQanswers > Applications Forums > RNA Sequencing



Similar Threads
Thread Thread Starter Forum Replies Last Post
assigning function to genomic coordinates Dessi Bioinformatics 4 06-29-2014 07:20 AM
RNA-Seq: VESPA: software to facilitate genomic annotation of prokaryotic organisms th Newsbot! Literature Watch 0 04-07-2012 03:33 AM
Exon-Junction mapping: re-assigning CDS-mapped reads to chromosomes sridharacharya RNA Sequencing 1 10-21-2010 04:07 PM
how to get annotation form genomic region? dwb0211 Bioinformatics 7 01-13-2010 09:35 PM

Reply
 
Thread Tools
Old 03-02-2016, 09:45 PM   #1
moldach
Member
 
Location: Vienna

Join Date: Jan 2016
Posts: 10
Default Assigning reads to genes in the absence of genomic annotation

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)?
moldach is offline   Reply With Quote
Old 03-03-2016, 12:51 AM   #2
SylvainL
Senior Member
 
Location: Geneva

Join Date: Feb 2012
Posts: 175
Default

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...
SylvainL is offline   Reply With Quote
Old 03-03-2016, 01:32 AM   #3
colindaven
Senior Member
 
Location: Germany

Join Date: Oct 2008
Posts: 399
Default

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.
colindaven is offline   Reply With Quote
Old 03-03-2016, 11:06 AM   #4
moldach
Member
 
Location: Vienna

Join Date: Jan 2016
Posts: 10
Default

Quote:
Originally Posted by SylvainL View Post
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 is offline   Reply With Quote
Old 03-03-2016, 11:44 AM   #5
moldach
Member
 
Location: Vienna

Join Date: Jan 2016
Posts: 10
Default

Quote:
Originally Posted by colindaven View Post
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.

Quote:
Originally Posted by colindaven View Post
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 is offline   Reply With Quote
Old 03-07-2016, 11:16 AM   #6
moldach
Member
 
Location: Vienna

Join Date: Jan 2016
Posts: 10
Default

Can anyone help?
moldach is offline   Reply With Quote
Old 03-07-2016, 10:35 PM   #7
SylvainL
Senior Member
 
Location: Geneva

Join Date: Feb 2012
Posts: 175
Default

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
SylvainL is offline   Reply With Quote
Old 03-07-2016, 11:27 PM   #8
colindaven
Senior Member
 
Location: Germany

Join Date: Oct 2008
Posts: 399
Default

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/
colindaven is offline   Reply With Quote
Old 03-09-2016, 01:20 PM   #9
shi
Wei Shi
 
Location: Australia

Join Date: Feb 2010
Posts: 235
Default

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

Tags
annotation, count, deseq, rnaseq

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 06:17 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2018, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO