I am doing RNA-Seq data analysis using data available on SRA (SRR6047326). Samples are from the pig and they are paired. So, I did alignment using HISAT2.I used pig reference genome available on ENSEMBL. Then, I sorted .sam output files obtained from HISAT2. Now I am trying to get the count file using sorted sam file and it is not working. I am using annotated GTF file for pig which I downloaded from Ensembl.

I get count files using following commands, but all the counts are 0.

htseq-counts -i no myfile.sam genes.gtf>output.count

What is wrong here ? Are there any alternatives to htseq-count?

  • 2
    $\begingroup$ Help us to help you. What does "not working" mean? Please post relevant error messages, if any. $\endgroup$ – fridaymeetssunday Apr 22 '18 at 8:57
  • $\begingroup$ I runs well and gives me the count file. But there are zero counts for all the genes. At the end of count files it scores most of the genes are on --alignment_not_unique. $\endgroup$ – L R Joshi Apr 22 '18 at 15:31
  • $\begingroup$ I tried kallisto and it seems to work fine with kallisto.I can see counts for the multiple genes in tsv file generated by kallisto. $\endgroup$ – L R Joshi Apr 22 '18 at 15:32
  • $\begingroup$ Is that the actual command you ran? Because the "no" in -i no myfile.sam is odd, and since kallisto is returning some reasonable counts we can assume the issue with the htseq command and not the alignments (which was my 1st thought). $\endgroup$ – fridaymeetssunday Apr 22 '18 at 18:07
  • $\begingroup$ Speaking of kallisto, it is also fine to use those results for your analysis, unless that there a specific reason to use htseq-count. $\endgroup$ – fridaymeetssunday Apr 22 '18 at 18:08

I m not sure about htseq-counts but I would suggest go for featureCounts, I will update my answer and share the script/command. Give it a try

the command i use

featureCounts -T 50 -C -p -t exon -g gene_id -a /run/media/data2/gencode.v21.annotation.gtf  -o PRIMARY_CELL.txt $(ls *.bam | sort )
  • $\begingroup$ Thanks. I haven't use featureCounts. I found that it could run on R. Hopefully it works. $\endgroup$ – L R Joshi Apr 22 '18 at 5:28
  • 1
    $\begingroup$ yes you can but it takes lot of memory $\endgroup$ – krushnach Chandra Apr 22 '18 at 5:52
  • $\begingroup$ I am running on cluster. But, I have a question here. I have .sam files obtained from HISAT2 which were generated from the paired fastq files. Can I use that .sam files directly for feature counts ? Do I need to sort the sam file and then use feature counts ? $\endgroup$ – L R Joshi Apr 22 '18 at 15:36
  • 1
    $\begingroup$ I used the following commands fcounts <- featureCounts(files=“reads.sam”,genome="genes.gtf"). And my results were: Unassigned_NoFeatures = 5305153, and Unassigned_Unmapped=49253221. So, none of the genes were mapped. Any idea what could have gone wrong here ? $\endgroup$ – L R Joshi Apr 22 '18 at 15:46
  • 1
    $\begingroup$ Yes, I found a lot of papers using STAR aligner. I will give it a try. Thanks. $\endgroup$ – L R Joshi Apr 23 '18 at 23:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.