Help with HTSeq-Count

Question

[emmachristison]

Hello,

This script will not work and I have no idea why, please help!

Thanks!

module add apps/STAR/2.5.2b
module add apps/samtools/0.1.19
module add apps/python27
module add libs/python/numpy
module add libs/python/matplotlib
module add apps/HTSeq/0.6.1

sample_dir=~/WORKING_DATA/Blakeley/Blakeley/ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP055/SRP055810/fq/STAR_out_3_Blakeley

samples=$( ls *.bam)
hg38=~/WORKING_DATA/genome_dir/bt2

for j in samples ; do

 echo "HTseq-count of $j"
        samtools view  ${j}fastqAligned.out.bam | htseq-count --stranded=no - 

${hg38}/Homo_sapiens.GRCh38.85.gtf > ${j}.counts.txt

done;

wait

Answer 1

[Michael]

Normally, I would not be able to help, because you forgot the error message which is essential. In this case it is easy to fix:

for j in samples ; do should be for j in $samples ; do

[emmachristison]

What a stupid mistake!

Thank you very much Michael, I'll try that!

[emmachristison]

It worked..kind of.....it went through them but all the output files are empty, does anyone know why?

[GenoMax]

Before going through a loop like this did you ensure that you are able to run one file (with the command line options that are created by the loop?). Fish one out of command line out of the logs and examine it. You may be able to spot the malformation (if one exists) even before running it standalone to diagnose what could be going wrong.