Entering edit mode
6.4 years ago
mccormack
▴
90
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
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.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 ?
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.
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.
Maybe they are stranded?
Please add the commands you used to obtain this result.
This is the command for S1.
That will print the output to the terminal right? You miss the part where you redirect the output to a file.