This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Why are my bed graphs and HTseq-count files not agreeing?

When I load my bed graphs into IGV I see that 40/96 cells have signal on Xist. However when I look at my counts output from HTseq-counts for the same data I only see about 20 cells with Xist signal. Any ideas why this is?

My HTseq-count command:

~/bin/htseq-count \
        --format=bam \
        --order=name \
        --stranded=no \
        --minaqual=0 \
        --type=exon \
        --mode=union \
        --idattr=gene_id \
        {outdir}/{s}_aligned_sortedByName.bam \
        {GTF} \
    > {outdir}/{s}.counts

My Bedtools command:

~/bin/bedtools \
        genomecov \
        -split \
        -bg \
        -ibam {outdir}/{s}_aligned_sortedByCoord.bam \
        -g {CHROMINFO} \
    > {outdir}/{s}.bedgraph

I should mention that when I load the bed graphs into IGV I auto scale all tracks.

htseq bedgraphs

1 answer

Most likely the difference is due to htseq-count excluding multimappers and secondary alignments, which will occur regardless of what the minimum MAPQ is set to.

Log in to answer this question.