You are viewing an old version of this page. View the current version.
Compare with Current
View Page History
« Previous
Version 15
Next »
Quality check the reads
Navigate to exercise location
Investigate reads. We paid the sequencing facility for 50,000 paired end reads.
How can be confirm we got this?
#we have two fastq files
ls -lh *.fastq
#take a look at the top of the files
less Lib01_R1.fastq
#how many lines are in each fastq file?
wc -l Lib01*.fastq
#look at how many reads there are (note that fastq generally files have 4 lines per read)
expr $(cat Lib01_R1.fastq | wc -l) / 4
expr $(cat Lib01_R2.fastq | wc -l) / 4
Paired end read files should have matching names read for read.
Double-check this is case for your reads
#this can be checked many ways, here is one option:
#search for lines that start with @ symbol (the read definition lines in our fastqs)
#and save the top 10 of them as a text file.
#Repeat for the R2 reads.
#Then line them up
grep "^@" Lib01_R1.fastq | head > first_10_R1_read_names.txt
grep "^@" Lib01_R2.fastq | head > first_10_R2_read_names.txt
paste first_10_R1_read_names.txt first_10_R2_read_names.txt
What if you wanted to repeat the above solution to check the bottom of the files as well?
grep "^@" Lib01_R1.fastq | tail > last_10_R1_read_names.txt
grep "^@" Lib01_R2.fastq | tail > last_10_R2_read_names.txt
paste last_10_R1_read_names.txt last_10_R2_read_names.txt
Based on the library preparation, we expect that our forward reads were cut with the restriction enzyme nlaIII (NLA3).
We expect that the vast majority of our reads should have the recognition sequence in them.
How can we check that our reads largely fit this expectation?
#Looking up nlaIII, we see it's cut site:
# CATG'
# 'GTAC
#check if we see this cut site in the forward reads
grep CATG Lib01_R1.fastq | wc -l
#compare with total read count
expr $(cat Lib01_R1.fastq | wc -l) / 4
Next run the program FastQC on the reads
module load fastqc
mkdir raw_Fastqc_Restults
fastqc -o raw_Fastqc_Restults/ -f fastq Lib01_R*.fastq
#scp the resulting .html files to your personal computer to see
The reads appear to tail off in quality toward the ends. This can be fixed with a program called cutadapt.
#run cutadapt without any adapter sequences to cut, but with a PHRED quality cutoff of 30
#(a PHRED score of 30 indicates an expected error rate of 1/1000)
cutadapt --pair-filter=any -q 30 -o Lib01_R1_qced.fastq -p Lib01_R2_qced.fastq Lib01_R1.fastq Lib01_R2.fastq
#now rerun FastQC to see how it worked
mkdir qced_Fastqc_Restults
fastqc -o qced_Fastqc_Restults/ -f fastq Lib01_R*_qced.fastq
Demultiplexing
Now demultiplex the reads. We will do this with process_radtags from the packages STACKS.
#look at the set of barcodes included for our samples
less barcodes_Lib1.tsv
#make a directory to put the resulting sample fastq files into
mkdir sample_fastqs
#look at the documentation for process_radtags
./process_radtags -h
#execute the command to process the rad data
./process_radtags -i 'fastq' -1 Lib01_R1.fastq -2 Lib01_R2.fastq -o ./sample_fastqs/ -b barcodes_Lib1.tsv --inline_index -e 'nlaIII' -r --disable_rad_check
B