Table of Contents |
---|
...
Often, the first thing you (or your boss) want to know about your sequencing run is simply, "how many reads do I haveare there?". For the $BI/gva_course/mapping/data/SRR030257_1.fastq file, the answer is 3,800,180. How can we figure that out?
The grep (or Global Regular Expression Print) command can be used to determine the number of lines which match some criteria as shown above. Above we searched used it to search for:
- anything from the group of ACTGN with the [] marking them as a group
- matching any number of times *
- from the beginning of the line ^
- to the end of the line $
Here, since we are only interested in the number of reads that we have, we can make use of knowing the 3rd line in the fastq file is a + and a + only, and grep's -c option to simply report the number of reads in a file.
Code Block | ||||||
---|---|---|---|---|---|---|
| ||||||
grep -c "^+$" $BI/gva_course/mapping/data/SRR030257_1.fastq |
...
Expand | ||
---|---|---|
| ||
The Per base sequence quality report does not look great. If I were making the call to trim based soley on this I'd probably pick 31 or 32 as the last base as this is the first base that the average quality score drops significantly. 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 Illumina runs 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
Cutadapt
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. Cutadapt provides a simple command line tool for manipulating fasta and fastq files. The program description on their website provides 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.
...
This is another example of different solutions giving the same final product, and how careful reading of documentation may improve your work. NGS 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 reproducible.
A note on versions
Info | |||||||||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
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? .
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.10 being released on April 22nd this year.
This won't be the last time we mention different program versions. |
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 in overall results.
...