Exome Sequencing Analysis

Required database indexes

Databases for Exome mapping can be either of the following

  1. Whole genome sequences
  2. Transcript/gene sequences with introns
  3. Make fasta files of these sequences and index them with novoindex
#Whole Genome database for Illumina/454 mapping
novoindex  hg19.nix hg19.fasta
#Whole Genome database for Illumina/454 mapping
novoindex -c  hg19.ncx hg19.fasta

For SOLiD mapping create a colorspace index

#Exon database for Colorspace mapping
novoindex -c hg19.exons.ncx exons.fasta

In the examples below we use the whole genome sequence database.

Exome command line examples

Some whole-exome targeting strategies may leave a 5′ or 3′ terminal adaptor on the end of the read. These can be removed by Novoalign. It is highly recommended that you install the  samtools package. We use novoalign’s read quality calibration to improve the quality of alignments.
Novoalign can be used with all major NGS platforms i.e. Illumina (paired,single and mate-pair) and Ion Torrent, as well as PacBio reads. PacBio reads may need to be split into 1Kb chunks to align more efficiently

A 3′ adaptor is present on all reads. Strip these off and align the read. Generate SAM output format:

#All alignment statistics are written to STDERR when the SAM output format is  specified
novoalign -k -d hg18.nix  -f Reads1.2.fastq Reads1.2.fastq  -a  -o SAM  "@RG\tID:1111\tSM:mySample\tPL:illumina" 2> stats.txt | samtools -bS -  > Reads1.bam

The code above aligns paired-end reads with the adaptor specified with “-a”. For Default Illumina adaptors use “-a” without any args.

Novoalign will auto-detect the FASTQ encoding in version 3 or higher. Otherwise use “-F STDFQ” for standard Sanger FASTQ format.

For single-end read applications e.g. Ion Torrent data, we can use Novoalign to align the reads. In the case of a 5′ adaptor present on the short read sequences:

#Create BAM output for all read alignments. Reads with the 5' adaptor string are preprocessed to remove the adaptor sequence
novoalign -k -d hg19.nix  -f IonReads.fastq -a -o SAM 2> stats.txt | samtools -bS -  > Reads134.bam

FASTA format is also accepted by Novoalign.

Report all alignments for every read in a given sample

#Use the "-r" option to control the reporting strategy to display all alignments. Speed performance is affected when using this option.
novoalign -d hg19_exons.nix  -f Reads134.fastq -a -r All -o SAM 2> stats.txt | samtools -bS -  > Reads134.bam

Sorting the BAM file with multiple cores

We use Novosort to perform an efficient sort of the BAM file and optionally remove or mark PCR duplicates at the same time. This feature replaces one or more Picard/Samtools calls into 1 command that reduces the number of steps needed for your pipeline. Be sure to use the latest version of Novosort for doing the following.

To sort without removing duplicates, using 24 cores and 14Gb of memory:

novosort  -m 14g -t . -c 24  Reads134.bam -i -o Sorted.bam

To sort and remove duplicates, using 24 cores and 14Gb of memory:

novosort  -m 14g -t . -c 24 --removeduplicates --keeptags Reads134.bam -i -o Sorted.bam

To sort-merge previously sorted files and remove duplicates…

novosort  -m 14g -t . -c 24 --removeduplicates --keeptags lane1.bam lane2.bam lane3.bam  -i -o Sorted.bam

Novosort also handles binary BAM-formatted input from STDIN.

Ion Torrent  (LifeTech) Alignment

Novoalign may be used to align Ion Torrent reads to a reference genome. These are typically single-end reads

novoalign -d hg19.nix  -f Ion_reads.fastq.gz -o SAM  -g 20 -x 6 | samtools view -bS – > Ion_reads.bam 

Align Ion Torrent  reads and write output in BAM format lowering gap open (“-g”)  and extend(“-x”) penalties.


Colorspace (SOLiD) Alignment

Align ABI SOLiD colorspace reads and write output in SAM format

novoalignCS -d hg18.ncx  -f SRR040807.fastq.gz -o SAM  > SRR040807.sam 2> stats.txt

Align ABI SOLiD colorspace reads and write output in SAM format

novoalignCS -d hg18.ncx  -f SRR040807.fastq.gz -o SAM  2> stats.txt | samtools view -bS - >SRR040807.bam

Remove PCR duplicates

It is recommended to remove or mark PCR duplicates before calculating gene/exon coverage. We use novosort to remove PCR and optical duplicates or reads with identical mapping coordinates.

novosort -m 8g -c 9 --removeduplicates --keeptags  SRR040807.bam -o SRR040807.deduped.bam

Finding alignments on target

If we have mapped the sequencing reads to a standard reference genome e.g. Human, Mouse, Arabidopsis then it’s quite likely that gene annotations are available for this genome. Depending on the exome capture method there should exist list of all targets using the same reference genome.  In our example below we use the Nimblegen seqcap annotations for the human genome UCSC hg18 version.

Convert the seqcap annotations to BED format for use with the BedTools package. We use the sorted and deduplicated BAM file generated above to calculate the coverage of exome capture targets.

# Note the BED file must be sorted. intersectBed  is part of the BedTools package
 intersectBed -wa -abam SRR040807.deduped.bam  -b Exome_annotations.bed > SRR040807.ontarget.bam

Note that the BED file used can be any annotation so long as it matches up to the reference sequence database that the reads were mapped to with novoalign.

Calculate gene/exon coverage

In cases where a reference genome is well annotated the annotations may be used to calculate gene coverage. The annotations below are examples of human Refseq gene annotations extracted from the UCSC genome database:

#BED-formatted gene annotations
chr1	1298972	1300443	NM_001127229	0	-	1299043	1299999	0	4	173,446,86,204,	0,270,975,1267,
chr1	1298973	1300425	NM_001127230	0	-	1299043	1299999	0	4	172,446,86,108,	0,269,974,1344,
chr1	1310954	1324581	NM_030937	0	-	1312473	1324549	0	11	1871,93,112,142,105,100,65,121,110,75,320,	0,2065,2242,4518,4747,5054,7684,9682,12521,12885,13307,

We can calculate coverage of our short reads using the BEDTools coverageBed program

#calculate gene coverage
coverageBed -abam SRR040807.ontarget.bam -b gene_annotations.bed > SRR040807.coverage.txt

Other programs such as Picard or GATK may be used to accomplish the same task.

Loading posts...