Running mutltiple Fastq files with Hisat2 in a loop

Question

[n.tear]

I am trying to loop paired end reads in hisat2 with an output of sorted bam using a pipe.

I am a novice at coding but am keen to improve my skills by learning how to loop.

I have created this loop using the answers from Run mutltiple Fastq files on Hisat2

 for f in `ls *.fastq.gz | sed 's/_[12].fastq.gz//g' | sort -u`
do
echo hisat2 -x /mnt/z/NathanHaffordTear/RNAseqproject/RNAseqanalysis/hg38_index -1 ${f}_1.fastq.gz -2 ${f}_2.fastq.gz | samtools sort -o ${f}.hisat.sorted.bam
done

For some reason the loop does work without the pipe but not with. How can i configure this loop to work with the piping as I do not want to have the intermediate files?

Many thanks in advance

[russhh]

Is the echo supposed to precede the call to hisat2?

[n.tear]

I added it in to check the code before running.

[cpad0112]

$ ls *.gz
test_1.fastq.gz  test_2.fastq.gz

--

for f in `ls *.fastq.gz | sed 's/_[12].fastq.gz//g' | sort -u`
do
echo "hisat2 -x /mnt/z/NathanHaffordTear/RNAseqproject/RNAseqanalysis/hg38_index -1 ${f}_1.fastq.gz -2 ${f}_2.fastq.gz | samtools sort -o ${f}.hisat.sorted.bam"
done

output:

hisat2 -x /mnt/z/NathanHaffordTear/RNAseqproject/RNAseqanalysis/hg38_index -1 test_1.fastq.gz -2 test_2.fastq.gz | samtools sort -o test.hisat.sorted.bam

OP is supposed to run:

for f in `ls *.fastq.gz | sed 's/_[12].fastq.gz//g' | sort -u`
do
echo "parising f...wait..." && hisat2 -x /mnt/z/NathanHaffordTear/RNAseqproject/RNAseqanalysis/hg38_index -1 ${f}_1.fastq.gz -2 ${f}_2.fastq.gz | samtools sort -o ${f}.hisat.sorted.bam 
done

OP can also try parallel:

parallel --plus --dry-run  "hisat2 -x /mnt/z/NathanHaffordTear/RNAseqproject/RNAseqanalysis/hg38_index -1 {} -2 {=s/1/2/=} | samtools sort -o {=s/_1.fastq.gz//=}.hisat.sorted.bam" ::: *_1.fastq.gz

output:

hisat2 -x /mnt/z/NathanHaffordTear/RNAseqproject/RNAseqanalysis/hg38_index -1 test_1.fastq.gz -2 test_2.fastq.gz | samtools sort -o test.hisat.sorted.bam

Answer 1

[ATpoint]

Remove the echo as you want to execute rather than print the command. If you want "status updates" do

echo $f && hisat2 (...)

[n.tear]

Hi thanks for the reply

Im aware of this but want to make sure the code works with echo before running

running my code I get the error:

[W::sam_read1] Parse error at line 1
samtools sort: truncated file. Aborting
[W::sam_read1] Parse error at line 1
samtools sort: truncated file. Aborting
[W::sam_read1] Parse error at line 1
samtools sort: truncated file. Aborting
[W::sam_read1] Parse error at line 1
samtools sort: truncated file. Aborting
[W::sam_read1] Parse error at line 1
samtools sort: truncated file. Aborting
[W::sam_read1] Parse error at line 1
samtools sort: truncated file. Aborting
[W::sam_read1] Parse error at line 1
samtools sort: truncated file. Aborting
[W::sam_read1] Parse error at line 1
samtools sort: truncated file. Aborting
[W::sam_read1] Parse error at line 1
samtools sort: truncated file. Aborting
[W::sam_read1] Parse error at line 1
samtools sort: truncated file. Aborting
[W::sam_read1] Parse error at line 1
samtools sort: truncated file. Aborting
[W::sam_read1] Parse error at line 1
samtools sort: truncated file. Aborting
[W::sam_read1] Parse error at line 1
samtools sort: truncated file. Aborting
[W::sam_read1] Parse error at line 1
samtools sort: truncated file. Aborting
[W::sam_read1] Parse error at line 1
samtools sort: truncated file. Aborting
[W::sam_read1] Parse error at line 1
samtools sort: truncated file. Aborting
[W::sam_read1] Parse error at line 1
samtools sort: truncated file. Aborting
[W::sam_read1] Parse error at line 1
samtools sort: truncated file. Aborting

I will try with the added echo $f && hisat2 (...) and let you know how it goes

thanks

[ATpoint]

Of course because hisat2 does not produce output and therefore samtools does not get any valid stdin. Remove the echo and maybe put a head like hisat2 | head -n 1000 | samtools sort to get an idea if it works properly (and it should, the command is ok).

[n.tear]

Thanks that makes perfect sense. doh

The command you suggested before seems to be working great, i like getting the updates