BACKGROUND
the advent of massively parallel sequencing has the potential to dramatically increase our understanding of genomic and epigenomic variation and of their interaction <cit> .
serving as markers of paternal and maternal chromosomes in heterozygous loci, single-nucleotide polymorphisms  have demonstrated utility to provide information about allele-specific histone marks  <cit>  and to identify differential cpg methylation due to imprinting <cit> .
our understanding of the functional consequences of snps is largely confined to the less than  <dig> % of the genome that codes for amino acid sequences.
increasing our understanding of epigenomically-mediated effects has the potential to elucidate functional consequences of genomic variation within the remaining  <dig> % of the genome  <cit> .
this requires integrative analyses of genomic and epigenomic variation.
pash  <dig>  enables such integrative analyses by achieving the speed required to map in acceptable time the high volumes of reads generated by massively parallel technologies while sensitively detecting dna-sequence level variation in mapped reads.
genome-wide epigenomic assays increasingly utilize massively parallel sequencing instead of microarrays <cit> .
one recent example involved whole-genome bisulfite sequencing to reconstruct two human methylomes  <cit> .
the project involved sequencing a total of  <dig>  billion reads, or  <dig> illumina lanes.
a naïve method would be to examine similarity between every basepair.
when mapping against the  <dig> ×  <dig> nucleotides of the human genome, a total of about  <dig> basepair comparisons would be required.
the gold-standard smith-waterman alignment algorithm <cit> , which performs such basepair-level comparisons, is therefore not practical even if run on the fastest processors.
the still dominant "seed-and-extend" paradigm for fast read mapping emerged during the early sanger sequencing era and has been implemented in comparison tools such as fasta <cit> , blast <cit> , blat <cit>  and ssaha <cit> .
these "seed-and-extend" tools perform filtering of potential similarities using k-mer level matches, called "seeds", and limit basepair-level comparisons to the areas around the seeds, thus reducing the total number of basepair-level comparisons while still performing at an acceptable sensitivity level.
a comprehensive review of early aligners can be found in  <cit> .
the large increase in the number of sequencing reads brought about by massively parallel sequencing required a further increase in comparison speed.
several new aligners such as maq <cit> , bowtie <cit> , bwa <cit> , and eland have initially improved the alignment speed by using one or a combination of heuristics, such as limiting comparison to short reads, performing ungapped alignment, or restricting the number of acceptable differences between the reads and reference genome.
these heuristics have had a generally negative impact on the ability to map reads onto the large fraction of the human genome that is semi-repetitive and to map reads that carry sequence variants not present in the reference sequence, either due to naturally occurring genomic variants, or due to modifications like bisulfite treatment.
newer versions of such aligners have overcome initial limitations, and are able to map long reads containing both basepair substitutions and indels.
for a comprehensive overview of next-generation aligners, we recommend a review by h li and n homer  <cit>
the length of illumina  <cit>  and  <dig>  <cit>  sequencing reads has nearly tripled over the past three years, opening opportunities to map more efficiently onto the large fraction of genomic dna that contains repetitive elements and segmental duplications.
these longer read lengths provide sufficient information for mapping onto polymorphic sites and for detection of sequence variation including indel polymorphisms.
the mapping of bisulfite-treated reads, which contain less information per basepair due to c-t conversion, also benefits from longer read lengths.
several specialized tools have been developed for the mapping of bisulfite treated reads, such as mrsfast <cit> , bsmap <cit> , rmap-bs <cit> , verjinxer <cit> , and brat <cit> .
we compare pash  <dig>  specifically to mrsfast, bsmap, and rmap-bs, and show that pash  <dig>  exhibits speed comparable to that of other bisulfite-seq mapping tools, and has higher sensitivity for mapping of real reads.
to fully exploit the benefit of read lengths exceeding  <dig> basepairs, it is necessary to account for indel-induced gaps, which requires the ability to perform gapped alignments.
instead of performing smith-waterman or other types of alignment by dynamic programming, which is not feasible, pash  <dig>  performs a k-mer level alignment.
in other words, instead of performing alignment by comprehensively examining similarity at every basepair, pash examines multi-basepair similarities, ignoring isolated basepair-level similarities.
specifically, pash examines similarities involving at least k basepairs and, analogous to the dynamic programming alignment algorithm, computes gapped alignments consisting of such similarities.
if the score is within a certain range, banded basepair-level dynamic programming alignment is performed to identify the best mapping.
pash  <dig>  implements k-mer level alignment using multi-positional hash tables.
this data structure generalizes a previous implementation of single-positional hashing  <cit> , which limited each hash table to exactly one sampling position.
by allowing the hashing of any number of sampling positions, the multi-positional hash tables allow greater adaptability to various hardware configurations, from single cpu machines with limited ram and disk space to large clusters consisting of multi-core cpus.
multi-positional hashing is particularly useful for reads-vs-reads and reads-vs-genome mapping, in contrast to the single-positional hashing method employed by pash  <dig>  and  <dig> , which is best suited for genome-to-genome comparisons.
implementation
positional hashing
the seed-and-extend paradigm features prominently in similarity search.
this paradigm originally emerged as a solution for the problem of searching a large database using short queries to detect remote homologies.
the key requirements for such applications were speed and sensitivity across large evolutionary distances.
initial seed-and-extend algorithms such as blast <cit>  and fasta <cit>  relied on query pre-processing, whereas second-generation algorithms such as blat <cit>  and ssaha <cit>  utilized in-memory indexing of genome-sized databases.
in contrast to the seed-and-extend paradigm, positional hashing groups matches between compared sequences.
matches involve multiple collinear basepairs that are obtained by sampling the compared sequences using gapped patterns.
specifically, pattern p is defined by sampled bases {x <dig> ...,xk}.
we say that a match-a gapped k-mer match- is detected between sequences s and t in respective positions i and j if s = t, ..., and s = t. a seed-and-extend method extends a match by local basepair alignment.
in contrast, pash  <dig>  <cit>  and pash  <dig>  <cit>  create a positional hash table for each sampled offset in the reads and group neighboring matches to produce a read mapping score.
a pash user has the option of specifying the gapped pattern.
by default, for certain pattern sizes, pash  <dig>  employs the best discriminating patterns obtained via extensive simulation and search in  <cit> .
while older versions of pash collate the matches, for improved accuracy, pash  <dig>  performs a heuristic alignment of k-mer matches.
efficient positional hashing for reads: multi-positional hash table
previous implementations of pash used distributed positional hash tables, building one separate hash table per sampling position within read.
this approach requires the i/o intensive step of hash table inversion, which involves writing to disk sets of reads that map onto a target genomic window.
for mapping applications using reads produced by next generation technologies, the storage requirement can approach 90- <dig> gb, with an overall disk activity of 1-2tb written and read during one pash  <dig> / <dig>  execution.
to avoid this problem, pash  <dig>  employs multi-positional hash tables, which in effect combine single-positional hashes at the k-mer level, as presented in figure 1a.
pash  <dig>  generalizes the implementation of positional hashing by removing the constraint that each positional hashing table corresponds to exactly one sampling position .
by allowing this multiplicity of offsets, pash performs better in specific applications  while better adapting to a diversity of hardware resources.
when performing read-to-genome mapping, the reads are hashed and pash streams the target reference against the multi-positional hash tables.
k-mer level alignment
a crucial step in pash  <dig>  and  <dig>  was the collation of neighboring matching k-mers.
pash  <dig>  performed ungapped collation of collinear k-mer matches which occur along the same diagonal in the comparison matrix.
pash  <dig>  achieved better sensitivity to indels by collating k-mer matches across a small fixed number of neighboring k diagonals using a greedy algorithm.
pash  <dig>  goes beyond simple collation and employs a heuristic involving the alignment of matching k-mers in a manner inspired by basepair-level dynamic programming, thus achieving far greater speed, as presented in figures 1b and 1c.
specifically, pash  <dig>  aligns matching k-mers between a read of size r and a genomic window of size 2w, with 2wr possible seeds as follows:
let s = k <dig> k <dig>  ...kp the matching k-mers between read r
and genomic window 2w, sorted by offset in the read
let t be the number of sampled bases in each k-mer
let ovl = number of common positions sampled by two k-basepair patterns that overlap by d bases
let m be the reward for a matching base, g the penalty for an indel base
let r = {emptyrun} the set of runs of matching k-mers
bestrun = emptyrun
runextended = false
foreach i =  <dig> ..,p
k = ki
bestextendscore = 0
foreach run r in r
if k extends r then
runextended = true
let d be the overlap length between k and the last k-mer
in run r
let g be the gap size between k and the last k-mer
in run r
newscore = score+)*m - g * gap_size
if  then
bestextendscore = newscore
bestrun = r
end
break
end
end
if  then
r = bestrun
append k to r
score = bestextendscore
end
end
the running time of the algorithm depends on the density of k-mer sampling.
let g be the offset between k-mer sampling positions.
in the worst case, when each k-mer starts a new run, the running time is o2).
note that for good matches, where most k-mers get added to the same run, the complexity becomes o. in practice, this k-mer-level alignment is significantly faster than a basepair-level smith-waterman alignment.
if the k-mer alignment score for a particular genome window and a read r exceeds 3/ <dig> of the current best mapping score for read r, then pash pursues further basepair-level analysis.
the first step is to estimate a bounding box for the alignment score.
for the best mapping of read r, pash evaluates bounds on an approximate alignment score using an affine function, termed a skeleton alignment score.
if the skeleton alignment score exceeds 3/ <dig> of the current best skeleton alignment score for read r, a banded dynamic programming alignment is performed.
finally, if the score obtained by dynamic programming is within a user-specified threshold from the best alignment, pash writes the mapping in a temporary output file.
at the end of the alignment, pash traverses the temporary alignment file and selects the best matches based on the alignment score for each read.
what distinguishes pash from other algorithms that do nominal seed collation, followed by basepair-level alignment, such as blat, is that pash performs the k-mer-level alignment to enable quick filtering of low quality matches.
k-mer-level alignment is a computationally efficient heuristic, since it operates in k-mer space, not at the basepair level.
mapping of bisulfite-treated reads
the k-mer-based framework of pash enables a straightforward extension to mapping of bisulfite-treated reads.
bisulfite sequencing is an accurate method of determining base-level dna methylation status  <cit> .
sample dna is treated with bisulfite, after which high-volume sequencing is employed.
methylated cytosine bases are preserved, but unmethylated ones are converted to uracil, and typically represented as ts in the raw sequence, as shown in figure 2a).
there are several challenges in mapping such reads back to a reference genome: 1) ts can map back to either ts or cs in the reference; 2) reads can come from either the forward or the reverse strand.
in figure 2a we show an example of reads containing both methylated and unmethylated bases, and the effects of applying the bisulfite treatment.
pash constructs all possible k-mers that may arise from bisulfite treatment by converting the ts into either cs or ts, as shown in figure 2b.
the k-mers are hashed in the read hash table, after which regular mapping occurs.
both the forward and the reverse complement strand of each chromosome are used as a reference for mapping.
tradeoffs between speed and sensitivity
pash can trade speed for sensitivity by modifying the sampling pattern , k-mer sampling density based on read length, and by ignoring high frequency k-mers.
pash comes with three built-in heuristics: high sensitivity, medium sensitivity, and low sensitivity; for every heuristic the sampling density is adjusted based on read size.
alternatively, a user can specify a specific sampling density.
another means that pash has for controlling speed is treatment of high-frequency  k-mers.
while in pash  <dig>  this was achieved by imposing a hard limit on the number of entries in any positional hash table bin, in pash  <dig>  a user can specify the percent of hashed k-mers to be used for collation.
typically, 99% of the lowest-frequency k-mers are used and the relatively small fraction of 1% high-frequency k-mers are discarded.
pash has another built-in mode, fast, in which the top 7% highest-frequency k-mers are discarded.
hardware platform and parameter settings
our experimental platform consisted of compute nodes with 8-core intel xeon x <dig> cpus,  <dig>  ghz, and  <dig> gb of ram running linux, kernel  <dig> .
<dig>  we benchmarked pash  <dig> , blat, bwa, bwa-sw, bsmap, mrsfast, rmap-bs, and ssaha <dig>  all experiments were run sequentially; when the input was split into multiple chunks, we reported total compute time.
for pash  <dig>  regular mapping we used the following pattern of weight  <dig> and span 21: 1110110110001101010; for bisulfite mapping we used the pattern of weight  <dig> and span  <dig>  <dig>  by low-level performance analysis, we determined that pash  <dig>  typically spends only 5-10% of its execution time performing basepair level smith-waterman alignment.
the most expensive step, consuming roughly 50-60% of its execution time, is the hash inversion step in which all possible k-mer matches between genomic windows and the reads stored in the multi-positional hash table are considered.
the memory requirements of pash were 2- <dig> gb.
blat used up to  <dig>  gb of ram, bwa and bwa-sw used up to  <dig>  gb of ram, mrsfast and rmap-bs used under  <dig> gb of ram, ssaha <dig> used  <dig> gb of ram, and bsmap used  <dig> gb of ram.
RESULTS
accurate alignment against the non-unique fraction of the genome
we first evaluated the performance of pash by mapping reads against the non-unique portion of the genome, which is increasingly accessible due to increasing read lengths produced by massively parallel technologies.
toward this goal, we developed a benchmark that accounts for segmental duplications, an important feature of mammalian genomes, and one particularly prevalent in the genomes of humans and other primates.
this benchmark was described in  <cit> .
briefly, the ud-csd benchmark models both the unique and the duplicated fraction of a genome; the duplicated reads go through coevolution, gradually acquiring uniqueness by accumulating differences due to independent mutations.
finally, we model the speciation; and then the divergence  of orthologous reads.
the divergence can be viewed either as a measure of species divergence in the case of cross-species alignments, or, in the cases of intra-species alignment, as a measure of the local mutation rate.
in both cases, any sequencing error would also contribute to the "divergence" rate.
the ud-csd benchmark is parameterized by the number of unique tags k; number of duplicated tags n; the rate of coevolution x; and the rate of divergence y.
using this benchmark, we compared the performance of pash  <dig>  to blat, a program well suited for fast mapping of longer reads onto the highly repetitive human genome or across close evolutionary distances .
for our experiment, we chose a read length of  <dig> bases, started with 90% unique reads and 10% repetitive, and varied the total number of reads from  <dig>  to  <dig> , <dig>  after the speciation and divergence simulation, we obtained two read sets, r and s, such that each pair of reads ri and si have a common  <dig> bp ancestor; ri and si have been evolved from the common ancestor by the process of coevolution, followed by speciation, and then independent divergence.
we then employed pash and blat to anchor the read set r onto s, by running each program and then filtering its output such that only best match for each read is kept.
any time a read ri is aligned against its counterpart si we count it as true positive; the overall ratio of true positives with respect to the total number of reads considered represents the true positive rate , and is used as a measure of sensitivity.
in addition, we determine the ratio of true positives out of the total numbers of unique mappings reported, termed positive predictive value , which is a measure of the mapping specificity.
in figure  <dig>  we present the execution times for pash and blat for 25% coevolution and 1% divergence, while in figure  <dig>  we present execution times for pash and blat for 25% coevolution and 5% divergence.
pash was run using a gapped pattern of weight  <dig> and span  <dig>  and a k-mer offset gap of  <dig>  while for blat we used the default settings.
in both cases, the agreement between pash and blat, obtained by comparing the numbers of reads mapped to the expected location by each tool, was 99%; pash and blat had similar ppv rates within 1%.
for up to  <dig> million reads, execution times for pash and blat are comparable.
when the number of reads increases to  <dig>   <dig>   <dig> million reads, however, pash outperforms blat by a factor of  <dig> to  <dig> .
accurate and scaleable alignment of whole genome shotgun sequencing reads
we next compared the performance of pash to an expanded set of aligners, including those specifically designed for massively parallel sequencing applications.
bwa <cit>  is a recently developed aligner that uses a burrows-wheeler transform to index the genome, then performs an index search for matches of reads; for fast execution, a limit on the number of mismatches is set based on the reads size.
bwa-sw  <cit>  is a variant of bwa that was tailored for mapping of long reads.
ssaha <dig> <cit>  is a seed-and-extend aligner, which first builds a "perfect" hash of the genome, and then uses it to seed alignments, and finalizes the alignments using cross_match <cit> .
for our first experiment, we used simulated human whole genome shotgun  reads.
we obtained them by randomly sampling the human genome  using a uniform distribution.
each basepair was then mutated with a probability of  <dig> %, corresponding to the expected human polymorphism rate; 90% of the mutations were basepair substitutions, and 10% were indels in the range 1- <dig> bp.
next, a uniform sequencing error rate of 2% was simulated.
each read dataset was then mapped onto the reference human genome.
to evaluate the performance of mapping, we considered only the reads mapping uniquely, and counted how many of those mapped to a correct location.
table  <dig> summarizes the execution time, true positive rate , and positive predictive value  for pash  <dig> , bwa, bwa-sw, ssaha <dig>  and blat.
execution time, in hours , and percent of reads uniquely and correctly mapped to the original location  are reported.
each dataset contains  <dig> million simulated reads.
the read size varies from  <dig> to  <dig> bp.
pash  <dig>  is run with three different sensitivity settings: high, medium, and low.
bwa and bwa-sw are run with default parameters; ssaha <dig> is run with the - <dig> flag, and blat is run with default parameters .
we report the true positive rate , defined as the overall percentage of correctly mapped reads relative to all the input reads, and the positive predictive value , which indicates what percent of all reported mappings are correct.
tpr is a measure of sensitivity, and ppv is a measure of specificity.
our simulation results show that for reads in the range 76- <dig> bp bwa is the fastest aligner, faster by a factor of  <dig> -8x compared to pash, but pash is slightly more sensitive.
for  <dig> bp, bwa and pash achieve comparable running time, while for reads of length  <dig> bp and  <dig> bp bwa-sw has comparable performance to pash and slightly higher sensitivity.
for long reads, pash is  <dig>  to  <dig>  times faster than ssaha <dig> while achieving comparable sensitivity.
pash outperforms blat up to 191x, and achieves higher sensitivity.
if time is a constraint, then one can use pash in fast mode; when accuracy and sensitivity of mapping are at a premium, pash in high mode should be used.
next, we used actual illumina data and  <dig> reads focusing on comparison between pash, bwa, ssaha <dig>  and blat.
we downloaded a dataset of  <dig>  million  <dig> bp illumina reads from the  <dig> genomes project website , run number srr013932_ <dig>  and mapped it with bwa, ssaha <dig>  and pash  <dig> .
for all runs, we determined the reads mapping uniquely.
bwa imposes a number of mismatches depending on the read size, and requires the entire read to match.
for pash and ssaha <dig>  we required reads to match over a length of at least  <dig> bp, with at most 10% difference between the reads and the reference.
the results are summarized in table  <dig>  bwa is faster by  <dig>  to  <dig>  times compared to pash, but pash maps 23-28% more reads.
note that 96-98% of the pash mappings agree with the ssaha <dig> mappings, whereas only 92% of the bwa mappings agree with the ssaha <dig> mappings.
for the purpose of mapping near-perfect short reads, or reads that are correctly trimmed to eliminate mismatches due to sequencing error, bwa is an extremely fast and accurate aligner, but relative performance deteriorates with increasing read length, and sequence mismatches due to sequencing error rates or the presence of genomic sequence variants.
pash  <dig>  is run with four different sensitivity settings: high, medium, low, and fast.
bwa is run with default parameters.
bwa is  <dig> - <dig> x faster than pash  <dig> , but pash  <dig>  maps 10-14% more of the input reads.
for  <dig> reads, we compared pash against bwa-sw, ssaha <dig> and blat, aligners better suited for long reads.
we used the  <dig> genomes sample srr <dig>  with  <dig>  million  <dig> reads.
in table  <dig> we present the mapping results, including execution time and percent of reads uniquely mapped with a sequence identity of at least 90%.
in high-sensitivity mode, pash runs  <dig>  faster than ssaha <dig>  and achieves comparable sensitivity, mapping only  <dig> % fewer reads.
pash exhibits an effective tradeoff between speed and sensitivity, translating into a graceful sensitivity degradation as the speed increases.
at medium sensitivity it maps 1% fewer reads, for a speed increase of  <dig> x, and at low sensitivity the number of mapped reads drops by  <dig> %, and the execution time decreases  <dig> -fold.
pash in fast mode runs  <dig> x faster than ssaha <dig>  at a sensitivity penalty of  <dig> % for each sensitivity setting of pash, 99- <dig> % of the pash mapping agree with the corresponding ssaha <dig> mappings.
compared to blat, pash is  <dig> to  <dig> times faster;  <dig> - <dig> % of the reads mapped by pash are mapped by blat at the same location.
finally, pash achieves a similar speed to bwa-sw on the  <dig> reads, and maps uniquely  <dig> % more of the input reads compared to bwa-sw.
execution time, in hours, and percent of reads uniquely mapped with at least 90% identity.
pash  <dig>  is run with four different sensitivity settings: high, medium, low, and fast.
ssaha <dig> is run with the - <dig> flag, and blat is run with default parameters .
pash is  <dig>  to  <dig>  times faster than ssaha, with a sensitivity loss of  <dig> - <dig> %.
pash is 15- <dig> times faster than blat.
accurate mapping of bisulfite reads that scales to high-coverage whole-genome bisulfite sequencing
unlike the mappers discussed above, pash  <dig>  also maps bisulfite-treated reads.
to evaluate the performance of pash  <dig> , we compared it with bsmap  <cit> , the first specialized bisulfite mapping software, mrsfast <cit> , and rmap-bs <cit> .
bsmap performs bisulfite sequence mapping by hashing the reference genome, and by hashing multiple seeds around the cpg sites, according to the possible methylation status.
our experiments indicate that bsmap requires a minimum of  <dig> gb of ram for the purpose of detecting cpg methylation.
we attempted to run bsmap by hashing multiple seeds around cs outside of cpg sites as well, but the memory required exceeded  <dig> gb, the maximum memory on the machines that we used for experiments.
because it does not rely on any assumptions about the location of methylated cytosines, pash does not have this limitation.
methylation detection outside of cpgs sites is encountered in embryonic stem cells, and can be also used to estimate bisulfite conversion rates.
mrsfast is another aligner capable of mapping bisulfite reads by performing ungapped mapping.
in contrast, pash is sensitive to gaps, which is important for accurate alignments, especially as the read length increases, and is also important for detecting sequence variants that may affect the methylation state or serve as markers to detect allele-specific epigenomic states.
rmap-bs performs gapped mapping of bisulfite-treated reads.
since bsmap version  <dig>  is limited to mapping reads up to  <dig> bp, a limitation absent in pash  <dig> , we used two datasets of  <dig> basepair reads for the purpose of comparison.
the first dataset was a simulated wgs dataset over the human genome build hg18:  <dig> bp reads were randomly picked from the genome, according to a uniform distribution, a  <dig> % mutation rate was simulated, with 90% of mutations being a basepair substitution and 10% indels in the 1- <dig> basepair range.
for cpg sites, we simulated methylation with probability of 50%; we did not simulate methylation in non-cpg sites due to the inability of bsmap to map against them within available ram.
the second dataset was a subset of the whole genome human bisulfite sequencing data set generated by lister et al.
<cit>  trimmed to  <dig> basepairs.
we ran bsmap using  <dig> basepairs seeds, and converting reference c's to t's only for c's in cpg contexts.
mrsfast and rmap-bs were run with default parameters for bisulfite treated reads.
pash was run using gapped seeds of length  <dig> and weight  <dig>  in table  <dig> we show the execution times for pash, bsmap, mrsfast, and rmap-bs in hours, as well as the true positive rate  and predictive positive value , for the simulated dataset.
pash is  <dig> - <dig> x faster than bsmap; bsmap maps 1% more reads.
mrsfast is  <dig> - <dig> x slower than pash, and maps  <dig> -3% fewer reads than pash.
pash is  <dig> x slower compared to rmap-bs, and rmap-bs maps  <dig> % more reads.
for the trimmed  <dig> basepairs dataset, pash is  <dig> - <dig> x faster than bsmap, and  <dig> - <dig> x faster than mrsfast; and  <dig> x slower compared to rmap-bs.
pash maps 8-11% more actual trimmed reads unambiguously compared to bsmap, mrsfast, and rmap-bs.
each dataset used for evaluation contains  <dig> million reads.
execution time, in hours is reported for the simulated dataset.
we report the true positive rate , defined as the overall percentage of correctly mapped reads relative to all the input reads, and the positive predictive value , which indicates what percent of all reported mappings are correct.
tpr is a measure of sensitivity, and ppv is a measure of specificity.
execution time, in hours, and percent of reads unambiguously mapped are reported for the real trimmed wgs bisulfite-treated reads.
pash  <dig>  is run with two sensitivity settings: high and fast.
pash is  <dig> - <dig> x faster than bsmap,  <dig> - <dig> x faster than mrsfast, and  <dig> x slower than rmap-bs.
on the simulated dataset bsmap and rmap-bs are more sensitive, but pash maps 8-11% more real trimmed reads at unique locations compared to the other aligners.
integrative analysis of genomic and epigenomic variation
epigenomic variation is superimposed on genomic variation, which complicates detection.
for example, apparent variation in methylation levels measured in two individuals at a particular basepair position may be either due to the actual variation of methylation levels of a cytosine in that position or, if a snp occurs in that position, it may be due to variation in genotypes between the two individuals.
to further evaluate the performance of pash  <dig>  when simultaneously mapping variation in the genomes and methylomes, we used the data from a recent study on the methylome of h <dig> and imr <dig>  <cit> .
the original publication reports using an in-house bisulfite-seq read mapping method that operates within a 3-letter mapping alphabet, and which relies on the bowtie aligner  <cit> .
we downloaded the same set of reads and mapped it using pash  <dig> .
only the reads mapping uniquely by either pash or the 3-letter alphabet method were selected for further analysis.
we next removed monoclonal reads and constructed a methylation map including c sites in the genome covered by at least  <dig> reads.
overall, pash  <dig>  led to methylation calls at an additional 9% of c sites.
for the sites where both pash and the 3-letter alphabet method made calls, the agreement was  <dig> %.
using bisulfite-mapped reads we next identified dna sequence variation in the h <dig> and imr <dig> genomes.
pash generates read mapping output in the sam/bam format  <cit> , which is gaining wide acceptance.
off-the-shelf tools to infer genomic variations from read mappings in the sam/bam format, such as samtools, are readily available.
we ran samtools to determine the single nucleotide polymorphism  in both the h <dig> and imr <dig> datasets.
we discovered  <dig>  million snps for h <dig>  out of which 88% were found in dbsnp  <cit> , and  <dig>  million snps for imr <dig>  out of which 85% were found in dbsnp.
we computed copy number variants  using the bisulfite-treated mapped reads.
we divided the genome  into  <dig> kbp-segments, then counted the number of reads that mapped to each segment of the genome.
in order to correct for the illumina platform's gc content bias, we calculated gc-content for each window, using the methylation map to adjust for the fact that unmethylated cs are converted to us by bisulfite treatment.
we then binned the gc-content values in increments of  <dig> % and applied loess correction in r. the data was then segmented into discrete copy-number regions using circular binary segmentation  <cit> .
windows that were more than 75% covered by repeat regions, or fall into gaps in the reference genome were removed.
we used the distribution of reads to set thresholds for homozygous or heterozygous deletion and amplification.
finally, we intersected the resulting putative cnvs with data from  <dig> normal individuals deposited in the database of genomic variants  <cit> .
for h <dig>  11% of the repeat masked genome is copy number variable, and  <dig> % of cnvs overlap with known variants; for imr <dig>  10% of the  genome is copy number variable, and 69% of cnv alterations overlap with known variants.
we next performed a two-step integrative analysis of methylation and snp variation.
first, we examined apparent methylation differences between h <dig> and imr <dig>  of all the c sites where h <dig> and imr <dig> differ in apparent methylation levels,  <dig> % were snps, a four times higher polymorphism rate than the expected  <dig> % rate of polymorphism.
in the second step, we focused on two groups of snps in c positions: those within  <dig> known imprinted regions and those outside of the imprinted regions, presented in table  <dig>  we hypothesized that due to the random placement of c/nonc alleles in heterozygotic differentially methylated maternal and paternal imprinted sites in the two cell lines, the ratio of frequency of agreement in methylation levels between h <dig> and imr <dig> in imprinted regions at c/nonc heterozygotic sites vs. c/c homozygotic sites will be significantly less than the same ratio observed in the rest of the genome.
as indicated in table  <dig>  that was in fact the case.
the significance of the ratio differences was  <dig>  by the chi-square test.
the comparison is performed for c/c sites and for c/nonc heterozygous sites.
the ratio of frequency of agreement in methylation levels between h <dig> and imr <dig> in imprinted regions at c/nonc heterozygotic sites vs. c/c homozygotic sites is significantly less than the same ratio observed in the rest of the genome
the methylation maps for h <dig> and imr <dig>  the snps for both cell lines, as well as copy number variants  can be downloaded from http://genboree.org/pash-supp/.
CONCLUSIONS
pash  <dig>  generally matches the accuracy and speed of niche programs for fast mapping of short reads, and matches or exceeds their performance on longer reads generated by new massively parallel sequencing technologies.
by using the full sequence information from increasingly long reads produced by massively parallel sequencing technologies, pash enables analyses of genomic and epigenomic variation within the semi-repetitive fraction of the human genome that will become increasingly accessible as the reads produced by massively parallel sequencing technologies increase in length.
pash  <dig>  uses positional hashing, a simple and transparent method that gives a user full access to all options that control performance and sensitivity.
it can thus be readily reconfigured to accommodate new originally unanticipated applications.
it is an accurate method that has a low error rate and is both efficient and scaleable, to match ever-increasing sequencing throughputs.
by applying k-mer level alignment, a non-seed comparison method, pash  <dig>  maps accurately to polymorphic sites, including indel polymorphisms.
pash  <dig>  also enables epigenomic applications by supporting accurate and fast mapping of bisulfite-treated reads.
we demonstrate the unique capability of pash to map bisulfite reads for the purpose of simultaneously determining variation in methylation levels and genomic sequence.
this capability opens the doors to investigations of allele-specific epigenomic states and also to investigations of the effects of genomic variation on the epigenomic state in trans and in cis.
software availability and requirements
pash  <dig>  runs on linux and other unix systems and is available at http://www.brl.bcm.tmc.edu/pash/pashdownload.rhtml.
the methylation maps for h <dig> and imr <dig>  the snps for both cell lines, as well as copy number variants  can be downloaded from http://genboree.org/pash-supp/.
authors' contributions
cc and am conceived the pash  <dig>  improvements, designed and interpreted computational experiments, and co-wrote the manuscript.
cc implemented pash  <dig>  and performed mapping performance evaluation.
fy contributed to the wgs mapping performance analysis.
zc and rah contributed to the bisulfite-mapping performance analysis.
cm performed cnv analysis and evaluation on the bisulfite data.
all authors have read and approved the final manuscript.

 the advent of massively parallel sequencing has the potential to dramatically increase our understanding of genomic and epigenomic variation and of their interaction. serving as markers of paternal and maternal chromosomes in heterozygous loci, single-nucleotide polymorphisms have demonstrated utility to provide information about allele-specific histone marks and to identify differential cpg methylation due to imprinting. our understanding of the functional consequences of snps is largely confined to the less than % of the genome that codes for amino acid sequences. increasing our understanding of epigenomically-mediated effects has the potential to elucidate functional consequences of genomic variation within the remaining % of the genome. this requires integrative analyses of genomic and epigenomic variation. pash enables such integrative analyses by achieving the speed required to map in acceptable time the high volumes of reads generated by massively parallel technologies while sensitively detecting dna-sequence level variation in mapped reads. genome-wide epigenomic assays increasingly utilize massively parallel sequencing instead of microarrays. one recent example involved whole-genome bisulfite sequencing to reconstruct two human methylomes. the project involved sequencing a total of billion reads, or illumina lanes. a naïve method would be to examine similarity between every basepair. when mapping against the × nucleotides of the human genome, a total of about basepair comparisons would be required. the gold-standard smith-waterman alignment algorithm, which performs such basepair-level comparisons, is therefore not practical even if run on the fastest processors. the still dominant "seed-and-extend" paradigm for fast read mapping emerged during the early sanger sequencing era and has been implemented in comparison tools such as fasta, blast, blat and ssaha. these "seed-and-extend" tools perform filtering of potential similarities using k-mer level matches, called "seeds", and limit basepair-level comparisons to the areas around the seeds, thus reducing the total number of basepair-level comparisons while still performing at an acceptable sensitivity level. a comprehensive review of early aligners can be found in. the large increase in the number of sequencing reads brought about by massively parallel sequencing required a further increase in comparison speed. several new aligners such as maq, bowtie, bwa, and eland have initially improved the alignment speed by using one or a combination of heuristics, such as limiting comparison to short reads, performing ungapped alignment, or restricting the number of acceptable differences between the reads and reference genome. these heuristics have had a generally negative impact on the ability to map reads onto the large fraction of the human genome that is semi-repetitive and to map reads that carry sequence variants not present in the reference sequence, either due to naturally occurring genomic variants, or due to modifications like bisulfite treatment. newer versions of such aligners have overcome initial limitations, and are able to map long reads containing both basepair substitutions and indels. for a comprehensive overview of next-generation aligners, we recommend a review by h li and n homer the length of illumina and sequencing reads has nearly tripled over the past three years, opening opportunities to map more efficiently onto the large fraction of genomic dna that contains repetitive elements and segmental duplications. these longer read lengths provide sufficient information for mapping onto polymorphic sites and for detection of sequence variation including indel polymorphisms. the mapping of bisulfite-treated reads, which contain less information per basepair due to c-t conversion, also benefits from longer read lengths. several specialized tools have been developed for the mapping of bisulfite treated reads, such as mrsfast, bsmap, rmap-bs, verjinxer, and brat. we compare pash specifically to mrsfast, bsmap, and rmap-bs, and show that pash exhibits speed comparable to that of other bisulfite-seq mapping tools, and has higher sensitivity for mapping of real reads. to fully exploit the benefit of read lengths exceeding basepairs, it is necessary to account for indel-induced gaps, which requires the ability to perform gapped alignments. instead of performing smith-waterman or other types of alignment by dynamic programming, which is not feasible, pash performs a k-mer level alignment. in other words, instead of performing alignment by comprehensively examining similarity at every basepair, pash examines multi-basepair similarities, ignoring isolated basepair-level similarities. specifically, pash examines similarities involving at least k basepairs and, analogous to the dynamic programming alignment algorithm, computes gapped alignments consisting of such similarities. if the score is within a certain range, banded basepair-level dynamic programming alignment is performed to identify the best mapping. pash implements k-mer level alignment using multi-positional hash tables. this data structure generalizes a previous implementation of single-positional hashing, which limited each hash table to exactly one sampling position. by allowing the hashing of any number of sampling positions, the multi-positional hash tables allow greater adaptability to various hardware configurations, from single cpu machines with limited ram and disk space to large clusters consisting of multi-core cpus. multi-positional hashing is particularly useful for reads-vs-reads and reads-vs-genome mapping, in contrast to the single-positional hashing method employed by pash and, which is best suited for genome-to-genome comparisons. implementation positional hashing the seed-and-extend paradigm features prominently in similarity search. this paradigm originally emerged as a solution for the problem of searching a large database using short queries to detect remote homologies. the key requirements for such applications were speed and sensitivity across large evolutionary distances. initial seed-and-extend algorithms such as blast and fasta relied on query pre-processing, whereas second-generation algorithms such as blat and ssaha utilized in-memory indexing of genome-sized databases. in contrast to the seed-and-extend paradigm, positional hashing groups matches between compared sequences. matches involve multiple collinear basepairs that are obtained by sampling the compared sequences using gapped patterns. specifically, pattern p is defined by sampled bases {x...,xk}. we say that a match-a gapped k-mer match- is detected between sequences s and t in respective positions i and j if s = t,..., and s = t. a seed-and-extend method extends a match by local basepair alignment. in contrast, pash and pash create a positional hash table for each sampled offset in the reads and group neighboring matches to produce a read mapping score. a pash user has the option of specifying the gapped pattern. by default, for certain pattern sizes, pash employs the best discriminating patterns obtained via extensive simulation and search in. while older versions of pash collate the matches, for improved accuracy, pash performs a heuristic alignment of k-mer matches. efficient positional hashing for reads: multi-positional hash table previous implementations of pash used distributed positional hash tables, building one separate hash table per sampling position within read. this approach requires the i/o intensive step of hash table inversion, which involves writing to disk sets of reads that map onto a target genomic window. for mapping applications using reads produced by next generation technologies, the storage requirement can approach 90- gb, with an overall disk activity of 1-2tb written and read during one pash / execution. to avoid this problem, pash employs multi-positional hash tables, which in effect combine single-positional hashes at the k-mer level, as presented in figure 1a. pash generalizes the implementation of positional hashing by removing the constraint that each positional hashing table corresponds to exactly one sampling position. by allowing this multiplicity of offsets, pash performs better in specific applications while better adapting to a diversity of hardware resources. when performing read-to-genome mapping, the reads are hashed and pash streams the target reference against the multi-positional hash tables. k-mer level alignment a crucial step in pash and was the collation of neighboring matching k-mers. pash performed ungapped collation of collinear k-mer matches which occur along the same diagonal in the comparison matrix. pash achieved better sensitivity to indels by collating k-mer matches across a small fixed number of neighboring k diagonals using a greedy algorithm. pash goes beyond simple collation and employs a heuristic involving the alignment of matching k-mers in a manner inspired by basepair-level dynamic programming, thus achieving far greater speed, as presented in figures 1b and 1c. specifically, pash aligns matching k-mers between a read of size r and a genomic window of size 2w, with 2wr possible seeds as follows:    let s = k k...kp the matching k-mers between read r       and genomic window 2w, sorted by offset in the read    let t be the number of sampled bases in each k-mer    let ovl = number of common positions sampled by two k-basepair patterns that overlap by d bases    let m be the reward for a matching base, g the penalty for an indel base    let r = {emptyrun} the set of runs of matching k-mers    bestrun = emptyrun    runextended = false       foreach i =..,p       k = ki          bestextendscore = 0       foreach run r in r          if k extends r then             runextended = true             let d be the overlap length between k and the last k-mer                in run r             let g be the gap size between k and the last k-mer                in run r             newscore = score+)*m - g * gap_size             if then                   bestextendscore = newscore                   bestrun = r             end             break          end       end       if then          r = bestrun          append k to r          score = bestextendscore       end    end the running time of the algorithm depends on the density of k-mer sampling. let g be the offset between k-mer sampling positions. in the worst case, when each k-mer starts a new run, the running time is o2). note that for good matches, where most k-mers get added to the same run, the complexity becomes o. in practice, this k-mer-level alignment is significantly faster than a basepair-level smith-waterman alignment. if the k-mer alignment score for a particular genome window and a read r exceeds 3/ of the current best mapping score for read r, then pash pursues further basepair-level analysis. the first step is to estimate a bounding box for the alignment score. for the best mapping of read r, pash evaluates bounds on an approximate alignment score using an affine function, termed a skeleton alignment score. if the skeleton alignment score exceeds 3/ of the current best skeleton alignment score for read r, a banded dynamic programming alignment is performed. finally, if the score obtained by dynamic programming is within a user-specified threshold from the best alignment, pash writes the mapping in a temporary output file. at the end of the alignment, pash traverses the temporary alignment file and selects the best matches based on the alignment score for each read. what distinguishes pash from other algorithms that do nominal seed collation, followed by basepair-level alignment, such as blat, is that pash performs the k-mer-level alignment to enable quick filtering of low quality matches. k-mer-level alignment is a computationally efficient heuristic, since it operates in k-mer space, not at the basepair level. mapping of bisulfite-treated reads the k-mer-based framework of pash enables a straightforward extension to mapping of bisulfite-treated reads. bisulfite sequencing is an accurate method of determining base-level dna methylation status. sample dna is treated with bisulfite, after which high-volume sequencing is employed. methylated cytosine bases are preserved, but unmethylated ones are converted to uracil, and typically represented as ts in the raw sequence, as shown in figure 2a). there are several challenges in mapping such reads back to a reference genome: 1) ts can map back to either ts or cs in the reference; 2) reads can come from either the forward or the reverse strand. in figure 2a we show an example of reads containing both methylated and unmethylated bases, and the effects of applying the bisulfite treatment. pash constructs all possible k-mers that may arise from bisulfite treatment by converting the ts into either cs or ts, as shown in figure 2b. the k-mers are hashed in the read hash table, after which regular mapping occurs. both the forward and the reverse complement strand of each chromosome are used as a reference for mapping. tradeoffs between speed and sensitivity pash can trade speed for sensitivity by modifying the sampling pattern, k-mer sampling density based on read length, and by ignoring high frequency k-mers. pash comes with three built-in heuristics: high sensitivity, medium sensitivity, and low sensitivity; for every heuristic the sampling density is adjusted based on read size. alternatively, a user can specify a specific sampling density. another means that pash has for controlling speed is treatment of high-frequency k-mers. while in pash this was achieved by imposing a hard limit on the number of entries in any positional hash table bin, in pash a user can specify the percent of hashed k-mers to be used for collation. typically, 99% of the lowest-frequency k-mers are used and the relatively small fraction of 1% high-frequency k-mers are discarded. pash has another built-in mode, fast, in which the top 7% highest-frequency k-mers are discarded. hardware platform and parameter settings our experimental platform consisted of compute nodes with 8-core intel xeon x cpus, ghz, and gb of ram running linux, kernel. we benchmarked pash, blat, bwa, bwa-sw, bsmap, mrsfast, rmap-bs, and ssaha all experiments were run sequentially; when the input was split into multiple chunks, we reported total compute time. for pash regular mapping we used the following pattern of weight and span 21: 1110110110001101010; for bisulfite mapping we used the pattern of weight and span by low-level performance analysis, we determined that pash typically spends only 5-10% of its execution time performing basepair level smith-waterman alignment. the most expensive step, consuming roughly 50-60% of its execution time, is the hash inversion step in which all possible k-mer matches between genomic windows and the reads stored in the multi-positional hash table are considered. the memory requirements of pash were 2- gb. blat used up to gb of ram, bwa and bwa-sw used up to gb of ram, mrsfast and rmap-bs used under gb of ram, ssaha used gb of ram, and bsmap used gb of ram.  accurate alignment against the non-unique fraction of the genome we first evaluated the performance of pash by mapping reads against the non-unique portion of the genome, which is increasingly accessible due to increasing read lengths produced by massively parallel technologies. toward this goal, we developed a benchmark that accounts for segmental duplications, an important feature of mammalian genomes, and one particularly prevalent in the genomes of humans and other primates. this benchmark was described in. briefly, the ud-csd benchmark models both the unique and the duplicated fraction of a genome; the duplicated reads go through coevolution, gradually acquiring uniqueness by accumulating differences due to independent mutations. finally, we model the speciation; and then the divergence of orthologous reads. the divergence can be viewed either as a measure of species divergence in the case of cross-species alignments, or, in the cases of intra-species alignment, as a measure of the local mutation rate. in both cases, any sequencing error would also contribute to the "divergence" rate. the ud-csd benchmark is parameterized by the number of unique tags k; number of duplicated tags n; the rate of coevolution x; and the rate of divergence y. using this benchmark, we compared the performance of pash to blat, a program well suited for fast mapping of longer reads onto the highly repetitive human genome or across close evolutionary distances. for our experiment, we chose a read length of bases, started with 90% unique reads and 10% repetitive, and varied the total number of reads from to, after the speciation and divergence simulation, we obtained two read sets, r and s, such that each pair of reads ri and si have a common bp ancestor; ri and si have been evolved from the common ancestor by the process of coevolution, followed by speciation, and then independent divergence. we then employed pash and blat to anchor the read set r onto s, by running each program and then filtering its output such that only best match for each read is kept. any time a read ri is aligned against its counterpart si we count it as true positive; the overall ratio of true positives with respect to the total number of reads considered represents the true positive rate, and is used as a measure of sensitivity. in addition, we determine the ratio of true positives out of the total numbers of unique mappings reported, termed positive predictive value, which is a measure of the mapping specificity. in figure we present the execution times for pash and blat for 25% coevolution and 1% divergence, while in figure we present execution times for pash and blat for 25% coevolution and 5% divergence. pash was run using a gapped pattern of weight and span and a k-mer offset gap of while for blat we used the default settings. in both cases, the agreement between pash and blat, obtained by comparing the numbers of reads mapped to the expected location by each tool, was 99%; pash and blat had similar ppv rates within 1%. for up to million reads, execution times for pash and blat are comparable. when the number of reads increases to million reads, however, pash outperforms blat by a factor of to. accurate and scaleable alignment of whole genome shotgun sequencing reads we next compared the performance of pash to an expanded set of aligners, including those specifically designed for massively parallel sequencing applications. bwa is a recently developed aligner that uses a burrows-wheeler transform to index the genome, then performs an index search for matches of reads; for fast execution, a limit on the number of mismatches is set based on the reads size. bwa-sw is a variant of bwa that was tailored for mapping of long reads. ssaha is a seed-and-extend aligner, which first builds a "perfect" hash of the genome, and then uses it to seed alignments, and finalizes the alignments using cross_match. for our first experiment, we used simulated human whole genome shotgun reads. we obtained them by randomly sampling the human genome using a uniform distribution. each basepair was then mutated with a probability of %, corresponding to the expected human polymorphism rate; 90% of the mutations were basepair substitutions, and 10% were indels in the range 1- bp. next, a uniform sequencing error rate of 2% was simulated. each read dataset was then mapped onto the reference human genome. to evaluate the performance of mapping, we considered only the reads mapping uniquely, and counted how many of those mapped to a correct location. table summarizes the execution time, true positive rate, and positive predictive value for pash, bwa, bwa-sw, ssaha and blat. execution time, in hours, and percent of reads uniquely and correctly mapped to the original location are reported. each dataset contains million simulated reads. the read size varies from to bp. pash is run with three different sensitivity settings: high, medium, and low. bwa and bwa-sw are run with default parameters; ssaha is run with the - flag, and blat is run with default parameters. we report the true positive rate, defined as the overall percentage of correctly mapped reads relative to all the input reads, and the positive predictive value, which indicates what percent of all reported mappings are correct. tpr is a measure of sensitivity, and ppv is a measure of specificity. our simulation results show that for reads in the range 76- bp bwa is the fastest aligner, faster by a factor of -8x compared to pash, but pash is slightly more sensitive. for bp, bwa and pash achieve comparable running time, while for reads of length bp and bp bwa-sw has comparable performance to pash and slightly higher sensitivity. for long reads, pash is to times faster than ssaha while achieving comparable sensitivity. pash outperforms blat up to 191x, and achieves higher sensitivity. if time is a constraint, then one can use pash in fast mode; when accuracy and sensitivity of mapping are at a premium, pash in high mode should be used. next, we used actual illumina data and reads focusing on comparison between pash, bwa, ssaha and blat. we downloaded a dataset of million bp illumina reads from the genomes project website, run number srr013932_ and mapped it with bwa, ssaha and pash. for all runs, we determined the reads mapping uniquely. bwa imposes a number of mismatches depending on the read size, and requires the entire read to match. for pash and ssaha we required reads to match over a length of at least bp, with at most 10% difference between the reads and the reference. the results are summarized in table bwa is faster by to times compared to pash, but pash maps 23-28% more reads. note that 96-98% of the pash mappings agree with the ssaha mappings, whereas only 92% of the bwa mappings agree with the ssaha mappings. for the purpose of mapping near-perfect short reads, or reads that are correctly trimmed to eliminate mismatches due to sequencing error, bwa is an extremely fast and accurate aligner, but relative performance deteriorates with increasing read length, and sequence mismatches due to sequencing error rates or the presence of genomic sequence variants. pash is run with four different sensitivity settings: high, medium, low, and fast. bwa is run with default parameters. bwa is - x faster than pash, but pash maps 10-14% more of the input reads. for reads, we compared pash against bwa-sw, ssaha and blat, aligners better suited for long reads. we used the genomes sample srr with million reads. in table we present the mapping results, including execution time and percent of reads uniquely mapped with a sequence identity of at least 90%. in high-sensitivity mode, pash runs faster than ssaha and achieves comparable sensitivity, mapping only % fewer reads. pash exhibits an effective tradeoff between speed and sensitivity, translating into a graceful sensitivity degradation as the speed increases. at medium sensitivity it maps 1% fewer reads, for a speed increase of x, and at low sensitivity the number of mapped reads drops by %, and the execution time decreases -fold. pash in fast mode runs x faster than ssaha at a sensitivity penalty of % for each sensitivity setting of pash, 99- % of the pash mapping agree with the corresponding ssaha mappings. compared to blat, pash is to times faster; - % of the reads mapped by pash are mapped by blat at the same location. finally, pash achieves a similar speed to bwa-sw on the reads, and maps uniquely % more of the input reads compared to bwa-sw. execution time, in hours, and percent of reads uniquely mapped with at least 90% identity. pash is run with four different sensitivity settings: high, medium, low, and fast. ssaha is run with the - flag, and blat is run with default parameters. pash is to times faster than ssaha, with a sensitivity loss of - %. pash is 15- times faster than blat. accurate mapping of bisulfite reads that scales to high-coverage whole-genome bisulfite sequencing unlike the mappers discussed above, pash also maps bisulfite-treated reads. to evaluate the performance of pash, we compared it with bsmap, the first specialized bisulfite mapping software, mrsfast, and rmap-bs. bsmap performs bisulfite sequence mapping by hashing the reference genome, and by hashing multiple seeds around the cpg sites, according to the possible methylation status. our experiments indicate that bsmap requires a minimum of gb of ram for the purpose of detecting cpg methylation. we attempted to run bsmap by hashing multiple seeds around cs outside of cpg sites as well, but the memory required exceeded gb, the maximum memory on the machines that we used for experiments. because it does not rely on any assumptions about the location of methylated cytosines, pash does not have this limitation. methylation detection outside of cpgs sites is encountered in embryonic stem cells, and can be also used to estimate bisulfite conversion rates. mrsfast is another aligner capable of mapping bisulfite reads by performing ungapped mapping. in contrast, pash is sensitive to gaps, which is important for accurate alignments, especially as the read length increases, and is also important for detecting sequence variants that may affect the methylation state or serve as markers to detect allele-specific epigenomic states. rmap-bs performs gapped mapping of bisulfite-treated reads. since bsmap version is limited to mapping reads up to bp, a limitation absent in pash, we used two datasets of basepair reads for the purpose of comparison. the first dataset was a simulated wgs dataset over the human genome build hg18: bp reads were randomly picked from the genome, according to a uniform distribution, a % mutation rate was simulated, with 90% of mutations being a basepair substitution and 10% indels in the 1- basepair range. for cpg sites, we simulated methylation with probability of 50%; we did not simulate methylation in non-cpg sites due to the inability of bsmap to map against them within available ram. the second dataset was a subset of the whole genome human bisulfite sequencing data set generated by lister et al. trimmed to basepairs. we ran bsmap using basepairs seeds, and converting reference c's to t's only for c's in cpg contexts. mrsfast and rmap-bs were run with default parameters for bisulfite treated reads. pash was run using gapped seeds of length and weight in table we show the execution times for pash, bsmap, mrsfast, and rmap-bs in hours, as well as the true positive rate and predictive positive value, for the simulated dataset. pash is - x faster than bsmap; bsmap maps 1% more reads. mrsfast is - x slower than pash, and maps -3% fewer reads than pash. pash is x slower compared to rmap-bs, and rmap-bs maps % more reads. for the trimmed basepairs dataset, pash is - x faster than bsmap, and - x faster than mrsfast; and x slower compared to rmap-bs. pash maps 8-11% more actual trimmed reads unambiguously compared to bsmap, mrsfast, and rmap-bs. each dataset used for evaluation contains million reads. execution time, in hours is reported for the simulated dataset. we report the true positive rate, defined as the overall percentage of correctly mapped reads relative to all the input reads, and the positive predictive value, which indicates what percent of all reported mappings are correct. tpr is a measure of sensitivity, and ppv is a measure of specificity. execution time, in hours, and percent of reads unambiguously mapped are reported for the real trimmed wgs bisulfite-treated reads. pash is run with two sensitivity settings: high and fast. pash is - x faster than bsmap, - x faster than mrsfast, and x slower than rmap-bs. on the simulated dataset bsmap and rmap-bs are more sensitive, but pash maps 8-11% more real trimmed reads at unique locations compared to the other aligners. integrative analysis of genomic and epigenomic variation epigenomic variation is superimposed on genomic variation, which complicates detection. for example, apparent variation in methylation levels measured in two individuals at a particular basepair position may be either due to the actual variation of methylation levels of a cytosine in that position or, if a snp occurs in that position, it may be due to variation in genotypes between the two individuals. to further evaluate the performance of pash when simultaneously mapping variation in the genomes and methylomes, we used the data from a recent study on the methylome of h and imr. the original publication reports using an in-house bisulfite-seq read mapping method that operates within a 3-letter mapping alphabet, and which relies on the bowtie aligner. we downloaded the same set of reads and mapped it using pash. only the reads mapping uniquely by either pash or the 3-letter alphabet method were selected for further analysis. we next removed monoclonal reads and constructed a methylation map including c sites in the genome covered by at least reads. overall, pash led to methylation calls at an additional 9% of c sites. for the sites where both pash and the 3-letter alphabet method made calls, the agreement was %. using bisulfite-mapped reads we next identified dna sequence variation in the h and imr genomes. pash generates read mapping output in the sam/bam format, which is gaining wide acceptance. off-the-shelf tools to infer genomic variations from read mappings in the sam/bam format, such as samtools, are readily available. we ran samtools to determine the single nucleotide polymorphism in both the h and imr datasets. we discovered million snps for h out of which 88% were found in dbsnp, and million snps for imr out of which 85% were found in dbsnp. we computed copy number variants using the bisulfite-treated mapped reads. we divided the genome into kbp-segments, then counted the number of reads that mapped to each segment of the genome. in order to correct for the illumina platform's gc content bias, we calculated gc-content for each window, using the methylation map to adjust for the fact that unmethylated cs are converted to us by bisulfite treatment. we then binned the gc-content values in increments of % and applied loess correction in r. the data was then segmented into discrete copy-number regions using circular binary segmentation. windows that were more than 75% covered by repeat regions, or fall into gaps in the reference genome were removed. we used the distribution of reads to set thresholds for homozygous or heterozygous deletion and amplification. finally, we intersected the resulting putative cnvs with data from normal individuals deposited in the database of genomic variants. for h 11% of the repeat masked genome is copy number variable, and % of cnvs overlap with known variants; for imr 10% of the genome is copy number variable, and 69% of cnv alterations overlap with known variants. we next performed a two-step integrative analysis of methylation and snp variation. first, we examined apparent methylation differences between h and imr of all the c sites where h and imr differ in apparent methylation levels, % were snps, a four times higher polymorphism rate than the expected % rate of polymorphism. in the second step, we focused on two groups of snps in c positions: those within known imprinted regions and those outside of the imprinted regions, presented in table we hypothesized that due to the random placement of c/nonc alleles in heterozygotic differentially methylated maternal and paternal imprinted sites in the two cell lines, the ratio of frequency of agreement in methylation levels between h and imr in imprinted regions at c/nonc heterozygotic sites vs. c/c homozygotic sites will be significantly less than the same ratio observed in the rest of the genome. as indicated in table that was in fact the case. the significance of the ratio differences was by the chi-square test. the comparison is performed for c/c sites and for c/nonc heterozygous sites. the ratio of frequency of agreement in methylation levels between h and imr in imprinted regions at c/nonc heterozygotic sites vs. c/c homozygotic sites is significantly less than the same ratio observed in the rest of the genome the methylation maps for h and imr the snps for both cell lines, as well as copy number variants can be downloaded from http://genboree.org/pash-supp/.  pash generally matches the accuracy and speed of niche programs for fast mapping of short reads, and matches or exceeds their performance on longer reads generated by new massively parallel sequencing technologies. by using the full sequence information from increasingly long reads produced by massively parallel sequencing technologies, pash enables analyses of genomic and epigenomic variation within the semi-repetitive fraction of the human genome that will become increasingly accessible as the reads produced by massively parallel sequencing technologies increase in length. pash uses positional hashing, a simple and transparent method that gives a user full access to all options that control performance and sensitivity. it can thus be readily reconfigured to accommodate new originally unanticipated applications. it is an accurate method that has a low error rate and is both efficient and scaleable, to match ever-increasing sequencing throughputs. by applying k-mer level alignment, a non-seed comparison method, pash maps accurately to polymorphic sites, including indel polymorphisms. pash also enables epigenomic applications by supporting accurate and fast mapping of bisulfite-treated reads. we demonstrate the unique capability of pash to map bisulfite reads for the purpose of simultaneously determining variation in methylation levels and genomic sequence. this capability opens the doors to investigations of allele-specific epigenomic states and also to investigations of the effects of genomic variation on the epigenomic state in trans and in cis. software availability and requirements pash runs on linux and other unix systems and is available at http://www.brl.bcm.tmc.edu/pash/pashdownload.rhtml. the methylation maps for h and imr the snps for both cell lines, as well as copy number variants can be downloaded from http://genboree.org/pash-supp/. authors' contributions cc and am conceived the pash improvements, designed and interpreted computational experiments, and co-wrote the manuscript. cc implemented pash and performed mapping performance evaluation. fy contributed to the wgs mapping performance analysis. zc and rah contributed to the bisulfite-mapping performance analysis. cm performed cnv analysis and evaluation on the bisulfite data. all authors have read and approved the final manuscript. 
