What is base calling?

Base-calling converts raw or processed data from a sequencing instrument into sequences and quality scores.

Base-calling usually refers to the conversion of intensity data into sequences and quality scores. Intensity information is extracted from images by the image analysis.

Base-calling has two aspects: Identifying the base-call and assigning a confidence estimate to the call

Quality scores quantify the probability that a base-call is correct (or wrong)

Base quality scores: Individual bases have quality scores which reflect the likelihood of the base being correct/incorrect

Alignment scores: Probability than an alignment to a given position in the reference genome is correct

Allele scores, SNP scores: Probability that a given allele, SNP was observed (often conditional on the alignment being correct)

Base and alignment scores: are single read scores; SNP scores are consensus scores

Consensus calls use information from multiple reads

Phred scores: A base quality score assigned by the phred software (or a program based on the phred which is the most prominent base-calling program for capillary sequencing.)

A quality score expressed on a logarithmic scale:

Q = -10 log10( probability of an error )

Example: Q20 = 1% error probability

The Phred method assigns quality scores to a base-call based on observed properties of the base (predictors)

Phred is a two-step process:

1. Training: Given a set of reads, labels as to which bases are correct, and a set of quality statistics for each base, produce a model that can predict error rates for unseen bases

2. Application: Given new reads and quality statistics, predict the quality for each of the bases.

Phred is essentially a big lookup table!

What goes into a quality table?

Quality predictors are numbers correlated with the quality of a base call, and attempt to quantify concepts such as:

–

“Is the signal for the called base much brighter than the others?”

–

“Did the spot get suspiciously dim, compared to the beginning?”

–

“Does the signal look clean in the next few cycles, and the previous few cycles?

A training data set

–

One of the fundamental questions is: How similar is the training data set to the data set the phred table is applied? Does it make sense for my data?

–

Run-to-run variance, type of sample etc.; need diversity in training

An alignment method and reference

–

phred is inherently based on the assumption that alignments are correct andthe reference is well known

–

Makes it intrinsically very hard to provide accurate scores for low quality bases!

–

Need high-quality references for training –what are the highest qualities we can get for references?

Pileup format

samtools

chr1 808922 G A 21 116 53 9 aAaaaaaAA 9=;:@?9:,

chr1 ---> name of this sequence is chr1

808922 ---> coordinate is 808922 (1-based)

G ---> reference base is 'G'

A ---> consensus base is 'A'

21 ---> consensus quality is 21

116 ---> SNP quality is 116

53 ---> maximum mapping quality is 53

9 ---> the number of reads covering the site is 9

aAaaaaaAA ---> read bases are 'aAaaaaaAA'

9=;:@?9:, --->base qualities are '9=;:@?9:,'

consensus quality (21) is the Phred-scaled probability that the consensus is wrong. SNP quality (116) is the Phred-scaled probability that the consensus is identical to the reference.

For SNP calling, SNP quality is of more importance.

For example, to further filter the SNP calls:

samtools.pl varFilter raw.pileup | awk '$6>=20' > final.pileup

SNP call

novoalign -o SAM -F ILMFQ -t120 -r0.2 -q5 -d hg19.ndx -f 7746x2_s_6_1_sequence.txt.gz 7746x2_s_6_2_sequence.txt.gz | grep chr > 7746x2.sam

samtools view -bS 7746x2.sam | samtools sort - 7746x2

samtools rmdup 7746x2.bam 7746x2.sorted.bam

samtools index 7746x2.sorted.bam

samtools view -u -q 20 7746x2.sorted.bam | samtools pileup -vcf /DataDisk/Analysis/Sun/fa/hg19/chr1.fa - | samtools.pl varFilter -D 100 -G 50 > 1.var

awk '($3=="*"&&$6>=50)||($3!="*"&&$6>30)' 1.var > 1.var.filtered

delete lines from a file

sed -e '2d' hello.txt # delete line #2 from hello.txt

sed -e '1,36d' hello.txt # delete line 1-36 from hello.txt

awk '!(NR>=1 && NR<=36){print $0}' hello.txt # delete line 1-36 from hello.txt

soft clipped

Clipped alignment. In Smith-Waterman alignment, a sequence may not be aligned from the first residue to the last one. Subsequences at the ends may be clipped off. We introduce operation ʻSʼ to describe (softly) clipped alignment. Here is an example. Suppose the clipped alignment is:

REF:AGCTAGCATCGTGTCGCCCGTCTAGCATACGCATGATCGACTGTCAGCTAGTCAGACTAGTCGATCGATGTG

READ: gggGTGTAACC-GACTAGgggg

where on the read sequence, bases in uppercase are matches and bases in lowercase are clipped off. The CIGAR for this alignment is: 3S8M1D6M4S.

A SNP Call pipeline

SNP and MicroIndel detection

Novoalign detects mismatches and short indels in short reads from next-generation sequencing platforms. The recipes below contain command line usage to generate list of single nucleotide polymorphisms SNPs and Indels for further analysis.

Generate sorted alignments in BAM format

Inititally we need to align our file(s) to the reference genome:

novoalign -d hg18.nix -f reads1.fastq -o SAM 2> stats.txt > novo.sam

samtools view -bS novo.sam | samtools sort - reads1

For paired-end or mate-pair reads adjust the parameters accordingly. If we have multiple alignment files these all need to be sorted individually and merged with either samtools or Picard.

Remove PCR duplicates

It is highly recommended to remove PCR duplicates before proceeding with the SNP and Indel calls. We use samtools as an example to remove PCR duplicates.

#remove PCR duplicates from single-end reads

samtools rmdup -s reads1.bam reads1.sorted.bam

#For Paired reads

samtools rmdup reads1.bam reads1.pe.sorted.bam

SNP and Microindel calling

In this example we use the samtools pileup SNP caller and accept reads with a minimum mapping quality of 20:

#The hg18.fasta reference genome sequence is required. samtools.pl is also required

samtools view -u -q 20 reads1.bam | samtools pileup -vc -f hg18.fasta - | samtools.pl varFilter -D 1000 -G 50 > variations.list

Filtering and final list generation

A simple awk script can be used to filter our list of variations where the indel SNP quality is greater than 50 or 60 respectively:

awk '($3=="*"&&$6>=50)||($3!="*"&&$6>30)' variations.list > variations.filtered.list

Exporting to other formats

The samtools distribution contains a helper script called sam2vcf.pl that may be used to convert to the Variant Call Format (VCF):

sam2vcf.pl variations.filtered.list -r hg18.fasta > variations.vcf

Walltime estimation

Input: ~2.5G *2 (pair-end)

Genome: ~4G

250MB reads and 4GB genome per node: ~4 hours