Question: HTSeq special counters occur at top and bottom of output file
0
gravatar for mccormack
7 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 • 413 views
ADD COMMENTlink modified 7 months ago • written 7 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 7 months ago by Macspider2.4k

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 7 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 7 months ago by genomax52k

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 7 months ago by mccormack20

Maybe they are stranded?

ADD REPLYlink written 7 months ago by Macspider2.4k

Please add the commands you used to obtain this result.

ADD REPLYlink written 7 months ago by WouterDeCoster30k

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 7 months ago by genomax52k • written 7 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 7 months ago by WouterDeCoster30k
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: 1350 users visited in the last hour