Has anyone with experience using both HISAT2 and STAR ever observed a large difference (0 vs 400+) in the counts for a gene depending on which aligner they used?
I processed some RNA-seq data first using STAR and used quantmode to get the counts. Afterwards I checked a few genes that are known to be regulated by the target protein of interest that was inhibited with a drug. Most of the genes looked fine (i.e. down in treated and up in vehicle). But, one gene called IFI27 had a 0 count across all samples, I confirmed this result using several other counting programs. I then aligned using hisat2 and this produced the expected result for IFI27 with more than 400 counts for the vehicles and 17 counts for the treated samples. Other genes had small differences but nothing close to the difference for IFI27.
So to me it makes sense to use hisat2 for my mapping. However, I am bothered by the fact that STAR is highly regarded but somehow could not handle this specific gene. I was curious to see if anyone else has encountered this or has any thoughts on what might be happening.
STARindex Command
module load star/2.7.3a
STAR --runThreadN 23 --runMode genomeGenerate --genomeDir /home --genomeFastaFiles /home/RNA/GRCh38.p13_genomic.fna --sjdbGTFfile /home/RNA/GRCh38.p13_genomic.gtf --sjdbOverhang 149
STARalign Command
STAR
--outSAMattributes All
--outSAMtype BAM SortedByCoordinate
--quantMode GeneCounts
--readFilesCommand gunzip -c
--runThreadN $Threads
--sjdbGTFfile $GTFFILE
--outReadsUnmapped Fastx
--outMultimapperOrder Random
--outWigType wiggle
--genomeDir $GENOMEDIR
--readFilesIn ${INPUTDIR}/${OUTPREFIX}_1.fq.gz ${INPUTDIR}/${OUTPREFIX}_2.fq.gz
--outFileNamePrefix $OUTPREFIX
HISAT2index Command
hisat2-build -p 23 -f /home/refFiles/GRCh38.p13_genomic.fna GRCh38
HISAT2align Command
hisat2 -p 8 -q -x $GENOMEDIR -1 ${INPUTDIR}/${OUTPREFIX}_1.fq.gz -2 ${INPUTDIR}/${OUTPREFIX}_2.fq.gz |