This is a test version of Biostars. For the public version, visit https://www.biostars.org.
looping through PE files ?

HI

Newbie to shell scripting, written a bash script loop through 100 PE files, unable to generate output. Can some one correct my code or any suggestions. 

#!/bin/bash
for f1 in *_R1_sic.fastq
do
  f2=${f1%%_R1.sic.fastq}"_R2.sic.fastq"
  /home/cent/anaconda3/envs/shannon_env/bin/shannon_cpp shannon -k 30 -l 150 -p /home/cent/data/cat_samples/fastqgz/trim/cent1/sickle_files/$f1 /home/cent/data/cat_samples/fastqgz/trim/cent1/sickle_files/$f2 -o /home/cent/data/cat_samples/fastqgz/trim/cent1/sickle_files/$f1$f2_shan -t 5 -u 5 -g 150 -m 150G > /home/cent/data/cat_samples/fastqgz/trim/cent1/sickle_files/$f1$f2_shan_log
done
rna-seq
• 2,523 views • link
updated by antonioggsousa antonioggsousa • written by sunnykevin97 sunnykevin97

unable to generate output

Use echo at start of the command /home/cent/anaconda3/envs/shannon_env/bin/shannon_cpp shannon to make sure command lines look ok. If not you should get an idea of what needs to be changed.

A better option is to grab sample names by doing:

name = $(basename $f1 _R1_sic.fastq)
# then reconstruct names as needed
f1=${name}_R1_sic.fastq
f2=${name}_R2_sic.fastq
${name}.log
• 0 views • link
written by GenoMax GenoMax 
for f1 in *_R1_sic.fastq
do
echo "/home/cent/anaconda3/envs/shannon_env/bin/shannon_cpp shannon -k 30 -l 150 -p /home/cent/data/cat_samples/fastqgz/trim/cent1/sickle_files/$f1 /home/cent/data/cat_samples/fastqgz/trim/cent1/sickle_files/${f1/_R1_/_R2_} -o /home/cent/data/cat_samples/fastqgz/trim/cent1/sickle_files/$f1${f1/_R1_/_R2_}_shan -t 5 -u 5 -g 150 -m 150G > /home/cent/data/cat_samples/fastqgz/trim/cent1/sickle_files/$f1${f1/_R1_/_R2_}_shan_log;"
done

Remove echo after checking the function.

However, your output directory and logs will have unnecessarily long names. Assuming that base name of the files is unique, you can try this:

for f1 in *_R1_sic.fastq
do
echo "/home/cent/anaconda3/envs/shannon_env/bin/shannon_cpp shannon -k 30 -l 150 -p /home/cent/data/cat_samples/fastqgz/trim/cent1/sickle_files/$f1 /home/cent/data/cat_samples/fastqgz/trim/cent1/sickle_files/${f1/_R1_/_R2_} -o /home/cent/data/cat_samples/fastqgz/trim/cent1/sickle_files/${f1%_R1_*}_shan -t 5 -u 5 -g 150 -m 150G > /home/cent/data/cat_samples/fastqgz/trim/cent1/sickle_files/${f1%_R1_*}_shan_log;"
done

Remove echo for execution after checking the function.

I suggest you to use parallel for this kind of work.

• 0 views • link
written by cpad0112 cpad0112

It not working. : No such file or directory I give all the input files. If I run interactively it works.

• 0 views • link
written by sunnykevin97 sunnykevin97

1 answer

Hi,

From your example, there is one inconsistency, the pattern name of the file is _R1_sic.fastq or _R1.sic.fastq?

I would suggest to use bash variables to define the path to the program and to the files. It simplifies your code. In addition you can use escape, i.e., \ to break down your code in multiple lines. It also helps to read the code.

The problem with your code is related with the f2=${f1%%_R1.sic.fastq}"_R2.sic.fastq" that is not doing what you're expecting.

Well assuming that is the second one, i.e., _R1.sic.fastq (if not, change accordingly), therefore do the following:

#!/bin/bash

PROGRAM=/home/cent/anaconda3/envs/shannon_env/bin/shannon_cpp
FILES=/home/cent/data/cat_samples/fastqgz/trim/cent1/sickle_files

