Trimmomatic is performed before genome assembly We have something called raw reads. If it is a paired-end sequencing, we have something called paired-end reads, consisting of an Read1 (R1) and read2 (R2) file. Can we directly assemble and reconstruct the genome from the raw read files? It is generally not advised to perform an assembly directly from these raw read files. These read files also contain low-quality reads that can not reliably be used to reconstruct the genome. These raw read files often are contaminated, therefore need to be subjected to a quality check up, before the assembly process. What is the major source of contamination? The major source of contamination is the sequencing adapters. Primer sequences can also be present within the reads. These adapter/ primer sequence-containing reads need to be processed for obtaining the valid read sequences. A FastQC analysis will generate a report, bearing information about the presence / absence of adapters in the R1 and R2 files. 💡 Here, you can see the adapter-ligated DNA fragment. In a SE sequencing, this blue-adapter-containing sequence will be read. In Paired-End sequencing, there is another round of cluster-generation. The second adapter-containing DNA fragment will be sequenced. How do adapter sequences get into the reads? After the library preparation, the library sequences will be loaded onto an Illumina Flow cell, where the adapter-ligated gDNA sequences (in case of de novo assembly) or RNA fragments (for transcriptome analysis) will be PCR-amplified, and further read. The library sequences will be read for a fixed number of cycles. In the case of Paired-end sequencing, a sequence will be read in the forward and reverse direction, for, let’s say, 150 cycles. This will result in the generation of reads of 150 bases. 💡 Here, the blue fragments, attached to the both ends of the insert, are called adapters. We can call them Read1 adapter and Read2 adapter. The Read1 adapter contains the Read1 primer-binding site, and the Read2 adapter contains the Read2 primer-binding site. These fragments will be read/ sequenced for a fixed number of cycles, in the forward and reverse direction. Sometimes there can be an inner distance that is left un-sequenced. Broadly speaking, this can happen when the read length (the number of read cycles) is less than the length of the DNA fragment/ insert. When there is shorter fragment, shorter than the length of the number of read cycles, the sequencing process can over-shoot to the adapter segment of the reverse read. This adapter contamination is called “adapter read-through”. This “adapter read-through” can result in the inclusion of either the full adapter or partial adapter being present in the read. How to remove adapters or other Illumina technical sequences like primers? It is very important to remove adapters, over represented sequences and other technical sequences like primers, from the reads. There are multitude of ways by which we can remove contaminating sequences. The popular bioinformatic pipelines that can potentially remove adapters are Trimmomatic, Trim_galore, fastp, etc. Removing adapters using Trimmomatic Adapters and other technical sequences can be removed using the Trimmomatic. The basic syntax of the Trimmomatic command is as follows. java -jar PE [-threads ] >] [-basein | ] [-baseout | The actual Trimmomatic program is a .jar file, and this needs JAVA to be installed in the machine. PE flag specifies that the read files are from a paired-end sequencing. -threads specifies the number of threads in a multi-core processor. -phred33 specifies the quality encoding of the base in the fastq nucleotide sequences. According to phred33 quality encoding, the quality score of each nucleotide base in the read is represented using ASCII encoding system. Normally, for each nucleotide base, the phred score is calculated as follows: phred = -10*log(probability that a base is wrong). This phred score is converted into an ASCII character via adding 33 to each phredscore value. -trimlog specifies a log file in text format. This log file contains information about the trimmed reads. The log file contains the following information: The read name The surviving sequence length. The amount of sequence trimmed from the start. The amount trimmed from the end. Next we need to specify the two input files. The first one is the forward fastq file, and the second one is the reverse fastq file. You may explicitly specify the fastq files. Otherwise, you may enter the name of the forward read followed by -basein flag, and the reverse read will be detected automatically. For example, if the name of the forward read file is Sample_Name_R1_001.fq.gz, the reverse file Sample_Name_R2_001.fq.gz will be detected automatically. If the forward file is Sample_Name.f.fastq, the program will look for Sample_Name.r.fastq for reverse file. Next, we need to specify the output files. There are four output files we need to specify. For paired forward read For unpaired forward read For paired reverse reads For unpaired reverse reads. The paired forward reads contain all the surviving forward reads. The unpaired forward reads contain all the reads whose partner reads were dropped. One may specify all the output files. Otherwise, a base-name can be provided following the --baseout flag, and the names of all output files will be created. For example, one may provide a base-name, mySampleFiltered.fq.gz, and the program will generate the output files, following a particular nomenclature pattern. o mySampleFiltered_1P.fq.gz - for paired forward reads o mySampleFiltered_1U.fq.gz - for unpaired forward reads o mySampleFiltered_2P.fq.gz - for paired reverse reads o mySampleFiltered_2U.fq.gz - for unpaired reverse reads All these parameters are called basic parameters, that are associated with setting up the work environment. Now we need to deal with the trimming modules. The first trimming module is called ILLUMINACLIP. This module accepts various parameters that are associated with finding and trimming the adapters. ILLUMINACLIP:::: Here, the parameters are: fastaWithAdaptersEtc: Here we need to specify the path to the fasta file containing the adapter sequence, primer sequence etc. seedMismatches: This is for a simple match strategy in Trimmomatic. Here, we need to specify the number of maximum mismatch count that would be taken into consideration as a full match. A higher value would result in the generation of more false positive data. palindromeClipThreshold: Here we need to specify how accurate the match between the two adapter-ligated reads must be for PE palindrome read alignment. simpleClipThreshold: Here, we need to specify how accurate the match between the adapter and the read must be. ILLUMINACLIP takes two additional parameters also. ILLUMINACLIP:::::: minAdapterLength: Here, we need to specify the length of the minimum length of the adapter, in which case the adapter will be removed. This is in addition to the palindromeClipThreshold value. The default value is 8. This means that as long as the minimum length of the adapter is 8 or greater, the adapter read-through will be trimmed. This helps to remove shorter adapter fragments. keepBothReads: This parameter applies to the palindrome match strategy. When an adapter read-through is found and the adapter removed, then the reverse read is the same as the forward one, except the fact that the latter is a reverse compliment of the former. In such a scenario, the reverse read will be dropped by the system. If keepBothReads is set to True, the reverse read will also be retained. The second module associated with Trimming is called SLIDINGWINDOW Here, a moving window of multiple-bases will be assessed for the average base quality. If the average base quality falls below a certain threshold, the sequence will be removed. So, even if the single base has a low quality, and the average base quality is still above the threshold, the sequence won’t be trimmed. The basic syntax of the SLIDINGWINDOW is: SLIDINGWINDOW:: Here, the is the number of bases to average across. is the threshold quality value The third module associated with Trimming: LEADING This module trims the base with a quality value below a certain threshold. The basic syntax is: LEADING: where the is the the value below which a base will be trimmed. The fourth important module associated with Trimming is, TRAILING This module takes a single argument of threshold quality value. If the last base in the read has a quality value below the set threshold, it will be removed. The basic syntax is: TRAILING: The fifth module associated with trimming the read: CROP With this module in action, bases, regardless of quality, can be trimmed from the end of the read, so that the read, maximally, has a length according to the set-value. CROP: where is the final length of the read after bases have been removed. The sixth module associated with Trimming the read: HEADCROP With this module in action, bases, regardless of their quality, will be trimmed from the beginning of the read, so that the read, maximally, has a length according to the set-value. HEADCROP: where is the number of bases to be removed from the beginning of the read The seventh module is associated with filtering the reads: MINLEN MINLEN takes a single argument. This specifies the minimum length of the read to be able to survive. After trimming, if the length of the read is less than the value stipulated in the module, the read will be dropped. TOPHRED33 This module is associated with the quality of the base. If the quality of the bases are calculated according to the ASCII encoding, TOPHRED33 is a valid module. More about the ILLUMINACLIP module ILLUMINACLIP does the finding and the clipping of the adapter. Technically, the adapter and other sequences like primer can appear in any location within the read, the most common mechanism of adapter contamination is when the sequencer reads a DNA fragment which is less than the length of the read. In this scenario, the read starts with a valid DNA sequence, and since the DNA fragment is shorter, the sequencer continues reading and “reads through” to the adapter of the opposite (reverse) read. This way, an adapter creeps in to the 3’-end of the read. The 5’-end starts with valid DNA sequence. This results in a partial or full adapter sequence in the 3’-end of the read. If the adapter sequence is full, it is easier to identify the adapter and remove it from the read. If it is a partial adapter, it is quite difficult to remove it. 💡 This is how adapters are ligated to the DNA fragments. 💡 How short inserts affect sequencing performance - Illumina Knowledge This picture shows very clearly how adapters are added to the DNA insert. The adapter complex is comprised of P5/P7 sequences, required to bind to the Flow Cell, Sample Indexes, Read1 and Read 2 primers. Now, these adapters are attached to the DNA insert. DNA inserts have Adenine over-hangs, and the adapters have Thymine over-hangs. As you see in the picture, Adenine and Thymine hybridize and the phosphate group binds with the 3’-OH, to form the backbone. Thus we have the sequencing reads. Now, the DNA inserts bind to the flow cell for sequencing. Primer binds to the primer-binding site on the read. How does Trimmomatic trim the adapter sequences In Single End sequencing method, where the DNA insert is sequenced only in one direction (there is no reverse pair), Trimmomatic uses a Simple strategy. Since the adapter sequence is known to us, continuous sections of each adapter, of 16bp length, will be matched with the read, in every possible positions. If there is a perfect match, or a close-match determined by the argument fed to the parameter, the entire alignment between the adapter and the read will be scored. If was fed with an argument of 2 (the default value), even if there is a mismatch in two bases between the “seed” and the read, it will be considered as a perfect match. Now, clipping of the adapter will happen based on the threshold value we feed, as an argument, to the parameter. For every match, a score will be calculated as follows: score = log10(probability that there is a random match of nucleotide base) score = log10(4) ~ 0.602 This, for each match of nucleotide base, the score increases by 0.602. If there is a mis-match, the score reduces by Q/10, where Q is the base-quality, that varies from 0 to 40. Therefore, the penalty value ranges from 0 to 4. The adapter will be trimmed as per the threshold score value we assign to the . If we provide a value of “7”, this is associated with a perfect match of 12 bases. If we want the adapter to be trimmed if there is a perfect match of 25 sequences, the threshold value is “15” This scenario can be understood well below. 💡 Trimmomatic: a flexible trimmer for Illumina sequence data | Bioinformatics | Oxford Academic (oup.com) The advantage of the SimpleClipThreshold method is that, it can detect any technical sequence located anywhere within the read. But as you can see in A and D, when there is only a short partial read-through/ match, detecting the adapter reliably is not possible. When there is a partial read-through at the 3’-end, Palindrome-Match is employed In palindrome mode, the adapter-ligated forward and the reverse reads are aligned, and the algorithm will calculate the score. If there is a matching score of 30, corresponding to a match of 50 bases, the adapters will be trimmed. 💡 Trimmomatic: a flexible trimmer for Illumina sequence data | Bioinformatics | Oxford Academic (oup.com) In the palindrome mode, the adapter-ligated forward and the reverse reads are aligned, and the global score will be calculated. The palindrome mode starts with an alignment of adapters and the start of each read, as you can see in A. If adapters directly joined, with no useful sequence information between them, such reads will be dropped. Then the testing process proceeds by moving the relative position of the reads backwards. When there is a match, both reads will be clipped. Even if the adapter is partial and reaches to the ligated-adapter, the reads will be clipped. The clipping will finish when there is no more over-hanging part that reaches to the adapter.

Demo files for Trimmomatic were downloaded from “bioinfogp.cnb.csic.es/files/samples/rnaseq/” The files were downloaded by using the following command: wget https://bioinfogp.cnb.csic.es/files/samples/rnaseq/RNA-Seq_Sample_Files.zip These are paired end-reads. For example, 1M_SRR9336468_1.fastq.gz and 1M_SRR9336468_2.fastq.gzarepaired end reads. Checking the number of reads in the fastq file We need to check the number of sequences in the fastq file. We have 1M_SRR9336468_1.fastq.gz and 1M_SRR9336468_2.fastq.gz. This can be achieved by the following command: zcat your_file.fastq.gz | wc -l | awk '{print $1/4}’ It turned out that there are 1000000 reads in the fastq file. For the purpose of testing, we need a simple fastq file, with a very few sequences. Let us cut down the fastq file to have, let’s say 10 reads in the file. Splitting the fastq file In order to cut down the fastq file, we will use seqkit split2 module. The following command can be used to achieve this task. seqkit split2 -1 1M_SRR9336468_1.fastq.gz -2 1M_SRR9336468_2.fastq.gz -s 10 -O out -f -e .gz This will generate 1 lakh parts of both forward and the reverse reads, containing 10 reads in each. For the sake of this demo-experimentation, we will have only very few fastq files. Now, let’s add forward and reverse adapter sequences to the respective read files. Appending the read files with adapter sequences We will use biopython for the purpose of appending the fastq reads with adapter sequences. Here is the code that we have crafted. from Bio.SeqIO.QualityIO import FastqGeneralIterator """Now we need to define the forward adapter and the reverse adapter""" forAdapter = "TACACTCTTTCCCTACACGACGCTCTTCCGATCT" revAdapter = "GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCT" """forward file path""" forPath = r"D:\website\wixSite\articles\trimmomaticDemo\inFiles\1M_SRR9336468_1.part_001.fastq." """Open the forward read file""" forHandle = open(forPath, "r", encoding="utf-8") """reverse file path""" revPath = r"D:\website\wixSite\articles\trimmomaticDemo\inFiles\1M_SRR9336468_2.part_001.fastq." """open the reverse file """ revHandle = open(revPath, "r", encoding="utf-8") #declare enmpty list forward newForFastq=[] #declare enmpty list reverse newRevFastq=[] """Iterate through each fastq block; declare itervariables to store header, sequence, and the quality""" for header, seq, qual in FastqGeneralIterator(forHandle): seqNewFor = seq.rstrip()+forAdapter newForFastq.append("@%s\n%s\n+\n%s\n" % (header, seqNewFor, qual)) """Iterate through each fastq block; declare itervariables to store header, sequence, and the quality""" for header, seq, qual in FastqGeneralIterator(revHandle): seqNewRev = seq.rstrip()+revAdapter newRevFastq.append("@%s\n%s\n+\n%s\n" % (header, seqNewRev, qual)) forPathOut = r"D:\website\wixSite\articles\trimmomaticDemo\outFiles\1M_SRR9336468_1.part_001.fastq" revPathOut = r"D:\website\wixSite\articles\trimmomaticDemo\outFiles\1M_SRR9336468_2.part_001.fastq" with open(forPathOut, "w")as forOut: for line in newForFastq: forOut.write(line) with open(revPathOut, "w")as revOut: for line in newRevFastq: revOut.write(line) We have added “TACACTCTTTCCCTACACGACGCTCTTCCGATCT” at the end of every sequence in the forward fastq file. And we have added GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCT at the end of every sequence in the reverse fastq file. Let us look at the first 10 lines of both forward and the reverse file. Forward ┌──(vijithkumar㉿LAPTOP-2MN1I195)-[~/nextGenSeq/RNA-Seq_Sample_Files/out] └─$ zcat 1M_SRR9336468_1.part_001.fastq.gz | head -n 10 @SRR9336468.1 1/1 NTGCATTTGAGAGGTTTGTGCACCTTGGAAAGTAGATTCATCGAATGGGTCGCCCACCTTGATGGATTCAGAAGCGGCTTTGAACTCTTCAATGAATTTGTCGTAAATAGATTCTTCAACATACACCCTCGAACCCGCACAACAGACCTC + #AAAFFAJJJJJJAAJJJ7FJJJJFJJJJJJJJJJJJ<FJJJJJJJJJJFJAJJJFJJ7<FJJJJJJJJFJJJFJJJJF7FFJJJJJJJJFFFAJJFFFFFFF-FJFAFFJ-7J7AFJFJJJFJFJJJJ--7FFAJ-F<JJJJJAFJ7-F @SRR9336468.2 2/1 NGTATCTGTTGTACTTTGGAATGTAATGCAAGTAAGCCCTTCTGATGACAATGGTACGGTGCATCTTGGTGGAGACGACGGTACCGGTCAAGATCTTACCACGGATGGAAACTAAACCAGTGAATGGACATTTCTTGTCAATGTAGGTAC + #AAFFJFJJJJJJJJJJJJFJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJFJFJJFFJFJFAJJJJJJJJFFJJJJFJJJJJJJJJJJJAJAJJJJJJJJJJJJFAJJJ7FFJJJF<AFAJFFFJJJJFFFJ-AF @SRR9336468.3 3/1 NGATGACATGTGGAAATATAACAGAGACAGAATATTTAATTTTATGTATGTAAGCATCGATTGGACACCAGGCTTATTGATGACCTTACTCGTCCAATTTGGCACGGACCGCTTTAACTTGCAAGTAGTTTTGTAAAGCATCAACAGACA Reverse ┌──(vijithkumar㉿LAPTOP-2MN1I195)-[~/nextGenSeq/RNA-Seq_Sample_Files/out] └─$ zcat 1M_SRR9336468_2.part_001.fastq.gz | head -n 10 @SRR9336468.1 1/2 NTTGGTATCTACTACAATTCTGGTGAGGTCTGTTGTGCGGGTTCAAGGGTGTATGTTGAAGAATCTATTTACGACAAATTCATTGAAGAGTTCAAAGCCGCTTCTGAATCCATCAAGGTGGGCGACCCATTCGATGAATCTACTTTCCAA + #AAAFJJJJJJJJJJJJJJJJJJJJJJJJFFJJJJJJJJJJJJJJJJJJJJJJJJJFJJFJJJJJJJJJJJJJJJJJJJFJJJ<AFJ-AF7FJJJJJJJJJ<AFJFJJ-77FFJJJAJFFFFFFFFFFJ7-<<JA7F<-AAFAAFJ-<FJ @SRR9336468.2 2/2 NTCCAAGAGAACCAAGAGATGGTACAAGAATGCCGGTTTGGGATTCAAGACCCCAAAAACCGCTATTGAAGGTTCCTACATTGACAAGAAATGTCCATTCACTGGTTTAGTTTCCATCCGTGGTAAGATCTTGACCGGTACCGTCGTCTC + #AAFFJFJFFJJJJJFJJJJJJJJJJJJJJJFJJJJFJJJJJJJFJJJFFJJJJJFJJJJJJJAJAJJJJJJFJJJJJJJJJJJJJJ-FJFFJJJJJJJFFJFJJJFFJJFFAJ<F7FFFJJ7JJF<JFJFAAAAAA-FFJ)--<<AA<F @SRR9336468.3 3/2 NACACCTCTAATATTAATACCGCCTTAAAAGTGGCTGATAGAGTTAATGCGGGTACGGTCTGGATAAACACTTATAACGATTTCCACCACGCAGTTCCTTTCGGTGGGTTCAATGCATCTGGTTTGGGCAGGGAAATGTCTGTTGATGCT Here, we can see that the length of the forward and the reverse read is 150bp. The forward and the reverse reads are not reverse complimentary to each other. Normally, when there is an adapter read-through, the forward and the reverse reads are complimentary to each other. Let’s try working with a single fastq file, for SE sequencing. Let’s run Trimmomatic. java -jar trimmomatic-0.39.jar SE -trimlog ~/nextGenSeq/TrimmomaticOutput/singleEnd1/Trimlog/trimlogFile.txt \ /mnt/d/website/wixSite/articles/trimmomaticDemo/outFiles/1M_SRR9336468_1.part_001.fastq \ ~/nextGenSeq/TrimmomaticOutput/singleEnd1/lane1_forward.fq.gz \ ILLUMINACLIP:/mnt/d/website/wixSite/articles/trimmomaticDemo/outFiles/TruSeq3.fa:2:30:10 \ MINLEN:36 So, let us look at the progress of the Trimmomatic. ┌──(vijithkumar㉿LAPTOP-2MN1I195)-[~/nextGenSeq/Trimmomatic/Trimmomatic-0.39] └─$ java -jar trimmomatic-0.39.jar SE -trimlog ~/nextGenSeq/TrimmomaticOutput/singleEnd1/Trimlog/trimlogFile.txt \ /mnt/d/website/wixSite/articles/trimmomaticDemo/outFiles/1M_SRR9336468_1.part_001.fastq \ ~/nextGenSeq/TrimmomaticOutput/singleEnd1/lane1_forward.fq.gz \ ILLUMINACLIP:/mnt/d/website/wixSite/articles/trimmomaticDemo/outFiles/TruSeq3.fa:2:30:10 \ MINLEN:36 TrimmomaticSE: Started with arguments: -trimlog /home/vijithkumar/nextGenSeq/TrimmomaticOutput/singleEnd1/Trimlog/trimlogFile.txt /mnt/d/website/wixSite/articles/trimmomaticDemo/outFiles/1M_SRR9336468_1.part_001.fastq /home/vijithkumar/nextGenSeq/TrimmomaticOutput/singleEnd1/lane1_forward.fq.gz ILLUMINACLIP:/mnt/d/website/wixSite/articles/trimmomaticDemo/outFiles/TruSeq3.fa:2:30:10 MINLEN:36 Automatically using 4 threads Using Long Clipping Sequence: 'AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC' ILLUMINACLIP: Using 0 prefix pairs, 1 forward/reverse sequences, 0 forward only sequences, 0 reverse only sequences Quality encoding detected as phred33 Input Reads: 10 Surviving: 10 (100.00%) Dropped: 0 (0.00%) TrimmomaticSE: Completed successfully Here, the report says that all the reads are surviving. No reads have been dropped. trimlog.txt; The log file of all the reads in the fastq file Let us look at the log report. SRR9336468.1 1/1 150 0 150 34 SRR9336468.2 2/1 150 0 150 34 SRR9336468.3 3/1 150 0 150 34 SRR9336468.4 4/1 150 0 150 34 SRR9336468.5 5/1 150 0 150 34 SRR9336468.6 6/1 150 0 150 34 SRR9336468.7 7/1 150 0 150 34 SRR9336468.8 8/1 150 0 150 34 SRR9336468.9 9/1 150 0 150 34 SRR9336468.10 10/1 150 0 150 34 Here, the columns are as follows. the read name the surviving sequence length the location of the first surviving base, aka. the amount trimmed from the start the location of the last surviving base in the original read the amount trimmed from the end The first column contains list of all the reads, prefixed by SRR. The second column contains the length of the surviving reads. The third column contains the position of the fist surviving base. It is zero for all the reads. The fourth column contains the position of the last surviving base. It is 150 for all the reads. The last column contains the amount trimmed from the end of the read. It is 34. Since MINLEN was set to 36 (means, drop the reads if the length of the read falls below 36), no reads were dropped. Setting MINLEN to 151 ┌──(vijithkumar㉿LAPTOP-2MN1I195)-[~/nextGenSeq/Trimmomatic/Trimmomatic-0.39] └─$ java -jar trimmomatic-0.39.jar SE -trimlog ~/nextGenSeq/TrimmomaticOutput/singleEnd1/Trimlog/trimlogFile.txt /mnt/d/website/wixSite/articles/trimmomaticDemo/outFiles/1M_SRR9336468_1.part_001.fastq ~/nextGenSeq/TrimmomaticOutput/singleEnd1/lane1_forward.fq.gz ILLUMINACLIP:/mnt/d/website/wixSite/articles/trimmomaticDemo/outFiles/TruSeq3.fa:2:30:10 MINLEN:151 TrimmomaticSE: Started with arguments: -trimlog /home/vijithkumar/nextGenSeq/TrimmomaticOutput/singleEnd1/Trimlog/trimlogFile.txt /mnt/d/website/wixSite/articles/trimmomaticDemo/outFiles/1M_SRR9336468_1.part_001.fastq /home/vijithkumar/nextGenSeq/TrimmomaticOutput/singleEnd1/lane1_forward.fq.gz ILLUMINACLIP:/mnt/d/website/wixSite/articles/trimmomaticDemo/outFiles/TruSeq3.fa:2:30:10 MINLEN:151 Automatically using 4 threads Using Long Clipping Sequence: 'AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC' ILLUMINACLIP: Using 0 prefix pairs, 1 forward/reverse sequences, 0 forward only sequences, 0 reverse only sequences Quality encoding detected as phred33 Input Reads: 10 Surviving: 0 (0.00%) Dropped: 10 (100.00%) TrimmomaticSE: Completed successfully Here, we can see that all the 10 reads have been dropped. SRR9336468.1 1/1 0 0 0 0 SRR9336468.2 2/1 0 0 0 0 SRR9336468.3 3/1 0 0 0 0 SRR9336468.4 4/1 0 0 0 0 SRR9336468.5 5/1 0 0 0 0 SRR9336468.6 6/1 0 0 0 0 SRR9336468.7 7/1 0 0 0 0 SRR9336468.8 8/1 0 0 0 0 SRR9336468.9 9/1 0 0 0 0 SRR9336468.10 10/1 0 0 0 0 Changing the simpleClipThresholdvalue to 21 simpleClipThreshold is a parameter of the ILLUMINACLIP module. So, once a “seed” of 16bp are aligned in every position possible with the read, the program will look for a match between the seed and the read. If a match is established, the full adapter will be aligned with the read and a score will be calculated. Score = log10(probability that there is random match) ~0.602. So, for 34 bases, the score is less than 21 Look at the log report: SRR9336468.1 1/1 184 0 184 0 SRR9336468.2 2/1 184 0 184 0 SRR9336468.3 3/1 184 0 184 0 SRR9336468.4 4/1 184 0 184 0 SRR9336468.5 5/1 184 0 184 0 SRR9336468.6 6/1 184 0 184 0 SRR9336468.7 7/1 184 0 184 0 SRR9336468.8 8/1 184 0 184 0 SRR9336468.9 9/1 184 0 184 0 SRR9336468.10 10/1 184 0 184 0 Here, we can see that no trimming has been done. Changing the simpleClipThresholdvalue to 20 Here, we have 34 sequences, and therefore the maximum permissible score is: 34 x log10(4) =0.602 x 34 =20.468 ~20 Let us look at the log report: SRR9336468.1 1/1 150 0 150 34 SRR9336468.2 2/1 150 0 150 34 SRR9336468.3 3/1 150 0 150 34 SRR9336468.4 4/1 150 0 150 34 SRR9336468.5 5/1 150 0 150 34 SRR9336468.6 6/1 150 0 150 34 SRR9336468.7 7/1 150 0 150 34 SRR9336468.8 8/1 150 0 150 34 SRR9336468.9 9/1 150 0 150 34 SRR9336468.10 10/1 150 0 150 34 See, here we can see that the trimming has been done! What it means is that, if you set 20 as the argument forsimpleClipThreshold , even if there is a couple of mismatches, the corresponding reads won’t be clipped. Changing two bases in the adapter sequence in the read and setting the simpleClipThreshold to 20 @SRR9336468.1 1/1 NTGCATTTGAGAGGTTTGTGCACCTTGGAAAGTAGATTCATCGAATGGGTCGCCCACCTTGATGGATTCAGAAGCGGCTTTGAACTCTTCAATGAATTTGTCGTAAATAGATTCTTCAACATACACCCTCGAACCCGCACAACAGACCTCAGATCGGAAGAGCACACGTCTGAACTCCAGTCGG + #AAAFFAJJJJJJAAJJJ7FJJJJFJJJJJJJJJJJJ<FJJJJJJJJJJFJAJJJFJJ7<FJJJJJJJJFJJJFJJJJF7FFJJJJJJJJFFFAJJFFFFFFF-FJFAFFJ-7J7AFJFJJJFJFJJJJ--7FFAJ-F<JJJJJAFJ7-FJJJJJJJJJJJJJJFFJJJJJJJJJJJJJJJJJJ @SRR9336468.2 2/1 NGTATCTGTTGTACTTTGGAATGTAATGCAAGTAAGCCCTTCTGATGACAATGGTACGGTGCATCTTGGTGGAGACGACGGTACCGGTCAAGATCTTACCACGGATGGAAACTAAACCAGTGAATGGACATTTCTTGTCAATGTAGGTACAGATCGGAAGAGCACACGTCTGAACTCCAGTCAC Look at the above fastq sequence. I have purposefully changed two bases in the adapter. This shows that only 32 out of 34 bases score and there are two penalties. So, let’s calculate the score. score = 0.602 x log10(4) =0.602 x 32 =19.264 ~19 So, when we set the simpleClipThreshold value to 20, the read whose adapter has got changes in two bases fail to meet the threshold that was set, therefore the adapter of the said read won’t get trimmed. Let’s look at the log report here. SRR9336468.1 1/1 184 0 184 0 SRR9336468.2 2/1 150 0 150 34 SRR9336468.3 3/1 150 0 150 34 SRR9336468.4 4/1 150 0 150 34 SRR9336468.5 5/1 150 0 150 34 SRR9336468.6 6/1 150 0 150 34 SRR9336468.7 7/1 150 0 150 34 SRR9336468.8 8/1 150 0 150 34 SRR9336468.9 9/1 150 0 150 34 SRR9336468.10 10/1 150 0 150 34 See here, the adapter hasn’t been clipped from the first read. When the read-through adapters are short and partial @SRR9336468.1 1/1 NTGCATTTGAGAGGTTTGTGCACCTTGGAAAGTAGATTCATCGAATGGGTCGCCCACCTTGATGGATTCAGAAGCGGCTTTGAACTCTTCAATGAATTTGTCGTAAATAGATTCTTCAACATACACCCTCGAACCCGCACAACAGACCTCAGATCGGAAGAGCACACGTCTGAACTCCA + #AAAFFAJJJJJJAAJJJ7FJJJJFJJJJJJJJJJJJ<FJJJJJJJJJJFJAJJJFJJ7<FJJJJJJJJFJJJFJJJJF7FFJJJJJJJJFFFAJJFFFFFFF-FJFAFFJ-7J7AFJFJJJFJFJJJJ--7FFAJ-F<JJJJJAFJ7-FJJJJJJJJJJJJJJFFJJJJJJJJJJJJJ @SRR9336468.2 2/1 NGTATCTGTTGTACTTTGGAATGTAATGCAAGTAAGCCCTTCTGATGACAATGGTACGGTGCATCTTGGTGGAGACGACGGTACCGGTCAAGATCTTACCACGGATGGAAACTAAACCAGTGAATGGACATTTCTTGTCAATGTAGGTACAGATCGGAAGAGCACACGTCTGAA + #AAFFJFJJJJJJJJJJJJFJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJFJFJJFFJFJFAJJJJJJJJFFJJJJFJJJJJJJJJJJJAJAJJJJJJJJJJJJFAJJJ7FFJJJF<AFAJFFFJJJJFFFJ-AFJJJJJJJJJJJJJJFFJJJJJJJJ @SRR9336468.3 3/1 NGATGACATGTGGAAATATAACAGAGACAGAATATTTAATTTTATGTATGTAAGCATCGATTGGACACCAGGCTTATTGATGACCTTACTCGTCCAATTTGGCACGGACCGCTTTAACTTGCAAGTAGTTTTGTAAAGCATCAACAGACAAGATCGGAAGAGCACACGT + #AAFFJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ<JJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJAJJFAJJJJFJJJJJJJJJJJFJJJJFJJJFJFJJFAJJJJFFJJJJF<7FJJ-7AFAFFFJAFJJFFJJJFFFFJFJJJJJJJJJJJJJJJJFFJJJ @SRR9336468.4 4/1 NATAGTTATAAAAATATATGTATATAAAGTATGTGAGTGTGTGAGTGTGTGAAAAGGATGTTATTAAATAGTATGTTATTATTTCTTGTCGTATCTGGAGTAGTATTTCTTTGAAATTATACCGTCCAGCGTCAAAAGACCCAGAGCGCCAGATCGGAAGAGCA + #AAFFJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJFJJJJJJJJJJJJJJJJ7JJJJJJJJJJJJJJJJJJJFJJJJJJJJJJJJJJJJJJJFJJJJFJJJJJJJFFJJJ-FJJJJJF-AFJ<7<<JJJJJ<JJJJJFJJJFJJJJJJJJJJJJJJJJ @SRR9336468.5 5/1 NGTTGAGTAACTTGACCTTCACACGTAGTTGCAAGTACCGATGCAAAATTTTTACACTGGTGACAACCTTAATGGTATAATCCAATTAGCGACATCACCGAGAAAGTATCCAAGTCGATTTGATGAATAAAATAATAACTTAAACATGCAAGATCGGAAGAGCACACGTCTGAACTCCAGTCAC Here, 5, 10, 15, 20 bases have been clipped from the first, second, third and the fourth fastq files. The simpleClipThreshold is set to 20. This quality score is associated with a perfect match of 34 bases. The following is the log report. SRR9336468.1 1/1 179 0 179 0 SRR9336468.2 2/1 174 0 174 0 SRR9336468.3 3/1 169 0 169 0 SRR9336468.4 4/1 164 0 164 0 SRR9336468.5 5/1 150 0 150 34 Adapters were not clipped from the first four reads. The simpleClipThreshold value was set to 15. This is associated with a perfect match of 25 bases. The log report is: SRR9336468.1 1/1 150 0 150 29 SRR9336468.2 2/1 174 0 174 0 SRR9336468.3 3/1 169 0 169 0 SRR9336468.4 4/1 164 0 164 0 SRR9336468.5 5/1 150 0 150 34 read 1 has 29 bases in the adapter, and it meets a quality score of 15. @SRR9336468.1 1/1 NTGCATTTGAGAGGTTTGTGCACCTTGGAAAGTAGATTCATCGAATGGGTCGCCCACCTTGATGGATTCAGAAGCGGCTTTGAACTCTTCAATGAATTTGTCGTAAATAGATTCTTCAACATACACCCTCGAACCCGCACAACAGACCTC + #AAAFFAJJJJJJAAJJJ7FJJJJFJJJJJJJJJJJJ<FJJJJJJJJJJFJAJJJFJJ7<FJJJJJJJJFJJJFJJJJF7FFJJJJJJJJFFFAJJFFFFFFF-FJFAFFJ-7J7AFJFJJJFJFJJJJ--7FFAJ-F<JJJJJAFJ7-F The simpleClipThreshold value was set to 10. This is associated with a perfect match of 16 bases. Since the first read has 29 bases, the match score will be higher than 10, so the adapter will be clipped. Since the second read has 24 bases, the match score will be still above 10, and therefore will be clipped. The third read has 19 bases, the match score is still above 10, and therefore the adapter will be clipped. The simpleClipThreshold value was set to 1. This is associated with a perfect match of 2 bases. But in such a scenario, there are chances of false positives being generated. Let’s look at the trim log report: SRR9336468.1 1/1 150 0 150 29 SRR9336468.2 2/1 150 0 150 24 SRR9336468.3 3/1 150 0 150 19 SRR9336468.4 4/1 148 0 148 16 SRR9336468.5 5/1 150 0 150 34

This post is about pyYAML library of Python to work with yaml files.