1
$\begingroup$

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?

Improve this question
$\endgroup$
6
  • 2
    $\begingroup$ Help us to help you. What does "not working" mean? Please post relevant error messages, if any. $\endgroup$
    – fridaymeetssunday
    Apr 22, 2018 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, 2018 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, 2018 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, 2018 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, 2018 at 18:08

1 Answer 1

Reset to default
2
$\begingroup$

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 )
Improve this answer
$\endgroup$
8
  • $\begingroup$ Thanks. I haven't use featureCounts. I found that it could run on R. Hopefully it works. $\endgroup$
    – L R Joshi
    Apr 22, 2018 at 5:28
  • 1
    $\begingroup$ yes you can but it takes lot of memory $\endgroup$
    – kcm
    Apr 22, 2018 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, 2018 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, 2018 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, 2018 at 23:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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