pash makes use of gapped k-mer alignment, a non-seed based comparison method, which is implemented using multi-positional hash tables. this allows pash to run on diverse hardware platforms, including individual computers with standard ram capacity, multi-core hardware architectures and large clusters. by exploiting longer read lengths, pash maps reads onto the large fraction of genomic dna that contains repetitive elements and polymorphic sites, including indel polymorphisms. we demonstrate the versatility of pash by analyzing the interaction between cpg methylation, cpg snps, and imprinting based on publicly available whole-genome shotgun bisulfite sequencing data. pash performs sequence comparison and read mapping and can be employed as a module within diverse configurable analysis pipelines, including chip-seq and methylome mapping by whole-genome bisulfite sequencing. pash generally matches the accuracy and speed of niche programs for fast mapping of short reads, and exceeds their performance on longer reads generated by a new generation of massively parallel sequencing technologies. massively parallel sequencing readouts of epigenomic assays are enabling integrative genome-wide analyses of genomic and epigenomic variation.
 massively parallel sequencing readouts of epigenomic assays are enabling integrative genome-wide analyses of genomic and epigenomic variation. pash performs sequence comparison and read mapping and can be employed as a module within diverse configurable analysis pipelines, including chip-seq and methylome mapping by whole-genome bisulfite sequencing.  pash generally matches the accuracy and speed of niche programs for fast mapping of short reads, and exceeds their performance on longer reads generated by a new generation of massively parallel sequencing technologies. by exploiting longer read lengths, pash maps reads onto the large fraction of genomic dna that contains repetitive elements and polymorphic sites, including indel polymorphisms.  we demonstrate the versatility of pash by analyzing the interaction between cpg methylation, cpg snps, and imprinting based on publicly available whole-genome shotgun bisulfite sequencing data. pash makes use of gapped k-mer alignment, a non-seed based comparison method, which is implemented using multi-positional hash tables. this allows pash to run on diverse hardware platforms, including individual computers with standard ram capacity, multi-core hardware architectures and large clusters. 