Overview
Before you start the alignment and analysis processes, it can be useful to perform some initial quality checks on your raw data. If you don't do this (or even if you do), you may notice later that something looks fishy in the the output: for example, many of your reads are not mapping or the ends of many of your reads do not align. Both can give you clues about whether you need to process the reads to improve the quality of data that you are putting into your analysis.
...
| Expand | ||
|---|---|---|
| ||
From the OPTIONS: section of the idev help output: -m minutes sets time in minutes (default: 30) -r reservation_name requests use of a specific reservation -A account_name sets account name (default: -A none) So you requested an idev node for 180 minutes, using the reservation named CCBB_Day_1, and asked that it be charged to the account named UT-2015-05-18. |
...
| Expand | |||||||
|---|---|---|---|---|---|---|---|
| |||||||
grep or Global Regular Expression Print can also be used to determine the number of lines which match some criteria. Since we know the 3rd line in the fastq file is a + and a + only, we can look for a line that only has a + in it, and use that number to determine the number of sequence blocks in the file.
the -c option tells grep to count the lines (rather than printing them all to the screen and tell you how many it found. The characters between the "" is what grep is looking for. The ^ symbol means, look for the beginning of the line, the $ symbol means look for the end of the line. Once again you see this returns 3800180 reads. |
While checking the number of reads a file has can solve some of the most basic problems, it doesn't really provide any direct evidence as to the quality of the sequencing data. To get this type of information before starting meaningful analysis other programs must be used.
...
FastQC is available from the TACC module system on lonestar. Interactive GUI versions are also available for Windows and Macintosh and can be downloaded from the Babraham Bioinformatics web site. We don't want to clutter up our work space so copy the SRR030257_1.fastq file to a new directory named BDIBGVA_fastqc_tutorial on scratch, use the module system to load fastqc, use fastqc's help option after the module is loaded to figure out how to run the program. Once the program is completed use scp to copy the important file back to your local machine (The bold words are key words that may give you a hint of what steps to take next)
| Code Block | ||||
|---|---|---|---|---|
| ||||
mkdir $SCRATCH/BDIBGVA_fastqc_tutorial cd $SCRATCH/BDIBGVA_fastqc_tutorial cp $BI/gva_course/mapping/data/SRR030257_1.fastq . module load fastqc fastqc -h # examine program options fastqc SRR030257_1.fastq # run the program |
...
| Code Block | ||||
|---|---|---|---|---|
| ||||
# on tacc terminal pwd # on new terminal of local computer scp <username>@ls5.tacc.utexas.edu:<pwd_results_from_other_window>/SRR030257_1_fastqc.html ~/Desktop # open the newly transferedtransferred file from from the desktop and see how the data looks |
...
| Expand | ||
|---|---|---|
| ||
The Per base sequence quality report does not look great, but more importantly, nearly 1.5% of all the sequences are all A's according to the Overrepresented sequences. This is something that often comes up in miseq data that has shorter insert sizes than the overall read length. Next we'll start looking at how to trim our data before continuing. |
FASTQ Processing Tools
FASTX Toolkit
...
Cutadapt
Cutadapt provides a set of simple command line tools tool for manipulating fasta and fastq files. The available modules are described program description on their website . They include a fast fastx_trimmer utility for trimming fastq sequences (and quality score strings) before alignment and fastx_clipper for trimming specific sequences (such as adapters).FASTX-Toolkit is available via the TACC module systemprovides good details of all the capabilities and examples for some common tasks. Cutadapt is also available via the TACC module system allowing us to turn it on when we need to use it and not worry about it other times.
| Code Block | ||
|---|---|---|
| ||
module spider fastx_toolkitcutadapt module load fastx_toolkitcutadapt |
Trimming low quality bases
...
Low quality base reads from the sequencer can cause an otherwise mappable sequence not to align. There are a number of open source tools that can trim off 3' bases and produce a FASTQ file of the trimmed reads to use as input to the alignment program, but fastx_trimmer has the but cutadapt has the advantage of being a module on TACC and therefore the easiest to use. By default To run the program reads its input data from standard input and writes trimmed sequences to standard output, use what you know about piping and printing text to the screen to trip the fastq file to 34 bases.
...
| title | Hint |
|---|
, you simply type 'cutadapt' followed by whatever options you want, and then the name of the fastq files without any option in front of it. Use the -h option to see all the different things you can do and see if you can figure out how to trim the reads down to 34 bases.
| Expand | ||
|---|---|---|
| ||
Type cutadapt -h to see program documentation. Look As there is a large number of options, look below the possible solution for more detailed information of what to focus on |
| Code Block | ||||
|---|---|---|---|---|
| ||||
cat SRR030257_1.fastq | fastx_trimmer cutadapt -l 34 >-o SRR030257_1.trimmed.fastq |
...
SRR030257_1.fastq |
...
- The -l 34 option says that base 34 should be the last base (i.e., trim down to 34 bases)
- The > redirects that output to the new file on the right side -o sets the output file, in this case SRR030257_1.trimmed.fastq
- Listing the input file without any option in front of it (SRR030257_1.fastq) is a common way to specify input files.
Exercise: compressing the fastx_trimmer outputtrimmed file
Compressed files are smaller, easier to transfer, and many programs allow for their use directly. How would you tell fastx_trimmercutadapt to compress (gzip) its output file?
| Expand | ||
|---|---|---|
| ||
Type fastx_trimmer cutadapt -h to see program documentation | ||
and look for information about compressed files |
| Expand | ||||||
|---|---|---|---|---|---|---|
| cat SRR030257_1.fastq | fastx_trimmer -l 34 -z >
| |||||
Above the citation you see a paragraph that starts: Input may also be in FASTA format. Compressed input and output is supported and auto-detected from the file name (.gz, .xz, .bz2) So simply by adding .gz to the output file name, cutadapt will compress it after it does the trimming. |
| Code Block | ||||
|---|---|---|---|---|
| ||||
cutadapt -l 34 -o SRR030257_1.trimmed.fastq.gz SRR030257_1.fastq |
| Code Block | ||||
|---|---|---|---|---|
| ||||
catgzip SRR030257_1.fastq | fastx_trimmer -l 34 | gzip > SRR030257_1.trimmed.fastq.gz.trimmed.fastq |
Both of the above solutions give the same final product, but are clearly achieved in different ways. This is done to show you that data analysis is a results driven process, if the result is correct, and you know how you got the result it is correct as long as it is reproducablereproducible.
Adapter trimming
As mentioned above, fastx_clipper cutadapt can be used to trim specific sequences, and based on our fastqc analysis, the sequence AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA is significantly overrepresented in our data. How would you use fastx_clippercutadapt to remove those sequences from the fastq file?
| Expand | ||
|---|---|---|
| ||
Type fastx_clipper -h to see program documentationAgain, we go back to the program documentation to find what we are looking for: cutadapt -h Look below the possible solution for more detailed information on what to focus on . | ||
| Code Block | ||
if you cant find what you are looking for . |
| Code Block | ||||
|---|---|---|---|---|
| ||||
fastx_clippercutadapt -i SRR030257_1.trimmed.fastq -o SRR030257_1.trimmed.depleted.fastq -a AAAAAAAAAAAAAAAAAAAA -l 34 -n m 16 SRR030257_1.fastq |
| Command portion | purpose | -i SRR030257_1.trimmed.fastq | use this file as input|
|---|---|---|---|
| -o SRR030257_1.trimmed.depleted.fastq | create this new output file | ||
-a AAAAAAAAAAAAAAAAA | remove bases containing this sequence | ||
| -l 34 | trim reads to 34 bases | ||
| -m 16 | discard any read shorter than 34 16 bases after sequence removed | ||
| -n | keep reads containing "N" bases in them. Consider how they are treated in downstream applications |
Other fastx toolkit programs
What other fastx and fastq manipulation programs are part of the fastx toolkit? The available modules are described on their website, but we can also learn some things from the command prompt using module spider fastx_toolkit.
| Expand | |||||
|---|---|---|---|---|---|
| |||||
Type fast then tab twice to see their names
|
...
| as these are more likely difficult to uniquely align to the genome | |
| SRR030257_1.fastq | use this file as input |
From the summary printed to the screen you can see that this removed a little over an additional 2.2M bp of sequence.
A note on versions
In our first tutorial we mentioned how knowing what version of a program you are using can be. When we loaded the the cutadapt module we didn't specify what version to load. Can you figure out what version you used, and what the most recent version of the program there is? .
| Expand | |||||||||
|---|---|---|---|---|---|---|---|---|---|
| |||||||||
try using the module system or the program's help files
|
Figuring out the most recent version is a little more complicated. Unlike programs on your computer like Microsoft Office or your internet browser, there is nothing in an installed program that tells you if you have the newest version or even what the newest version is. If you go to the programs website (easily found with google or this link), the changes section lists all the versions that have been list with v2.3 being released on April 25th this year.
| Expand | ||
|---|---|---|
| ||
The biggest reason is that someone at tacc has to go through a process of noticing that there is a new version, figuring out if all of the changes are compatible with tacc, installing it, and then fielding questions and problems from users who were used to using the old version and have problems with the new version somehow. The next reason is that the existing version works, and if you read through some of the recent changes, are very small and do not effect the function of the program very much. Together, this is why I encourage you to make note of what version of the programs you use when you use them (primarily by loading modules complete with the versions in your .bashrc file), and to consider installing programs yourself when appropriate (as is discussed in the advanced trimmomatic tutorial for read trimming). |
Optional Exercise: Improve the quality of R2 the same way you did for R1.
Unfortunately we don't have time during the class to do this, but as a potential exercise in your free time, you could improve R2 the same way you did R1 and use the improved fastq files in the subsequent read mapping and variant calling tutorials to see the difference it can make.
Return to GVA2017 GVA2019 course pageGlobal regular expression print