I'm trying to compute a gene expression profile for an organism. I have gene nucleotide sequences of the mentioned organism stored in a fasta file and a set of paired reads stored in two separate files with the same fasta format. Now I want to compute the coverage of every gene by the given reads. That is, on average how many reads cover a given segment of the sequence.
Is there any tool for this specific purpose? I'm trying to use the velvet sequencer but I have encountered a few problems.
It's unclear if your paired-end reads are actually in fasta format, I'll presume that they're in fastq instead.
The easiest tool to use is salmon, which nicely deals with things like multimapping. If you're trying to judge the quality of an assembly, then I recommend having a look at transrate, which uses some related methodology for assessing contig quality.
It's not possible to compute absolute expression from RNASeq reads if they are processed in the usual way, where a sequencer produces the same number of reads regardless of the input RNA amount. At best, RNASeq will give you an indication of proportional expression within a single sample. For this reason, relative expression (i.e. that used by differential expression tests) is easier to determine than absolute expression.
The closest approximation to absolute expression is to generate an expression relative to the average expression of a set of housekeeping genes, but I don't think there's a universal set that has been decided on. Gene expression, even for common housekeeping genes, can vary depending on the environmental conditions of the cell. For example, GAPDK is involved in immune cell activation.
However, as long as experimental conditions are similar, and you're not planning on looking for statistical significance, the proportional expression can still give qualitative insights into how cell populations behave in relation to other populations. DESeq2 provides a variance-stabilising transformation function that minimises variation for small-count genes, assuming that each sample has roughly the same total expression. I have found that I get better outcomes / comparisons from this transformation when carrying out a further adjustment to account for gene length (i.e. divide by the length of the longest transcript for each gene). See our Th2 paper, section "Read mapping and differential expression analysis" for more information. The "transcripts-per-million" values produced by Kallisto and Salmon provide measures similar to this.
If, on the other hand, you were able to modify the experimental design, single cell sequencing (or "known cell count" sequencing) can be used for determining absolute expression: use a spiked-in transcript that is added in proportion to the cell count, so that results can be compared in proportion to the expression for that transcript.
