3.1 Phred Quality

Essentially, they quality is the probability that the base calling is wrong. If you are an Audio addict, you can also see it as the signal-noise ratio measured in decibels. The first program to use this measure was Phred, so these are called Phred qualities.

In this model, the quality value is the integer part of -10*log(Probability(error)). That is, you take the probability of error, you calculate its logarithm (in base 10), you multiply the result by -10, and then you take the integer part. The following table shows some typical quality values:

P(correct)	P(error)	log10(P(error))	-10 times
90.0%	0.1	-1	10
95.0%	0.05	-1.3	13
98.0%	0.02	-1.7	17
99.0%	0.01	-2	20
99.9%	0.001	-3	30
99.99%	0.0001	-4	40
99.999999%	0.00000001	-8	80

The sequence of bases and their quality values are then stored in a new file. Several file formats can be used. Today the most common one is FASTQ.

3.2 FASTQ format

The result of the base-calling step has to be stored in a file. There are several formats to store this information, including the original phred format. Today the most common way to store predicted bases and their qualities is FASTQ.

Files in FASTQ format are text files. Thus, they are easy to read and understand. One file can contain several reads. Each read is represented by a record in the file. This format is based on the FASTA format, but with some important differences. Let’s look at the next example:

@SEQ_ID
GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT
+
!''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC65

All FASTQ records follow the same pattern, which is this:

@ ID [description]
SEQUENCE
+
quality

As in FASTA files, each record stats with a marker character. Instead of starting with >, FASTQ records start with @. Like in FASTA, this starter character is immediately followed by the read identification, and optionally a description. Then, like in FASTA, the second line has the nucleotide sequence. In theory, the nucleotide sequence can wrap and use more than one line, but in most cases all the sequence is written in a single line.

The first important difference between FASTA and FASTQ is that the last ones includes qualities. For that, the third line of a FASTQ record is a single + character and the fourth line has the qualities as a text of the same length as the sequence. Each character represents a quality, according to a simple rule: their position in the following lines

!"#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHIJKLMNOPQRSTUVWXYZ[\]^_`abcdefghijklmnopqrstuvwxyz{|}~

Thus, a quality 0 will be represented by !, while a quality of 25 will be represented by :.

3.3 Clipping and Trimming

The last step in this quality control stage is to discard the wrong nucleotides. There are two reasons for doing this.

The first one is to discard low quality letters at the beginning and the end of the read. As we saw in the Figure 3, the first few and the last several bases have lower quality, so they are more noise than signal. It is recommended to clip the reads and remove low quality nucleotides from the 5’ and the 3’ extreme. A reasonable quality value to use as clipping threshold is 10.

The second step is to see if part of the cloning vector or adaptor has been accidentally sequenced. As we know, different sequencing method use different molecular techniques. Most of them use PCR and require specific primer binding sites. These sites are called adaptors. In the Sanger method they are part of the cloning vector. Clearly, they are not part of the genome we want to sequence. Indeed, they are a technical artifact and they should not be used in the next step. Otherwise we risk to confuse the assembler and produce chimeric sequences.

To do vector trimming we need to compare every read to the sequences of known adaptors. We use a small database with the sequences of standard vectors and adaptors.

4.1 Overlap Threshold

To detect that two reads are part of the same genome region, assemblers check if the last nucleotides of one read match the first nucleotides of another read. When this happens, we say that there is an overlap.

It is easy to see that this can induce errors, unless we discard very short overlaps. If the overlap size is small, it is easy to find the same short sequence in two reads just by (bad) chance. If that happens, the genome assembler may erroneously put together reads that belong to different regions of the genome, that are not close in reality. When that happens, we say that there is a chimeric assembly.

Without any additional biological knowledge, we can estimate the probability that two reads overlap on T letters as \[P(\text{overlap length T})=4^{-T}.\] When we have \(N\) reads, we have to compare every one to each other. That is, we have to check \(N(N-1)\approx N^2\) pairs. Combining this two formulas, we can find the expected number of chimeric assemblies for an overlap threshold \(T\) as \[E(\text{number of chimeric assemblies})\approx N^2\cdot 4^{-T}\] For example, if there are \(10^5\) reads, the expected number of chimeras, depending on threshold value \(T\) is

T	Chimeras
5	9765625
10	9537
15	9.313
20	0.009

4.2 Contigs

A contig is a set of reads that overlap. The name comes from the fact that they represent a contiguous part of a chromosome. The assembly program uses these reads to determine their relative position, called layout. Combining reads and layout the program produces a consensus sequence: the bases of each position of the contig, and their quality.

4.3 Scaffolds

In most cases each DNA fragment is sequenced twice, one time for each extreme. We get two reads for each fragment, oriented forward and reverse. Since the reads belong to the same fragment, this gives us some extra information: these two reads are physically close.

When the forward read is located in a different contig than the reverse read, this extra information allows us to know the relative position of each contig. We say that the two contigs are linked by one or more fragments. Following this idea we can connect several contigs into a bigger structure called scaffold.

Finding scaffolds is good for several reasons. First, it contain contigs that belong to the same chromosome, so we get an idea of the chromosome structure. We can use extra information, such as genetic markers or gene linkage to validate or extend the scaffold. The scaffold structure can be validated by doing PCR between contigs that are in theory side by side.

Second, a scaffold is a kind of “big contig with gaps”. The regions between contigs are called gaps. We do not know the value of the bases there, but we know its length, more or less. This allows us to design PCR primers on each side of the gap, inside the contigs, and then amplify and sequence the gap. Normally “gap closing” is the last stage of a sequencing project. Nevertheless many projects, including the human genome, are declared complete even when some gaps still exist in the scaffold.

Some parts of the chromosomes are harder to sequence, such as telomers and centromers. In those cases researchers are happy with an assembly that has one scaffold for each chromosome, and over 99% of the bases are known. The unknown bases are represented by N, as you can see in this example:

>NC_005108.4 Rattus norvegicus strain mixed chromosome 9, Rnor_6.0
TCCACTAACTCTGTTACAAGGAGGGGAGAGAACAGTAAATTGCTTCTGTTACTGAGATCATCTGGAAGTT
GAACTATTTCCTACAGGTGTAACTTTGTTTAGCCCACTGGGAAAGAGATGCTAGATCAGCTCTTGGAAAC
CTTAGCCATAGAGAGAAACCTCAGTGTGAAACTCCAGCTGATGCTTAAACTTCCTCAGTTTCTTGAAAAT
GGTGTGAGCACAGGATCTCCTGATGTTATCTCTTTCCTACCAGTCCTTTACTTGGAGGAATACTCTGCTC
AGCCTGGGTATTTTATGAGTAAACTGACTGGCTATCCCAATGGAGATTCATCAGTAGTTTTTGCTAGATG
TTTGTGAGTCCTCCCCTTGTAATATTATCTTTCTACTGTGATTTGTTCCTTGACCTGGAAGCCGTGCCCA
GGGACACACTATGCTCTGCCTGGGAATCTCAGTGCTGGGAACACAAAGTGTCAGATTCAAGAGGGCTTGT
CCAGGCAATTTCTGCTTGAGGACCTGATTTGTNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNAAGAAATACAGCAGACTTCTCACCAGATTCTATGAATGCCAGAAGGTCCGGGGCAGAT