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?
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.
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.
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
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.
I think HTSeq does almost that. It outputs a matrix of read counts per gene given a fastq sample and annotation file