Atlassian uses cookies to improve your browsing experience, perform analytics and research, and conduct advertising. Accept all cookies to indicate that you agree to our use of cookies on your device.
Atlassian uses cookies to improve your browsing experience, perform analytics and research, and conduct advertising. Accept all cookies to indicate that you agree to our use of cookies on your device. Atlassian cookies and tracking notice, (opens new window)
As mentioned in the introduction tutorial as well as the read processing tutorial, read processing can make a huge impact on downstream work. While cutadapt which was introduced in the read processing tutorial is great for quick evaluation or dealing with a single bad sample, it is not as robust as some other trimmers in particular when it comes to removing sequence that you know shouldn't be present but may exist in odd orientations (such as adapter sequences from the library preparation). This tutorial is adapted from the 2021 trimmomatic tutorial which sought to do the same basic things as fastp: get rid of adapter sequences first and foremost, ideally even before fastQC so you can make any quality or length based improvements on actual data not artifacts. The #1 biggest reason why fastp is now the instructor's preferred trimming program is this box taken from the trimmomatic tutorial:
A note on the adapter file used here
The adapter file listed here is likely the correct one to use for standard library preps that have been generated in the last few years, but may not be appropriate for all library preps (such as single end sequencing adapters, nextera based preps, and certainly not appropriate for PacBio generated data). Look to both the trimmomatic documentation and your experimental procedures at the bench to figure out if the adapter file is sufficient or if you need to create your own.
The more collaborative your work is, the less confidence you will have in picking the correct adapter file with trimmoatic, and while thanks to conda installations it can be pretty easy to test multiple different adapter filees, fastp does all the guess work for you, and can generate some interesting graphs itself.
Learning objectives:
Install fastp
Remove adapter sequences from some plasmids and evaluate effect on read quality, or assembly.
Installing fastp
fastp's home page can be found on github and has links to the paper discussing the program, installation instructions for conda, and information on each of the different options available to the program. This is far above the quality the average programs will have as most will not have a user manual (or not nearly so detailed), may not have been updated since originally published (or may not have been published), etc. It having been updated since the publication is one thing that makes it such a good tool as the more who use it the more likely problems are found, and having a group who is going to actively improve the program will significantly increase its longevity.
There actually are not a lot of "wrong" answers here at least from the theoretical side. As read processing takes place upstream of basically all other analysis steps it makes sense to put it in almost every environment. Practically, though that means that it should be installed in every environment which starts to defeat the purpose of having any different environments at all. As will be discussed on Friday, you might want to start thinking about grouping programs into chunks. Almost no matter what analysis you do, you are going to want to trim adapters (fastp), check the quality (fastqc) and likely compare to other similar samples (multiqc). So putting all these programs into a single "read pre-processing" environment seems like a good grouping.
At this point in the class you can start making your own calls about what environments you want to put programs in, or what names you want to give them. While you can keep using the same names and groupings I suggest, last year there was feedback that having to make the choices of how to modify commands based on different environments was helpful.
number of additional processors (48 max on stampede2)
--detect_adapter_for_pe
automatically detect adapter sequence based on paired end reads, and remove them
-j <LOG.json>
json file with information about how the trim was accomplished. can be helpful for looking at multiple samples similar to multiqc analysis
name of json file you want to use
-h <LOG.html>
html file with infomration similar to the json file, but with graphs
name of html file you want to use
All of the above has been put together from the help fastp --help command.
Trimming a single sample
Get some data
set up directories and copy files
mkdir -p $SCRATCH/GVA_fastp_1sample/Trim_Reads $SCRATCH/GVA_fastp_1sample/Raw_Reads
cd $SCRATCH/GVA_fastp_1sample
cp $BI/gva_course/plasmid_qc/E1-7* Raw_Reads
The ls command should show you 2 gzipped fastq files. You may notice that here that we used a wildcard in the middle of our copy path for the first time. This is done so that you can grab both R1 and R2 easily without having to type out the full command. Double tab will help tell you when you have a sufficiently specific base name to only get the files you are after.
According to the --help information: "-p, --parents no error if existing, make parent directories as needed" so it is allowing us to make nested directories rather than having to make them 1 at a time. Additionally we use a ::space:: to create 2 directories at the same time. Almost every command has more information about it that can be read at the command line We have used -h and --help and tried calling commands without any options and mentioned the 'man' command throughout the course for the various programs we have installed, but here we see we can actually use that same framework to access more information about even the most basic of commands without even needing the internet.
Trim the fastq files
The following command can be run on the head node. Like with FastQC if we are dealing with less than say 1-2Million reads, it is reasonable to run the command on the head node unless we have 100s of samples in which case submitting to the queue will be faster as the files can be trimmed all at once rather than 1 at a time. Use what you have learned in the class to determine if you think this command should be run on the head node. (this was covered in more detail in the first part of the evaluating and processing read quality tutorial.)
zgrep -c "^+$" Raw_Reads/*.fastq.gz
Example command for trimming illumina paired end adapters
Most likely cause here is that you forgot to activate your new conda environment if you have another issue, you will likely want to ask a question.
Evaluating the output
Using everything you have learned so far in the class, can you answer the following questions?
E1-7_S187_L001_R2_001.trim.fastq.gz and E1-7_S187_L001_R1_001.trim.fastq.gz These were created with the -o and -O options, they are in the Trim_Reads folder, and you likely found them using the ls command
fastp.html and fastp.json These are log files created by default since we didn't specify their names. This is part of why -j and -h were discussed above with the general command. While the json file can be evaluated in the terminal (cat less more head tail), the html file has to be transferred back to your computer to view.
5884 paired end reads 11768 total reads You likely found this out from using the zgrep command, or from the following blocks that printed as the command ran:
Read1 after filtering:
total reads: 5884
total bases: 791763
Q20 bases: 782948(98.8867%)
Q30 bases: 765510(96.6842%)
Read2 after filtering:
total reads: 5884
total bases: 791763
Q20 bases: 711414(89.8519%)
Q30 bases: 658164(83.1264%)
Filtering result:
reads passed filter: 11768
reads failed due to low quality: 2014
reads failed due to too many N: 0
reads failed due to too short: 0
reads with adapter trimmed: 3970
bases trimmed due to adapters: 193972
From the information generated while the command ran we see: Insert size peak (evaluated by paired-end reads): 171 This tells us that the average peak size was 171 bases, and that it was estimated by looking at the overlap between the read pairs. It is potentially inaccurate as reads which do not overlap each other can not estimate the size. If you transferred the .html file back to your laptop, you would see this relevant histogram:
The general section of the summary at the top of the html tells us that the average insert size was 171, while the histogram tells us that 50% of our data is <18 or >272 bases
If you only look at the information that printed to the screen, you probably answer "No" you likely see the following block and think this is the end of the answer: Detecting adapter sequence for read1...
No adapter detected for read1
Detecting adapter sequence for read2...
No adapter detected for read2 A more fuller answer might be "maybe" or "probably" or "I'm not sure" as: 1. Not finding any adapter would be super rare 2. If 45% of our reads have an insert size of 171 bases, and we did 151bp PE sequencing, we should be able to find adapter sequences 3. in the filtering results we see: Filtering result:
reads passed filter: 11768
reads failed due to low quality: 2014
reads failed due to too many N: 0
reads failed due to too short: 0
reads with adapter trimmed: 3970
bases trimmed due to adapters: 193972 If you look at the html file you probably answered "yes" There is a section for Read1 and Read2 adapters which show a growing stretch of DNA which recreates the illumina adapter sequences.
Like we saw in our fastqc reports (over represented sequences having "no hit" and adapter content staying at bottom of graph), for something to be classified as an "adapter" in the first section of the printed information, it has to meet certain criteria that in this (and many other instances) is perhaps a bit too stringent.
This is pretty open ended, take a look at the html file in particular, see what of it does or doesn't make sense and consider asking a question if you would like to know more. Of several things that you may stand out to you is large fraction of reads end with stretches of "G" on the end. There are 2 things to note with this: 1. 2 color sequencing on illumina (detailed information here) reads "no color" as "G", 2. This library is very fragmented and contains adapter dimers meaning that in some cases there are only ~40bp downstream of the sequencing primer location leaving 60 cycles that have no template available. If you look at the help for fastp the following options may stand out to you as a way to deal with this: -g, --trim_poly_g force polyG tail trimming, by default trimming is automatically enabled for Illumina NextSeq/NovaSeq data --poly_g_min_len the minimum length to detect polyG in the read tail. 10 by default. (int [=10]) -G, --disable_trim_poly_g disable polyG tail trimming, by default trimming is automatically enabled for Illumina NextSeq/NovaSeq data -x, --trim_poly_x enable polyX trimming in 3' ends. --poly_x_min_len the minimum length to detect polyX in the read tail. 10 by default. (int [=10]) Consider rerunning the fastp command while adding "-g" to the command line and see how the results differ.