Linux for loop for 2 inputs and 4 outputs for Trimmomatic fastq quality trimming

Question

[Rashedul Islam]

I need help to write a for loop to run Trimmomatic tool for quality trimming of paired end fastq files.

Input: file01_R1.fastq, file01_R2.fastq.

The first file contains first read of paired end reads and the second file contains second reads. I have 100 input files in the same directory (50 files contain *e.R1g.fastq and another 50 files, file02_R1.R2fastq, file02_R2.fastq suffixes) and so on.

Output: file01_R1_PE.fastq, file01_R1_SE.fastq, file01_R2_PE.fastq, file01_R2_SE.fastq

PE = paired end
SE = single end

I need to write a for loop so that I can run an executable for all 100 files. Any help please! Thanks!

[Alex Reynolds]

Hello Rashedul Islam!

We believe that this post does not fit the main topic of this site.

This is not a question about bioinformatics, but about programming loops. You may find answers on Stack Overflow. I'd recommend also specifying which language or shell you are working with.

For this reason we have closed your question. This allows us to keep the site focused on the topics that the community can help with.

If you disagree please tell us why in a reply below, we'll be happy to talk about it.

Cheers!

[Michael]

Can we reopen this one? It's still a beginner question but relation to bioinformatics is specific enough. Running Trimmomatic on both PE and SE data has its own complications (need to differentiate PE from SE and chose correct adapter set).

[Deepak Tanwar]

Opening this question again for answers would be good.

[Rashedul Islam]

I have edited my question. Please give a look.

Answer 1

[Giovanni M Dall'Olio]

Don't use for loops in bash, or possibly in any programming language.

Have a look at GNU/parallel:

find yourfastqfolder -iname 'file*_R1.fastq' | parallel -j 4 'trimmomatic_command {} {.}_R2.fastq'

This will launch all the jobs in parallel, in bunches of 4. Have a look at Gnu Parallel - Parallelize Serial Command Line Programs Without Changing Them for more examples.

[cschu181]

Thanks, this might come in handy. I am wondering, though, what is wrong with for loops? Care to elaborate?

Answer 2

[Rashedul Islam]

Finally, I used this code to run Trimmomatic in a loop in linux.

for file in /path/to/*_R1.fastq
do
withpath="${file}"
filename=${withpath##*/}
base="${filename%*_*.fastq}"
echo "${base}"
java -jar trimmomatic-0.33.jar PE /path/to/"${base}"*_R1.fastq /path/to/"${base}"*_R2.fastq /path/to/"${base}"_R1.trimmed_PE.fastq /path/to/"${base}"_R1.trimmed_SE.fastq /path/to/"${base}"_R2.trimmed_PE.fastq /path/to/"${base}"_R2.trimmed_SE.fastq ILLUMINACLIP:/path/to//Trimmomatic-0.33/adapters/TruSeq3-PE-2.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36 HEADCROP:14
done

[bioinformaticssrm2011]

I know its a bit late. But above script won't be able to distinguish between lets say sample1_R1.fastq and sample12_R1.fastq.

Above script might get confused whether to select sample1 and sample12.

[Michael]

I think that you are now running all the files - SE and PE - through PE processing with PE adapter trimming only, but the the SE reads should be trimmed using SE processing using TruSeq-SE.fa adpater file. Therefore you should build in a conditional processing step depending on if the filename contains PE or SE.

Answer 3

[pmiller]

This is a late response but for anyone else who gets stuck, I finally found a way to do this for myself. The use of Gnu/parallel with "find" was a big fail for me - would not run. Then I finally realized after much searching that Trimmomatic does not have any kind of batch mode and that was the problem. I don't know why I did not understand this, but I found a workaround with a loop that is awesome at - http://www.datacarpentry.org/2015-08-24-ISU/lessons/06-readQC.html. I used it for a different tool in Trimmomatic (CROP) and it worked beautifully.

$ for infile in *.fastq
do
    outfile=$infile\_trim.fastq
    java -jar ~/Trimmomatic-0.32/trimmomatic-0.32.jar SE $infile $outfile SLIDINGWINDOW:4:20 MINLEN:20
done