6
$\begingroup$

I hope this is a good place to ask such question. I have to do some data analysis on RNA-seq data from human cells. I am currently searching for tools to help me with that. Specifically, I would need some tools to analyze the gene expression from the data. Something to help me plot the expression of selected genes in each fastq file and compare the differences in the expression with the possibility to export the results or some command line interface for scripting. Basically I need something where I can put a fastq file and perhaps also a human genome annotation file as input and get gene expression as output. I have looked at bioconductor and it's packages and on Wikipedia's List of RNA-Seq bioinformatics tools. I suppose some of these tools have to be able to do what I need, but I have been unable to find out which one and how should they be used to achieve that. Could someone please give me some advice?

$\endgroup$
6
$\begingroup$

You will likely need a tool to "map" the reads on the reference genome. You may find such a reference genome, together with annotations, here: ftp://ussd-ftp.illumina.com/.

Mapping tools such as bowtie2 or bwa take fastq files and reference genomes and output mapping results in a format called sam.

You then have a lot of options to estimate gene expression.

  • You can write your own algorithm to parse sam format and estimate normalized read counts on each gene.

  • You can combine more or less low-level tools such as samtools, pysam, htseq with some scripting to do this.

  • You can use tools that do the counting (like bedtools ot htseq-count) and differential expression analysis (like deseq2).

In the last case, I would advice to start from the documentation of the final tool to find out what are the tools you need to generate the output of the preceding step.

It is very likely you will use some R or Python, or use the web platform galaxy for some of the steps.

Edits

As mentioned by @scribaniwannabe in this answer, the paper about the Tuxedo suite of tools gives a good example of the steps to carry out an RNA-seq analysis using recent tools (as of October 2016).

As @Student T reminds in this answer, RNA-seq data contain reads that can come from exon-exon junctions, so the read mapper has to be set up in such a way as not to discard reads not mapping continuously on all their length on the genome. To my knowledge, HISAT2 and CRAC do this by default. Bowtie2 needs special settings.

$\endgroup$
  • 1
    $\begingroup$ TopHat is the part of the tuxedo suite that does splice-aware mapping. $\endgroup$ – swbarnes2 Oct 21 '16 at 16:34
  • $\begingroup$ @swbarnes2 You are right concerning the older tuxedo suite (the new one uses HISAT2). Curiously, TopHat uses bowtie2 internally, and it seems that it doesn't use the "--local" mode (look for "local" in the tophat executable, which happens to be a complex python program). I don't know how it achieves splice-awareness. $\endgroup$ – bli Oct 21 '16 at 16:43
3
$\begingroup$

While I also agree @bli that R and Python (in particular Bioconductor) have more than enough packages for you to compare gene expression. You shouldn't align your reads with bwa or bowtie because they don't take introns into consideration. You should use TopHat or STAR.

$\endgroup$
2
$\begingroup$

The answer @bli gave is great. I thought I would point out that Johns Hopkins also recently upgraded their tuxedo suite. Looks promising and has great instructions for use.

Also, I've begun to grow quite fond of the GeneTrail 2 tool for my RNA-Seq secondary analysis. Gives great results for enrichment analyses.

Hope this is helpful.

$\endgroup$
  • $\begingroup$ Lastly, I also like to use [HOMER][3] for de novo motif finding up/downstream of the promoter regions of deferentially expressed genes. [3]: homer.salk.edu/homer/motif $\endgroup$ – aaiezza Oct 10 '16 at 15:25
  • 1
    $\begingroup$ Thanks for the link about the tuxedo suite. The explanations in the article seem useful. $\endgroup$ – bli Oct 21 '16 at 13:12
2
$\begingroup$

I think that STAR is the preferred splice-aware aligner nowadays. STAR can output counts by gene or by transcript. Assuming you have Illumina data, you can try using the tools on Illumina's BaseSpace. RNASeq might be one of the things that you can do for free there.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.