Newbie question regarding mapping of RNA-seq data

[analog900]

Hi all,
I'm stuck with a total newbie problem here. I'm analyzing RNA-seq data from mouse, I mapped the (paired-end) sequences using TopHat against mm9 (using bowtie1) but when I look at the SAM output files, the hits list chromosomes as map targets, when instead I'm interested in gene IDs. I'm assuming I missed something trivial?

[GenoMax]

Since you did the mapping against genome you need to summarize the alignments using a program like featureCounts or HTSeq-count along with an annotation file that will translate the alignments you have into counts per gene/exon (any features included in the annotation file).

You could have also provided that annotation file to TopHat (when you ran it) if you only wanted to look at the transcriptome (instead of the whole genome).

[dpryan]

That's what should happen. Your next step is to get counts of aligned fragments per gene, for which you can use featureCounts or htseq-count. Both of those expect exactly what you have as input.

Edit: Genomax beat me by a minute. I should note that mapping against the transcriptome with tophat still produces alignments in genomic coordinates.

[analog900]

Thanks guys!! Appreciate it! I'll give it a try.
Thanks again!

[analog900]

Using featureCounts gives me nice summary and counts text files. However, my SAM and BAM files still contain the original, genomic, annotations (obviously). Ideally, I would like to convert the annotations in the BAM/SAM files so that I can further process them.

This leads me to a more broader question: what reference (for mouse rna-seq) do people use when they want gene_ids instead of genomic targets?. I noticed that reference files such as mRNA.fa or refMrna.fa only contain accession numbers, but not gene ids.
Thanks in advance

[Brian Bushnell]

Gene IDs, names, and numbers vary depending on the database in question. You can either get a translation table, or try find a fasta file already named with the identifiers you want to use.

[GenoMax]

Using featureCounts gives me nice summary and counts text files. However, my SAM and BAM files still contain the original, genomic, annotations (obviously). Ideally, I would like to convert the annotations in the BAM/SAM files so that I can further process them.

What is "further processing" referring to here? Most downstream analysis is going to use the counts files (unless you are going to call SNPs from this data) and will always refer to the gene names contained in that file.

[analog900]

What is "further processing" referring to here? Most downstream analysis is going to use the counts files (unless you are going to call SNPs from this data) and will always refer to the gene names contained in that file.

I've been loosely following the "simple fool's guide for rna seq" by the group of Stephen Palumbi (http://sfg.stanford.edu/guide.html). They parse their SAM output files with a series of python scripts to obtain similar summary statistics like the ones I can now get with featureCounts. Then, they use DESeq for functional enrichment (which I would really like to do in order to compare my different samples).

[dpryan]

I would recommend ignoring that guide. If you want to use DESeq (use DESeq2), just directly use the counts from featureCounts. This would be the standard and accepted pipeline and there's no reason to use any kludgy scripts.

[analog900]

I would recommend ignoring that guide. If you want to use DESeq (use DESeq2), just directly use the counts from featureCounts. This would be the standard and accepted pipeline and there's no reason to use any kludgy scripts.

Thank you. Really appreciate it! Can you recommend any other standard/accepted pipelines downstream of featureCounts?

[shi]

We use limma/voom and edgeR in downstream analyses to discover differentially expressed genes. The link below is a short tutorial for using our pipeline for analyzing RNA-seq data which you might find helpful:

http://bioinf.wehi.edu.au/RNAseqCaseStudy/

[Michael Love]

for DESeq2 you would use the DESeqDataSetFromMatrix function to start the analysis, using the counts matrix returned by featureCounts. Example of starting from count matrix is in the DESeq2 vignette.

[analog900]

Thanks so much guys!
Working through the DESeq2 vignette now and learning new stuff... really excited!
Thanks again!