MaxSSmap: A GPU program for mapping divergent short reads to genomes with the maximum scoring subsequence
Abstract:
Background: Programs based on hash tables and Burrows-Wheeler are very fast for mapping short reads to genomes but have low accuracy in the presence of mismatches and gaps. Such reads can be aligned accurately with the Smith-Waterman algorithm but it can take hours and days to map millions of reads even for bacteria genomes.
Results: We introduce a GPU program called MaxSSmap with the aim of achieving comparable accuracy to Smith-Waterman but with faster runtimes. Similar to most programs MaxSSmap identifies a local region of the genome followed by exact alignment. Instead of using hash tables or Burrows-Wheeler in the first part, MaxSSmap calculates maximum scoring subsequence score between the read and disjoint fragments of the genome in parallel on a GPU and selects the highest scoring fragment for exact alignment. We evaluate MaxSSmap's accuracy and runtime when mapping simulated Illumina {\it E.coli} and human chromosome one reads of different lengths and 10\% to 30\% mismatches with gaps to the {\it E.coli} genome and human chromosome one. We also demonstrate applications on real data by mapping ancient horse DNA reads to modern genomes and unmapped paired reads from NA12878 in 1000 genomes.
Conclusions: We show that MaxSSmap attains comparable high accuracy and low error to fast Smith-Waterman programs yet has much lower runtimes. We show that MaxSSmap can map reads rejected by BWA and NextGenMap with high accuracy and low error much faster than if Smith-Waterman were used. On short read lengths of 36 and 51 both MaxSSmap and Smith-Waterman have lower accuracy compared to at higher lengths. On real data MaxSSmap produces many alignments with high score and mapping quality that are not given by NextGenMap and BWA. The MaxSSmap source code in CUDA and OpenCL is freely available from \url{http://www.cs.njit.edu/usman/MaxSSmap}.

Contact: usman@cs.njit.edu

• Supplementary Material
• MaxSSmap source code
• Simulated data:
  • E.coli simulated reads of length 251 (with simulated base qualities): ecoli_reads.tgz
  • E.coli simulated reads of lengths 36, 51, 76, 100, and 150 (with simulated base qualities): ecoli_alllengths.tgz
  • Human chromosome one simulated reads of length 250 (with simulated base qualities): human_chr1_reads.tgz
• Real data:
  • E.coli genome: ecoli_K12_MG1665.fasta.gz
  • Human chromosome one: hs_ref_GRCh37.p13_chr1.fa.gz
  • Horse genome: Equus_caballus.EquCab2.75.dna.chromosome.1to31.fa.gz
  • Human genome: hs_ref_GRCh37.p13_chr1toY.fa.gz
  • E.coli genome in 512 length fragments (for CUDA-SW++): ecoli_K12_MG1665.fasta.fixed_512fragments.gz
  • Ancient horse bone DNA: SRR111892.fastq.len76.first100000.gz
  • Paired reads from NA12878 in 1000 genomes: SRR016607_1_2.tgz
• Perl scripts:
  • BWA, NextGenMap, and MaxSSmap accuracy and error from Stampy simulated reads
  • SSW accuracy and error from Stampy simulated reads

Citation: Turki Turki and Usman Roshan, MaxSSmap: A GPU program for mapping divergent short reads to genomes with the maximum scoring subsequence, (accepted to BMC Genomics, PDF)