Hi I ran three genecount software tools (ht-seq, RSEM, Kallisto) to calculate genecount of RNA-seq data. For Ht-seq, i used STAR aligned Transcriptomesortedcordinate.bam file and defautl MAPQ score with intersection_nonempty mode. For RSEM, i used STAR aligner (used .gtf for building reference) to calculate expression, and Tximport to convert TPM into genecount. Kallisto does pesudoalignment (created index using transcriptome reference) and quantified the abundance using it and converted TPM into gene count using a transcript2gene map (created it from annotation file) file.
Ideally, thes counts from all three tools be approximately same in sample for a particular gene, i think. But they are extremely different. Please see the image.
I am going to visulaize the aligned sequences of that particular gene in a given region to find the reason if i can. In the meantime, any feedback, concerns or suggestions are most welcome!!
I don't think HTSEQ-count is nearly as smart at dealing with ambiguous reads as the other two. In fact, I think it mostly tosses them. So it's no surprise that it isn't registering many counts. As for the other two, they are different algorithms, and this is a hard case, splitting the reads between 4 rather similar genes, so it's not a surprise that there isn't a solid consensus here. Both are widely accepted, and I imagine if you have replicates, that each would treat the replicates more or less the same, and in a lot of cases, that's what you care about.
Agree!!Thanks