...

for r1R1 in Raw_Reads/*_R1_*001.fastq.gz; do r2R2=$(echo $r1$R1| sed 's/_R1_/_R2_/'); name=$(echo $r1$R1|sed 's/_R1_001.fastq.gz//'|sed 's/Raw_Reads\///'); echo "fastp PE $r1 $r2-i $R1 -I $R2 -baseouto Trim_Reads/$name.trim.R1.fastq.gz ILLUMINACLIP:TruSeq3-PE-2.fa:4:30:10 MINLEN:30";done > trim_commands

for R1 in Raw_Reads/*_R1_001-O Trim_Reads/$name.trim.R2.fastq.gz; do R2=$(echo $R1| sed 's/_R1_/_R2_/'); name=$(echo $R1|sed 's/_R1_001.fastq.gz//'|sed 's/Raw_Reads\///'); echo "fastp -i $R1 -I $R2 -o Trim_Reads/$name.trim.R1.fastq.gz -O Trim_Reads/$name.trim.R2.fastq.gz -w 4 --detect_adapter_for_pe -j 04_Trim_Logs/$name.json -h 04_Trim_Logs/$name.html &> 04_Trim_Logs/$name.log.txt";done > fastp.commands

Use the head and wc -l to see what the output is and how much there is of it respectively.

Again as we discussed in the multiqc tutorial, running this number of commands is kind of a boarder line case, there are not a lot of total reads, but there are a large number of samples and our command does request additional processors so we should not be on the head node. 

-w 4 --detect_adapter_for_pe -j Trim_Logs/$name.json -h Trim_Logs/$name.html &> Trim_Logs/$name.log.txt";done > fastp.commands

Use the head and wc -l to see what the output is and how much there is of it respectively.

Again as we discussed in the multiqc tutorial, running this number of commands is kind of a boarder line case of being better to run in idev or as a submitted job, there are not a lot of total reads, but there are a large number of samples. Because our command does request additional processors we should not run on the head node, if we wanted to use only a single processor, the job would take so long to run on the head node that even waiting in the queue system would be faster.

Do only one of the following, but do read through both options as there are different discussions about the process in each.

Submitting as a job

cp /corral-repl/utexas/BioITeam/gva_course/GVA2022.launcher.slurm trim.slurm
nano trim.slurm

Again while in nano you will edit most of the same lines you edited in the in the breseq tutorial. Note that most of these lines have additional text to the right of the line. This commented text is present to help remind you what goes on each line, leaving it alone will not hurt anything, removing it may make it more difficult for you to remember what the purpose of the line is

The changes to lines Line 17 being set to -n 17 allows 17 jobs to run at the same time, since our command uses -w 4 (4 threads) this job will use all 68 threads available. The changes to lines 22 and 23 are optional but will give you an idea of what types of email you could expect from TACC if you choose to use these options. Just be sure to pay attention to these 2 lines starting with a single # symbol after editing them.

...

Again use ctl-o and ctl-x to save the file and exit.

sbatch trim.slurm

The job should take less than 10 minutes once it starts if everything is working correctly, the showq -u command can be used to check for the job finishing.

Evaluating the output

ls Trim_Reads | wc -l
grep -c "TrimmomaticPE: Completed successfully" Queue_job.o*
grep -c "100.00%" Queue_job.o*

The above 3 commands are expected to show 1088 and 272 and 0 respectively, if you see other answers it suggests that something went wrong with the trimming command itself. If so remember I'm on zoom if you need help looking at whats going on. The most common error that I expect will be that you ran out of ram by trying to process too many samples at once (such as if you used a -n 68 option in your .slurm file). 

Beyond the job finishing successfully, the best way to evaluate this data would actually be to move back to the multiqc tutorial and repeat the analysis there that was done on the raw files for the trimmed files here.

Personal experience

The for loop above focuses just on generating the trim commands. In my experience that is only half of the job, the other half is capturing the individual outputs so you can evaluate how many reads were trimmed in what way for each sample. Perhaps you will find this command helpful in your own work: 

mkdir trimLogs; mkdir Trim_Reads; for r1 in Raw_Reads/*_R1_*.fastq.gz; do r2=$(echo $r1|sed 's/_R1_/_R2_/'); name=$(echo $r1|sed 's/_R1_001.fastq.gz//'|sed 's/Raw_Reads\///'); echo "trimmomatic PE $r1 $r2 -baseout Trim_Reads/$name.fastq.gz ILLUMINACLIP:TruSeq3-PE-2.fa:4:30:10 MINLEN:30 >& trimLogs/$name.log"; done > trim_commands

After the job completes, the following command is useful for evaluating its success:

echo -e "\nTotal read pairs:";wc -l trim_commands;echo "Successful trimmings:"; tail -n 1 trimLogs/*|grep -c "TrimmomaticPE: Completed successfully";echo "Potential trimming errors:"; cat trimLogs/*|grep -c "100.00%"

This gives me a quick readout of if all of the individual commands finished correctly.

...

mkdir Trim_Reads Trim_Logs
sbatch trim.slurm

The job should take less than 10 minutes once it starts if everything is working correctly, the showq -u command can be used to check for the job finishing.

Running on idev

The alternative to running all the commands as a submitted job is of course to run on an idev node. In the multiqc tutorial, we were able to tell fastqc to analyze all samples just by using the "*" wildcard as the only required input to fastqc is the filename. Here our command is much more intricate which may seem like it precludes us from being able to run interactively as we never would type 272 nearly identical commands. Obviously there is a trick for this.

chmod +x <filename>
./<filename>

Before we jump to making our commands file executable and executing it, we want to change it to be slightly different. Specifically, above we used -w 4 to specify we wanted to use 4 processor for every command. While this worked great when we also were launching 17 processes at the same time as it used all 68 processes, when executing a commands file from the command line without the help of the queue system, only 1 sample at a time will launch so you likely think we need to increase to 68 processors. 

#1 change the for loop:
for R1 in Raw_Reads/*_R1_001.fastq.gz; do R2=$(echo $R1| sed 's/_R1_/_R2_/'); name=$(echo $R1|sed 's/_R1_001.fastq.gz//'|sed 's/Raw_Reads\///'); echo "fastp -i $R1 -I $R2 -o Trim_Reads/$name.trim.R1.fastq.gz -O Trim_Reads/$name.trim.R2.fastq.gz -w 68 --detect_adapter_for_pe -j Trim_Logs/$name.json -h Trim_Logs/$name.html &> Trim_Logs/$name.log.txt";done > fastp.commands

#2 use sed to do an in file replacement (something new you haven't seen before in this class)
sed -i 's/ -w 4 / -w 68 /g' fastp.commands

Note that if you use the sed command above, you need to be very careful in what you choose to match to. If you just choose "4" and replace with "68" then the commands file will then change any file name that has a 4 into 68 and all those samples will fail. When using sed to do replacements, always make sure you have a unique handle, when you don't, and when you don't need one.

chomd +x fastp.commands
./fastp.commands

Once the command is started continue reading below.

Comparing different run options

In previous years it has been common to question what the fastest way of getting a large set of samples analyzed is with respect to threads and Nodes and tasks. Here we hav an opportunity to do just that, and have some surprising results. Since we have been working with idev sessions all along we'll start with the following:

This brings us to our first question: how can using fewer processors speed things up. I am not completely sure but I do have some hypotheses, and expect the answer to be somewhere in the middle:

1. We have so few reads (10,000s) that trying to spread them in 68 different groups, read them, modify them, write them back to the growing trimmed files, perform statistics on all the reads. That the computation required to spread them out (the overhead) actually exceeds whatever speed bonus is available with the extra threads.

2. If you used this command and look in the .log.txt file that was generated thanks to our pipe "&>", you will find a interesting note:

   No Format
   WARNING: fastp uses up to 16 threads although you specified 68

   this means that despite requesting fastp use 68 threads, it is only capable of, and only used 16 so there is an even smaller difference in number of threads available. Unfortunately this is something that doesn't seem to be documented on their website or in the help information.

Given that threads only effect the speed of a single sample I also attempted to trim all the read files at once by telling the idev node to run them in the background (this was done by adding a single & as the last character on each line of the commands file)

While both of these errored after ~1/5 and ~1/2 of the total samples making them a non viable single pass method, they did finish in <10 seconds. Using grep we could generate a list of samples which did finish, store those in a file, then use another grep command with the -v and -f options to remove the finished samples from analysis. This is something that can be very helpful if you have a large number of samples and your run timed out as you can focus on just samples that have not yet finished instead of having to rerun all samples, or manually edit the commands file to remove samples which finished. Since this run failed due to lack of memory, not lack of time its less reasonable to try to tackle it in this manner.

That brings us to our last set of runs: those working with the commands as a submitted job.

Based on what we have already seen, it is probably not surprising that using 68 (really 16) threads and only evaluating 1 sample at a time took approximately the same amount of time as it did when running on an idev node as those conditions are functionally equivalent. Whaat may be surprising is the lack of improvement despite running 4x more samples at the same time. Again hypothesies:

1. The amount of overhead the job launcher has in checking when a command finishes and starting the next command is similar to the amount of overhead in splitting reads into 4 sets
2. The biggest contributor to time is not something that can be improved with additional threads

Evaluating the output

The first thing you always want to check when working with a lot of commands simultaneously is if they all finished correctly (note above how running all 272 jobs in the background at basically the same time did not finish correctly). Typically this is where the log files we generate with "&>" come in handy. if you look at the tail of any 1 of the files you are likely to see something like:

Duplication rate: 6.5207%

Insert size peak (evaluated by paired-end reads): 80

JSON report: Trim_Logs/E2-1_S189_L001.json
HTML report: rim_Logs/E2-1_S189_L001.html

fastp -i Raw_Reads/E2-1_S189_L001_R1_001.fastq.gz -I Raw_Reads/E2-1_S189_L001_R2_001.fastq.gz -o Trim_Reads/E2-1_S189_L001.trim.R1.fastq.gz -O Trim_Reads/E2-1_S189_L001.trim.R2.fastq.gz -w 68 --detect_adapter_for_pe -j Trim_Logs/E2-1_S189_L001.json -h Trim_Logs/E2-1_S189_L001.html 
fastp v0.23.2, time used: 3 seconds

Typically what you should look for is some kind of anchor that you can pass to grep that is as far down the file as possible. Sometimes you will be lucky and the program will actually print something like "successfully complete". In our case the last line looks promising, "fastp v0.23.2, time used:" seems likely to be printed as the last step in the program.

ls Trim_Logs/*.log.txt|wc -l
wc -l fastp.commands
tail -n 1 Trim_Logs/*.log.txt|grep -c "^fastp v0.23.2, time used:" 

The above 3 commands are all expected to return 272. If so remember I'm on zoom if you need help looking at whats going on. 

Beyond the job finishing successfully, the best way to evaluate this data would actually be to actually use it. That could me running it through multiqc to see how the trimming has effected overall quality, or assembling the reads, or mapping them to a reference (not applicable in this case since they are from a variety of sources and you have not been provided reference files.

Optional next steps:

1. Consider moving back over to the multiqc tutorial use multiqc to determine well the trimming worked. 
2. The reads could then further be assembled in the Spades tutorial as they are all plasmid samples.

...