Question: HTSeq special counters occur at top and bottom of output file
0
gravatar for mccormack
10 months ago by
mccormack20
United States
mccormack20 wrote:

I have a bunch of HTSeq output files that have the special counters at the top and bottom of the file, and the numbers for each are different. For example, at the top of a file:

__no_feature    5272304
__ambiguous 0
__too_low_aQual 0
__not_aligned   0

And, at the bottom of the file:

__no_feature    308479
__ambiguous 745
__too_low_aQual 0
__not_aligned   0
__alignment_not_unique  9422976

Has anyone seen this before and have an explantation ?

This is the command:

htseq-count -f bam --stranded=no   /PHShome/data/HT/S1.bam  /apps/lab/Brassica_napus.AST_PRJEB5043_v1.35.gtf

This will print to an output file.

This is run on a cluster and in the script, see below, is a command, BSUB -o HTseq_S1_HTSeq.txt, to send the output to a text file named HTseq_S1_HTSeq.txt.

The whole script would be:

!/bin/bash   
module use /apps/modulefiles/test   
module load python/2.7.3   
module load numpy/1.9   
module load matplotlib/1.4.3   
module load htseq/0.6.1   
module load pysam/0.9.1.4   
cd /PHShome/data/HT   
BSUB -L /bin/bash   
BSUB -J HTseq_S1   
BSUB -sla e_sc   
BSUB -q big-multi   
BSUB -o HTseq_S1_HTSeq.txt  
BSUB -e HTseq_S1_HTSeq.err  
BSUB -N  
htseq-count -f bam --stranded=no   /PHShome/data/HT/S1.bam /apps/lab/Brassica_napus.AST_PRJEB5043_v1.35.gtf
htseq • 495 views
ADD COMMENTlink modified 10 months ago • written 10 months ago by mccormack20
1

I strongly suspect there is a >> in your command, which would append the results and not overwrite. Please put the whole command including the way you direct it to a file.

ADD REPLYlink written 10 months ago by Macspider2.6k

The BSUB -o command sends the output to the file named HTseq_S1_HTseq.txt. But, even if it were appending, why would it append only the special counters and not the count data ?

ADD REPLYlink written 10 months ago by mccormack20
1

Not answering your question but wanted to make a suggestion. Can you try featureCounts to do the counting? It is very fast and should hopefully not produce any errors.

ADD REPLYlink written 10 months ago by genomax57k

We could do this, but the HTSeq was fine for 16 of the 24 samples. It was only 8 samples that had this problem and I do not see any difference in the commands (other than file names) between samples that worked and samples that resulted in the special counters at the top and bottom of the output.

ADD REPLYlink written 10 months ago by mccormack20

Maybe they are stranded?

ADD REPLYlink written 10 months ago by Macspider2.6k

Please add the commands you used to obtain this result.

ADD REPLYlink written 10 months ago by WouterDeCoster32k

This is the command for S1.

htseq-count -f bam --stranded=no   /PHShome/data/HT/S1.bam    /apps/lab/Brassica_napus.AST_PRJEB5043_v1.35.gtf
ADD REPLYlink modified 10 months ago by genomax57k • written 10 months ago by mccormack20
1

That will print the output to the terminal right? You miss the part where you redirect the output to a file.

ADD REPLYlink written 10 months ago by WouterDeCoster32k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 772 users visited in the last hour