Newbie question here, sorry. I am trying to generate a counts matrix for paired-end non-stranded data.  I feel pretty good about the steps I have taken so far, but this is the first time I have used this particular type of data.  The counts matrix I have generated has allele specific reads which I wasn't anticipating, and the counts seem vary off, especially between the two alleles.  When generating a counts matrix with allele specific information, is there any special steps I need to consider?  Also is there a way to get non-allele specific read counts?

I think a potential issue maybe with samtools or ht-seq.  Feedback appreciated.  I am willing to switch to bioconductor/R tools also.  I don't have too much experience outside with RNA-seq outside of these programs though. 

    samtools sort -n accepted_hits.bam  sample_sorted
    samtools view sample_sorted.bam >sample.sam
    
    htseq-count -s no --idattr=Parent sample.sam organism.gff > sample.counts

What would be helpful for me to know:
Does it look like I'm on the right track?  If yes, why are my allele reads so different?  If no, what steps did I overlook, and do you have suggestions?
WARNING! Phased genome :)

As you have correctly noticed, C. albicans is one of few species which has a phased genome, i.e. sequence from both parental chromosomes (but keep in mind that in particular case of C. albicans we don't know if all A chromosomes are from one parental, and B chromosomes from another, here is the [paper][1] that generated this assembly is you are curious). 

To my knowledge, today most of the software and pipelines (especially mappers and variant callers) are generally not well adapted in dealing with phased data.  So if you are not really interested in allele-specific expression, its better to use reference genome and annotations from one parental (just subset fasta file and gff). If you are interested in allele-specific expression and will use phased genome, then you should expect to see a huge number of multimapped reads, because parental are very similar in DNA sequence. To avoid this at some extent you can restrict the number of allowed mismatches in settings of your mapper for example to 1 or 0 (in latter case you will miss the reads with sequencing errors though). In STAR mapper, these options are `--outFilterMismatchNmax`, `--outFilterMismatchNoverLmax`. More details are in [manual][2]. 

And to speed up things, I would recommend you to use STAR mapper with option `--quantMode GeneCounts`, it will automatically generate count table of reads identical to HTSeq. 


  [1]: https://genomebiology.biomedcentral.com/articles/10.1186/gb-2013-14-9-r97
  [2]: https://github.com/alexdobin/STAR/blob/master/doc/STARmanual.pdf