Before you start the alignment and analysis processes, it is useful to perform some initial quality checks on your raw data. Here we will assume you have data from GSAF's Illumina HiSeq sequencer.

When following along here, please start an idev session (within another ssh session) for running any example commands:

idev -r RNASeq-Tue -m 60 -q normal -A UT-2015-05-18

Leave out -r RNASeq-Mon option if you are not part of the reservation

Illumina sequence data format (FASTQ)

GSAF gives you paired end sequencing data in two matching fastq format files, contining reads for each end sequenced -- for example Sample_ABC_L005_R1.cat.fastq and Sample_ABC_L005_R2.cat.fastq. Each read end sequenced is representd by a 4-line entry in the fastq file.

A 4-line fastq file entry looks like this:

A four-line FASTQ file entry representing one sequence
@HWI-ST1097:104:D13TNACXX:4:1101:1715:2142 1:N:0:CGATGT
GCGTTGGTGGCATAGTGGTGAGCATAGCTGCCTTCCAAGCAGTTATGGGAG
+
=<@BDDD=A;+2C9F<CB?;CGGA<<ACEE*1?C:D>DE=FC*0BAG?DB6
  1. Line 1 is the read identifier, which describes the machine, flowcell, cluster, grid coordinate, end and barcode for the read. Except for the barcode information, read identifiers will be identical for corresponding entries in the R1 and R2 fastq files.
  2. Line 2 is the sequence reported by the machine.
  3. Line 3 is always '+' from GSAF (it can optionally include a sequence description)
  4. Line 4 is a string of Ascii-encoded base quality scores, one character per base in the sequence. For each base, an integer quality score = -10 log(probabilty base is wrong) is calculated, then added to 33 to make a number in the Ascii printable character range.

See the Wikipedia FASTQ format page for more information.


Get your data

To try out exercises, we've provided some sample data on stampede2.

So, go get it!

Get set up for the exercises
cds
cd my_rnaseq_course
cd day_1_partA/fastqc_exercise

ls 
ls data

Exercise: Examine the 2nd sequence in a FASTQ file

What is the 2nd sequence in the file Sample1_R1.fastq?

Use the head command.

head data/Sample1_R1.fastq

If this doesn't work, check what directory you are in currently (pwd) and that you've provided the right path to the file. Tab is your friend!

Exercise: Get the first 10 read IDs in the fastq file

grep for lines starting with @HWI since our reads start with that.

For one of GSAF illumina machine names starts with HWI, so if your data was generated on that machine, the reads should start with @HWI

grep '^@HWI' data/Sample1_R1.fastq |head


Counting sequences

One of the first thing to check is that your fastq files are the same length, and that length is evenly divisible by four. The wc command (word count) using the -l switch to tell it to count lines, not words, is perfect for this:

Using wc -l to count lines
wc -l data/Sample1_R1.fastq

Exercise: Counting FASTQ file lines

How many sequences are in the FASTQ file above?

The wc -l command says there are 16000000 lines. FASTQ files have 4 lines per sequence, so the file has 16,000,000/4 or 4,000,000 sequences.

grep '^@HWI' data/Sample1_R1.fastq |wc -l


Lets move on to assessing the quality of this data...

  • No labels

Confluence Documentation | Web Privacy Policy | Web Accessibility