for f1 in *_R1.sic.fastq
do
        f2=${f1%%_*}_R2.sic.fastq
        $PROGRAM shannon -k 30 -l 150 \
                -p $FILES/$f1 $FILES/$f2 \
                -o $FILES/${f1%%_*}_shan \
                -t 5 -u 5 -g 150 -m 150G > $FILES/${f1%%_*}_shan_log
done

Let me know if the code worked.

António

• 0 views • link
written by antonioggsousa antonioggsousa

When I run the script it overwrites the R1 files. All the R1 replaced with text files. 

-rwxr-xr-x 1 cent cent  926 Jul 27 23:10 10_R1.fastq
-rwxr-xr-x 1 cent cent 5.2G Jul 26 20:20 10_R2.fastq
-rwxr-xr-x 1 cent cent  926 Jul 27 23:10 11_R1.fastq
-rwxr-xr-x 1 cent cent 5.4G Jul 26 20:14 11_R2.fastq
-rwxr-xr-x 1 cent cent  926 Jul 27 23:10 12_R1.fastq
-rwxr-xr-x 1 cent cent 5.9G Jul 26 20:18 12_R2.fastq
-rwxr-xr-x 1 cent cent  926 Jul 27 23:10 13_R1.fastq
-rwxr-xr-x 1 cent cent 4.9G Jul 26 20:17 13_R2.fastq
-rwxr-xr-x 1 cent cent  926 Jul 27 23:10 14_R1.fastq
-rwxr-xr-x 1 cent cent 5.9G Jul 26 20:18 14_R2.fastq
-rwxr-xr-x 1 cent cent  926 Jul 27 23:10 15_R1.fastq
-rwxr-xr-x 1 cent cent 6.5G Jul 26 20:10 15_R2.fastq
-rwxr-xr-x 1 cent cent  926 Jul 27 23:10 16_R1.fastq
-rwxr-xr-x 1 cent cent 6.6G Jul 26 20:22 16_R2.fastq
-rwxr-xr-x 1 cent cent  926 Jul 27 23:10 17_R1.fastq
-rwxr-xr-x 1 cent cent 6.5G Jul 26 20:15 17_R2.fastq
-rwxr-xr-x 1 cent cent  926 Jul 27 23:10 18_R1.fastq
• 0 views • link
written by sunnykevin97 sunnykevin97

Hi,

Yes, I understand why. This is because I did not noticed in your code that $f1$f2_shan is actually the same as $f1, that's why is overwriting the R1 files. I corrected the code to appear in the output only the name of the sample: ${f1%%_*}, i.e., the pattern before the first _ underscore.

#!/bin/bash

PROGRAM=/home/cent/anaconda3/envs/shannon_env/bin/shannon_cpp
FILES=/home/cent/data/cat_samples/fastqgz/trim/cent1/sickle_files

for f1 in *_R1.sic.fastq
do
        f2=${f1%%_*}_R2.sic.fastq
        $PROGRAM shannon -k 30 -l 150 \
                -p $FILES/$f1 $FILES/$f2 \
                -o $FILES/${f1%%_*}_shan \
                -t 5 -u 5 -g 150 -m 150G > $FILES/${f1%%_*}_shan_log
done

Let me know if this solves your problem,

António

• 0 views • link
written by antonioggsousa antonioggsousa

It works great. I make similar script for preprocesing analysis any suggestions.

for f1 in *_R1.fastq
  do
          f2=${f1%%_*}_R2.fastq
          time java -Xmx12G -jar $PROGRAM PE \
                 $FILES/$f1 $FILES/$f2 \ 
                 $FILES/${f1%%_*}.trim.fastq  $FILES/${f1%%_*}.untrim.fastq \
                 $FILES/${f2%%_*}.trim.fastq  $FILES/${f2%%_*}.untrim.fastq \
                 SLIDINGWINDOW:4:20 MINLEN:20 ILLUMINACLIP:/home/cent/data/softwares/Trim/adapters/TruSeq3-PE.fa:2:30:15 > $FILES/${f1%%_*}_trim_log
done
• 0 views • link
written by sunnykevin97 sunnykevin97

I don't have further suggestions. The code seems good to me.

I just want to add that I changed the answer above, to correct the problem related with the overwriting of R1 files.

António

• 0 views • link
written by antonioggsousa antonioggsousa

Log in to answer this question.

More posts like this