Dear Sourav, I am a bit surprised by your question, since RNA-seq is a relatively well covered area in bioinformatics and several tools are quite collectively accepted depending on what output you want to get. If you're only interested in expression in terms of whole genes, I would consider using your bam file as an input for counting scripts like featureCounts (subread.sourceforge.net/) or HTseq-counts ( https://htseq.readthedocs.io/ ). That'll give you raw counts data (integers) that can be analyzed for differential expression between conditions through R (https://www.r-project.org/) packages like DESeq2 (https://bioconductor.org/packages/release/bioc/html/DESeq2.html).

If you're more interested in transcripts detection and RPKM / FPKM values, then you need to consider transcripts sizes (i.e. sum of exons sizes) but there are already some pipelines dealing with that, particularly to detect which isoforms are the more likely to be expressed, the more popular being Tuxedo (bowtie / tophat / cufflinks): http://cole-trapnell-lab.github.io/cufflinks/.

I advise you to also have a look on basic introductions to RNA-seq analysis here: bioinformatics.ucdavis.edu/docs/.../Th_MB_RNASeq_Intro.pdf or here: https://www.rna-seqblog.com/introduction-to-rna-sequencing-and-analysis/

Good luck with your analyses, hope that'll help!

Dear Sourav, I totally agree with Fabrice. Currently one of the best way to analyse RNA-seq data is to use STAR for the mapping then featureCounts for counting the mapped reads per gene. You have the option to count the reads per exon or per gene with featureCounts and it usually gives the same result plus or minus few reads. The easiest to have a clean output is to count per gene.

For TPM/FPKM, it seems that TPM has a better reputation than FPKM/RPKM these days (i.e. some reviewers can specifically ask for TPM).

One of the easy solution would be to use already made and validated RNA-seq pipelines such as the one available here (from fastq files to DE, MAplots and so on):

https://github.com/maxplanck-ie/snakepipes

The code is transparent so if you have a doubt about how this or that step is made, you can check it and even bring your own modifications to the pipeline.

As an alternative, to identify DEG, there's this very nice and comprehensive pipeline under R using Limma and EdgeR (starts from gene counts):

https://f1000research.com/articles/5-1408/v2

Dear Fabrice and Gautier,

Thanks a lot for your answers. Actually the organisms we are working on is a bit unexplored through RNA-Seq so many things we are performing manually. Probably that is why we are not using any pipeline but designing our own because we have a couple of check points. It is true that we have tried couple of tools for counting and found results differs. But I believe I need to explore literature more as you suggested. Thanks for the links again.

In RNA-seq provide the level to understanding the RPKM / FPKM /TPM is the way to normally of the reading bias during the mapping but need to understand the read not really represent the gene expression. It is depend on what tools are you use for stats the gene count. So carefully if you normalize not do double normalize via the DEseq2 or EDGR tools because in both tools already had the script to calculate the size factor and normalizes.

Gautier Richard Thank you! Your answer really helps!