ACM transactions on information systems (January)

145 Pages·2005·3.068 MB·English

Introduction to Genomic Information Retrieval
Genomic information retrieval (GIR) is an emerging area of intersection
betweeninformationretrieval(IR)andbioinformatics.Bioinformaticianshave
defined their discipline as the “research, development, or application of com-
putational tools and approaches for expanding the use of biological, medical,
behavioral or health data, including those to acquire, store, organize, archive,
analyze, or visualize such data.”1 This definition clearly overlaps with that
of information retrieval. Indeed, IR practitioners can bring substantial exper-
tisetobioinformaticsbyapplyingandinnovatingwithtechniquesusedintext
and multimedia retrieval for storing, organizing, analyzing, and visualizing
data.
The application domain of GIR is broad. Expertise can be applied to the
retrieval and management of collections of genomic data (such as DNA and
proteinsequences),medicalandbiologicalliterature,microarraydata,orother
bioinformaticsdomains.Problemsthatcanbesolvedincludeunderstandingthe
structure,function,androleofproteins,discoveringevidenceforthecausesof
disease, and annotating medical and biological artifacts. However, the skills
neededtoworkonsuchproblemsarebroad:sufficientunderstandingofbiology
and genetics, information retrieval expertise, and skills in fundamental com-
puter science techniques such as compression, string matching, data mining,
anddatastructuresandalgorithms.
Thisspecialsectionservestobringtogetherarticlesthatshowhowdiverse,
important bioinformatics problems can be solved with GIR. As editor of this
section, I was pleased by the response to the call for papers, and particularly
impressedbythebreadthofworksubmittedforreview.
ArticlesInThisSection
Thisspecialsectioncontainstwoarticlesthatdescribesolutionstoverydiffer-
entproblemsinGIR.
Korodi and Tabus describe, to my knowledge, the most effective DNA com-
pressionschemetodate.CompressionofDNAisinterestingforbothpractical
reasons (such as reduced storage and transmission costs) and functional rea-
sons(suchasinferringstructureandfunctionfromcompressionmodels).DNA
consists of four symbols, and so compaction of the ASCII representation to
two bits per symbol is trivial. However, compression to less that two bits per
symbol is much harder: for example, compression with the well-known bzip2
algorithm is ineffective, leading typically to expansion. Prior to the work de-
scribedinKorodiandTabus’sarticle,themosteffectiveschemes(describedin
detail in the article) achieved compression to between 1.41 and 1.86 bits per
symbol.Thisarticledescribesaschemethatisalwaysmoreeffectivethanother
approaches, achieving between 1.36 and 1.82 bits per symbol. Their approach
1ThisisaworkingdefinitionproposedbytheUSNationalInsitutesofHealth.Theworkingdocu-
mentcanbefoundashttp://grants1.nih.gov/grants/bistic/CompuBioDef.pdf.
ACMTransactionsonInformationSystems,Vol.23,No.1,January2005,Pages1–2.
2 • Introduction
isbasedonanormalizedmaximumlikelihoodmodel,withheuristicsthatmake
theschemefastondesktophardware.
ThearticleofSanderetal.proposesandexaminestechniquesforanalyzing
Serial Analysis of Gene Expression (SAGE) data. SAGE data is a new but
importantresourcethatallowsthecomparisonofnormalhumantissues,cells,
and cancers. The work in this article shows that SAGE data is robust, that
is, that data produced from different laboratories is comparable and can be
usedtogetherforlargestudies.Importantly,theworkalsosuggestsinteresting
similaritiesbetweenbrainandovariancancer,andshowsthatbrainandbreast
cancer have strong gene expression signatures. This may be an interesting
direction for further biomedical research. The GIR techniques described and
usedincludemethodsforpreprocessingtheSAGEdata,clusteringandselecting
datasubsetsofinterest,andmethodsforpresentingandunderstandingSAGE
experiments.
ACKNOWLEDGMENTS
Ithanktherefereesfortheirhardworkinreviewingandprovidingfeedbackon
thepapers,insomecases,toverytightdeadlines.Ialsothanktheauthorsofall
submittedpapers.Last,IamgratefultotheEditor-in-Chief,GaryMarchionini,
forsupportingtheGIRspecialsection,andtoAssociateEditorJavedMostafa
forhishands-onsupportandmentoringthroughoutthereviewing,discussion,
andfinalizingofthespecialsection.IhopethatTOISreadersenjoythespecial
sectiononGenomicInformationRetrieval.
HUGHE.WILLIAMS
RMITUniversity
Melbourne,Australia
ACMTransactionsonInformationSystems,Vol.23,No.1,January2005.
An Efficient Normalized Maximum Likelihood
Algorithm for DNA Sequence Compression
GERGELYKORODIandIOANTABUS
TampereUniversityofTechnology
This article presents an efficient algorithm for DNA sequence compression, which achieves the
bestcompressionratiosreportedoveratestsetcommonlyusedforevaluatingDNAcompression
programs.Thealgorithmintroducesmanyrefinementstoacompressionmethodthatcombines:
(1) encoding by a simple normalized maximum likelihood (NML) model for discrete regression,
throughreferencetoprecedingapproximatematchingblocks,(2)encodingbyafirstordercontext
coding and (3) representing strings in clear, to make efficient use of the redundancy sources in
DNA data, under fast execution times. One of the main algorithmic features is the constraint
on the matching blocks to include reasonably long contiguous matches, which not only reduces
significantlythesearchtime,butalsocanbeusedtomodifytheNMLmodeltoexploittheconstraint
forgettingsmallercodelengths.ThealgorithmhandlesthechangingstatisticsofDNAdatainan
adaptivewayandbypredictivelyencodingthematchingpointersitissuccessfulincompressinglong
approximatematches.ApartfromcomparisonwithpreviousDNAencodingmethods,wepresent
compressionresultsfortherecentlypublishedhumangenomedata.
CategoriesandSubjectDescriptors:E.4[CodingandInformationTheory]—Datacompaction
andcompression;G.3[ProbabilityandStatistics]—Correlationandregressionanalysis;Markov
processes;J.3[LifeandMedicalSciences]—Biologyandgenetics
GeneralTerms:Algorithms,Theory
AdditionalKeyWordsandPhrases:Approximatesequencematching,DNAcompression,normal-
izedmaximumlikelihoodmodel
1. INTRODUCTION
Beingthedefinitivecodebookforthetraitsandbiologicalfunctionalbehavior
ofeverylivingorganism,DNAsequencesarethemostchallenginginformation
sourcesthathumanstrytodecipheranduse.Anenormousquantityofdifferent
DNAsequencesexists,thesizeofeachsequencevaryingintherangeofmillions
tobillionsofnucleotides.Inbothscientificandcommercialcommunitiesthereis
intenseactivitytargetedatsequencingtheDNAofmanyspeciesandstudying
thevariabilityofDNAbetweenindividualsofthesamespecies,whichproduces
Authors’address:InstituteofSignalProcessing,TampereUniversityofTechnology,P.O.Box553,
FIN-33101Tampere,Finland;email:{gergely.korodi,ioan.tabus}@tut.fi.
Permissiontomakedigitalorhardcopiesofpartorallofthisworkforpersonalorclassroomuseis
grantedwithoutfeeprovidedthatcopiesarenotmadeordistributedforprofitordirectcommercial
advantageandthatcopiesshowthisnoticeonthefirstpageorinitialscreenofadisplayalong
withthefullcitation.CopyrightsforcomponentsofthisworkownedbyothersthanACMmustbe
honored.Abstractingwithcreditispermitted.Tocopyotherwise,torepublish,topostonservers,
toredistributetolists,ortouseanycomponentofthisworkinotherworksrequirespriorspecific
permissionand/orafee.PermissionsmayberequestedfromPublicationsDept.,ACM,Inc.,1515
Broadway,NewYork,NY10036USA,fax:[phone],[e-mail].
(cid:1)C 2005ACM1046-8188/05/0100-0003$5.00
ACMTransactionsonInformationSystems,Vol.23,No.1,January2005,Pages3–34.
4 • G.KorodiandI.Tabus
huge amounts of information that needs to be stored and communicated to a
large number of people. Therefore there is a great need for fast and efficient
compressionofDNAsequences.
The traditional way to tackle the lossless compression of data is to make
use of statistically significant repetitions of the symbols or groups of symbols
in the data. Many of the commonly used general purpose methods eliminate
the redundancies of exact replicas of a certain subsequence by using efficient
referencingtoitsclosestoccurrenceinthepast.
Unfortunately, these general-purpose compression algorithms typically fail
toattainanysizereductionwhenappliedtoDNAsequences,whichimpliesthat
thesesequencesrarelypossessthekindofregularities(likehighpredictability
ofthenextsymbolbasedonitscontext)thataretypicalfor,sayhumancreated
text.However,itiswellknownthatredundancyiscommonlypresentinDNA,
butitinvolvesmoresubtleforms,suchaspalindromematchesandapproximate
repetitions.TheeffectivehandlingofthesefeaturesisthekeytosuccessfulDNA
compression.
The direct practical advantages of compressing DNA sequences stem from
the already large, and rapidly increasing amount of available data. Although
in recent years mass storage media and devices became easily affordable and
ubiquitous, data compression now has an even more important role in reduc-
ing the costs of data transmission, since the DNA files are typically shared
and distributed over the Internet. Space-efficient representation of the data
reduces the load on FTP service providers, as the transmissions are done
faster, and it also saves costs for clients who access and download the se-
quences. Since the price of transmission is proportional to the sizes of ac-
cessed files, even savings of 10–20% are useful. In this article we demon-
stratethatthissizereductioncanbeeasilyaccomplishedfortypicaleukaryotic
data.
The compression of DNA has an intrinsic value in itself, but additionally it
provides many useful clues to the nature of regularities that are statistically
significantintheDNAdata,indicatinghowdifferentpartsofthesequenceare
inrelationwitheachother,howsensitivethegenomeistorandomchangessuch
ascrossoverandmutation,whattheaveragecompositionofasequenceis,and
wheretheimportantcompositionchangesoccur.Duetoitsmanypotentialuses
asastatisticalanalysistoolforDNAsequences,anyDNAcompressionmethod
isfacedwithallofthechallengesthatmodelingtoolsarefacing:maximizingthe
statistical relevance of the model used, and choosing between various models
basedontheoverallperformance-complexitymerit.
One example of a statistical question easily solved by an efficient compres-
sionschemeiswhetherornotacertainapproximatematchoftwosubsequences
isstatisticallyrelevant.Tosolvethis,onesimplyhastocheckwhetherencod-
ing one subsequence based on a reference to the other subsequence and also
includingthedescriptionofthemismatchesisabletocompresstheoverallDNA
sequenceornot.Tomaximizethelikelihoodofacorrectanswer,thecompression
programthatisusedhastobethebestknowncompressor,otherwisestatisti-
callyrelevantmatcheswillremainundetected.Thisshowsthatdesigningthe
mostefficientcompressorisagoalofparamountimportanceforstatisticaland
ACMTransactionsonInformationSystems,Vol.23,No.1,January2005.
AnEfficientNMLAlgorithmforDNASequenceCompression • 5
biologicalinference[Lietal.2001;Lietal.2003],evenifthegainsoverother
compressors are just marginal when judging only from the storage savings
benefits.
The complexity of the models that are implicitly or explicitly used by each
specificlosslesscompressionmethodforcapturingtheregularitiesinthedata
should be weighted carefully when comparing various compression methods,
due to its impact on various aspects such as: statistical relevance of model
parameters, simplicity of parameter estimation, simplicity of using the model
forencoding,andattainablecompressionrate.
Thepresentarticlediscussesacompressionmethodaimedatachievingvery
efficientcompressionwithafastexecutiontime.Thefoundationsforthiswork
havebeenlaiddowninTabusetal.[2003a],whichisinturnbasedontheNor-
malizedMaximumLikelihoodmodelpresentedinTabusetal.[2003b].Whereas
inTabusetal.[2003a]thesimplicityofthemodelwasoneofthestipulatedobjec-
tives,herewerelaxthisrestrictioninfavorofincreasedcompressionefficiency
and practical feasibility. As a result, in Section 2 we introduce the framework
forouralgorithm,morecomplexthantheoneconsideredinTabusetal.[2003a],
followedbyadetaileddiscussionandanalysisinSection3regardingimplemen-
tationissues,solutionsandperformanceevaluations.
InSection4wepresenttheempiricalevaluationofourmethod,concentrating
on constraints such as compression efficiency, processing speed and memory
requirements.Weshowthatatthetimeofthiswritingourprogramsurpassesin
termsofcompressibilityallpreviouslypublishedresultsonadefactosetoftest
sequences.Thepracticalfeasibilityofourprogramisdemonstratedbyencoding
theentirehumangenome,achievingacompressionrateof1.66bits/base.
1.1 ConnectionswithPreviousWorks
Thefirstcompressionalgorithmdedicatedexclusivelytonucleicacidsequences,
andhencesignificantlysurpassinggeneral-purposemethods,waspublishedby
GrumbachandTahi[1993].Theirprogram,calledBiocompress,usedaLempel-
Ziv style substitutional algorithm [Ziv and Lempel 1977] that, apart from the
exactdirectmatches,wasalsoabletodetectcomplementarypalindromes,which
are exact matches of an invertible transform of the subsequence, obtained by
takingthesequenceinreversedorderandreplacingeachbasewithitscomple-
ment(A↔T andC ↔G).Suchpalindromematchesareingeneralasfrequent
asdirectmatches,beingaredundancysourceworthexploitingforcompression.
Inalatermethod,calledBiocompress-2GrumbachandTahi[1994],improved
Biocompress by adding an order-2 context coder followed by arithmetic cod-
ing[Rissanen1976].
A new generation of DNA compression algorithms appeared when practi-
cal use was made of the data redundancy involved in approximate repeti-
tions[Chenetal.2001,2002].Especiallyefficientaretwoalgorithmsintroduced
inChenetal.[2001]byChen,KwongandLi:thefirstversion,GenCompress-1,
againusedtheLempel-Zivalgorithmandanorder-2contextcoder,butherethe
substitutional method only used replacement operations when matching two
subsequences. In the second version, GenCompress-2, insertion and deletion
ACMTransactionsonInformationSystems,Vol.23,No.1,January2005.
6 • G.KorodiandI.Tabus
operations were allowed in addition to the replacements, when looking for an
approximate match. However, this did not prove to offer consistent improve-
mentsoverGenCompress-1,accordingtothereportedresults[Chenetal.2001],
leavingopentheusefulnessofconsideringinsertionsanddeletionsinaddition
to substitutions in their scheme. In both GenCompress-1 and GenCompress-2,
alongwiththepointertothematchedsubsequence(includingoffsetandlength)
theencoderhastoprovidecodesfortheeditoperationsinvolvedintheapprox-
imate match, in order that compression would be lossless. Establishing codes
fortheeditoperationwasanempiricalprocess,andassuch,subjectedtoover-
trainingwithrespecttothetestdatathatwasused.
Anearliersubstitutionalmethod,Cfact[Rivalsetal.1995]usedatwopass
coding strategy: in the first pass the exact matches worth providing a gain
in compression were detected in the overall file, and in the second pass the
encodingwasperformed.Assuch,itcanbethoughtofasaprecursortooneof
thebestDNAcompressionalgorithms,DNACompress,introducedbyChen,Li,
Ma,andTromp [Chenetal.2002],whichalsoemploysatwo-passstrategy,and
isagainbasedonasubstitutional(Lempel-Zivstyle)compressionscheme.Inthe
firstpassaspecializedprogramcalledPatternHunterisusedasapreprocessor
for finding significant approximate repetitions. The second pass then encodes
these by a pointer to their previous occurrence, and it also encodes the edit
operations to correct the approximations of the match, while the rest of the
sequencesarecodedusinganorder-2contextarithmeticcoder.
ThesuccessofDNACompressshowstheimportanceofthemethodusedfor
findingsignificantapproximatematches,whichinonewayoranothershouldbe
partofasubstitutionalcompressionalgorithm.Webrieflymentionafewmajor
methods for aligning sequences, because they are well related to the approxi-
mate matching problem, and furthermore because their principle of search is
worthexploitinginacompressionalgorithmfortheapproximatematchsearch
task.
Locatinginadatabasethesequencesthataresignificantlywellalignedwith
a given sequence was a task first solved in the seventies by exact computing
methodsbasedondynamicprogramming,butthemethodswereslowforsearch-
inginlargedatabasesinapracticaltime.WilburandLipmandevelopedseveral
algorithmsforaligningtwonucleotidesequences,suchasinWilburandLipman
[1983], which was based on matching fixed-length seeds, that is, consecutive
strings of nucleotides. This approach became popular and it was later refined
infurtheralgorithmswiththesameobjective,suchastheFASTAalgorithmby
PearsonandLipman[1988],andthefamilyofBLASTprograms[Altschuletal.
1990], which became very popular ways of comparing one sequence against
all the sequences of a database. A more recent search algorithm attaining
significantly faster speed than BLAST programs at the same level of sensi-
tivity is PatternHunter by Ma, Tromp and Li [Ma et al. 2002]. In contrast
to earlier algorithms, it uses strings of nonconsecutive symbols as seeds for
search.
Allofthementionedmethodsarewelloptimizedforthetaskoffindingalign-
mentscoresforsequencesofarbitrary(notnecessarilythesame)lengths,and
theyoperateintwosteps:firstlocatingapotentialgoodmatchusingaseed,then
ACMTransactionsonInformationSystems,Vol.23,No.1,January2005.
AnEfficientNMLAlgorithmforDNASequenceCompression • 7
subsequently trying to extend the seed sequence at both sides and evaluating
after each extension the matching (or alignment, if gaps are included) score,
to get an approximate match as long as possible. This second part makes the
processevenmorecomputationallyexpensive.Inourcompressionprogramthe
approximatesearchtaskismoreparticular:itneedsonlytocompareblocksof
the same size, and although we use the same principle of “seed” based accel-
eration search, the subsequent operations are less expensive than those used
by the BLAST family of programs, and we optimize them to get a very fast
encodingalgorithm.Atthedecoderside,thesearchprocessisnotneeded,and
hence the decoding time is much faster than the encoding time, which is true
forallsubstitutionalcompressionmethods.
Several other methods that do not use the substitution principle have also
been proposed for DNA compression. Some of them only offer estimates of
the achievable compression, or estimates of the entropy [Loewenstern and
Yianilos 1997; Lanctot et al. 2000]; others are full compression programs.
From the latter category we mention a combination of context tree weight-
ing(CTW)methodwithLempel-ZivcodingcalledCTW+LZ[Matsumotoetal.
2000],whichachievesverygoodcompressionresults,butitisveryslow,anda
Burrows-WheelerTransform-basedcompressor[Adjerohetal.2002],whichat-
tainsonlymodestcompressionresults.Thoughseveralfundamentallydifferent
algorithmshavealsobeenproposedinthe10yearsfollowingtheintroduction
of Biocompress, Lempel-Ziv style algorithms always represented the state-of-
the-artinDNAcompression.
A departure from Lempel-Ziv style algorithms is the method NMLComp
presented in Tabus et al. [2003a], which replaces the Lempel-Ziv method by
encodingwiththeNormalizedMaximumLikelihoodmodelfordiscreteregres-
sion[Tabusetal.2003b].InTabusetal.[2003a]itwasdemonstratedthatthe
NMLmodeliswellsuitedforencodingapproximatematchesbasedonreplace-
mentoperationsandalsothatitcanbeimplementedwithverylowalgorithmic
andcomputationalrequirements.
In our present article we introduce refinements and algorithmic improve-
ments for the method introduced in Tabus et al. [2003a], directed towards
speedinguptheexecutionprocessandobtainingbettercompressionresults.
The kernel of the method was presented in Tabus et al. [2003a] and it con-
sistsofblockparsingtheDNAsequenceandusingoneofthethreeoptionsfor
encoding: the first uses a reference to a previously encoded DNA segment for
conditionally encoding the current block, by means of a simple NML model,
the second uses an order-1 context coder, and the third uses “in clear” repre-
sentation of the current block. The present article refines the kernel scheme,
proposingimprovementsinthewaythereferencetothepreviousmatchisen-
coded,howtheveryexpensivesearchofthepreviousmatchcanbeaccelerated
byconstrainingthematches,withjustasmallpriceinthecodingperformance,
andhowthenewmatchconstraintscanbeusedtoderiveamoreefficientNML
model,whichisaseasytohandlealgorithmicallyaswastheoneinTabusetal.
[2003a].Additionally,moreadaptivepoliciesareintroduced,suchastheselec-
tion of the block size used in parsing, and scalable forgetting factors for the
memory model. As a result of these and other algorithmic improvements, our
ACMTransactionsonInformationSystems,Vol.23,No.1,January2005.
8 • G.KorodiandI.Tabus
program achieves the best compression ratios published to date for a de facto
standardtestset,withveryfastexecution.
1.2 OrganizationoftheArticle
In Section 2 we start reviewing the three encoders that compete for encoding
a block of data, and present all components of Tabus et al. [2003a] which are
necessarytomakethispaperself-contained.TheNMLencoderprincipleisdis-
cussedandsomeimportantissuesareelaboratedinmoredetailthaninTabus
et al. [2003a], for example, the analysis of the threshold for deciding which
stringstoincludeinthenormalizationsum,inSections2.1.3and2.1.4.These-
lectionoftheblocksizebecomesamorerefinedprocess,beingintertwinedwith
decisionsoverotherimportantparameters,forexample,forgettingmemory,in
anoverallalgorithm,whichispresentedinSection2.3.Section3containsother
major algorithmic novelties in the article, first addressing the acceleration of
theapproximatematchsearch,andthenconsideringimprovementsintheper-
formance,whichcanbeachievedbyexploitingtheparticularsearchalgorithm
used. A new version of NML encoding is derived, that uses an additional con-
straint in the block model, namely that in each approximate match found, a
contiguousrunofmatchedcharactersthathasatleasta(given)fixedlengthis
guaranteedtoexist.Thealgorithmicimprovementsareassessedoverartificial
randomsequences(whereexpectedvaluescanbecomputed)andoverabench-
mark of DNA sequences. In Section 4 we compare our results with previous
DNAcompressionprogramsoversomeDNAsequences,thenevaluatevarious
options in our algorithm by tuning it to different constraints, and finally we
presenttheresultsofourprogramwhencompressingallofthe3billionbases
ofhumanchromosomes,revealinganaveragecompressibilityof1.66bits/base
forthewell-specifiedbases.
2. ALGORITHMDESCRIPTION
In this section we present an algorithm with the objective of attaining good
compression of DNA sequences while being practically low-cost, especially on
the side of the decoder. In order to efficiently process various types of DNA
sequences,webuildupourencoderfromthreemodules,eachwithseveralpa-
rametersandtargetingittoaddressdifferenttypesofredundancyintheinput
sequence.
In Sections 2.1 and 2.2 we present these three algorithms, while the final
encoderisdescribedinSection2.3.
2.1 TheNMLModelforDiscreteRegression
TheNMLmodelfordiscreteregressionisapowerfultooltodetecthiddenreg-
ularities in the input data [Tabus et al. 2003b]. Such regularities relevant to
DNAcompressionaretheapproximaterepetitions,bywhichonesequencecan
beefficientlyrepresentedusingknowledgeofasequencethatoccurredearlier,
by replacing, deleting or inserting characters in the earlier one to match the
later[Chenetal.2001].Ourstrategytohandleapproximaterepeatsistofocus
on blocks of relatively small length and study their encoding based on substi-
tution operations alone (not considering insertions and deletions). Thus our
ACMTransactionsonInformationSystems,Vol.23,No.1,January2005.
AnEfficientNMLAlgorithmforDNASequenceCompression • 9
firstobjectiveistomatchtwosequencesonlybypossiblysubstitutingacertain
numberofbasesinthefirstsequence,andtoderiveanefficientrepresentation
for the operations needed to recover the second sequence from the first. This
demand is efficiently accomplished by a simplified version of the NML model,
firstdescribedinTabusetal.[2003a].Theideabehindthisapproachistosplit
theinputsequenceintofixed-lengthblocks,andtofindforeachblockanappro-
priate regressor sequence in the already encoded—thus known—data to code
thatblock.TheregressorsequenceshouldhaveminimalHamming-distanceto
thecurrentblock,amongallpossiblesequencesinthepast.
Byourindependentaddressingoftheblocks,wecanmodellongapproximate
repeats well, since the block concatenation can model substitution operations
(insideblocks),andinsertionsanddeletions(byshiftingtheregressorblocks).
AmorecomplexNMLmodelwouldenableustocaptureandeliminatemore
sophisticated redundancies from the DNA sequence. However, such a model
wouldrequireadifferentchoiceoftheregressordata,whichnaturallyimposes
largercostsforfindingthebestmatchingsequence,bothintermsofspeedand
memoryrequirements.Sinceinpracticalapplications,thedatatobeprocessed
canhaveenormoussizes(atypicalchromosomeofaeukaryoteorganismvaries
in the range of several tens up to several hundreds of millions of base pairs)
the advantage of using complex NML models which may offer only slightly
better compression rates, but which at the same time increase computational
requirementsconsiderablyisnotyetclear.Forthisreasoninthefollowingwe
baseouralgorithmonthesimplerNMLvariantfirstdescribedinTabusetal.
[2003a]. In the present section we give an overview of this algorithm, as well
as improvements that have not been discussed in Tabus et al. [2003a] and a
suitablewayofintegratingthemintoanefficientimplementation.
2.1.1 The Probability Model. From the statistical point of view DNA se-
quencesaremessagess=s0,s1,s2,...,sNs−1 emittedbyasourcewithalphabet
size of M = 4 symbols [Grumbach and Tahi 1993]. In biological literature
thesesymbols,calledbases,arenotedbytheletters A,C,G andT,butforcon-
venience we associate to them the numbers 0,1,2 and 3 in respective order;
thus s ∈ {0,1,2,3}. In our method we split the sequence s into successive,
k
nonoverlappingblocksoflengthn,whereonegenericblocks(cid:1)n,...,s(cid:1)n+n−1 will
be referred to as yn = y0,..., yn−1. For convenience, we drop the reference to
therunningindex(cid:1)fromtheblocknotation,althoughwealwaysconsider ynto
bethe“current”block,with y0 =s(cid:1)n,..., yn−1 =s(cid:1)n+n−1forsomeinteger(cid:1).Each
block ynisencodedwiththehelpofsomeregressorblocksk,...,sk+n−1 written
as xn = x0,...,xn−1 inthealreadyencodedstream,whichdoesnotnecessarily
startatapositionthatisamultipleofn.Wedenotebyθ,theprobabilityofthe
event that the symbol at any position in the current block matches the corre-
sponding symbol in the regressor block. If the symbols do not match, there is
an additional choice of 1 in M −1 to specify the correct base y , and thus the
i
probabilitymodelis
(cid:1)
θ if y = x
P(y |x ;θ)= i i , (1)
i i ψ if y (cid:4)= x
i i
ACMTransactionsonInformationSystems,Vol.23,No.1,January2005.
10 • G.KorodiandI.Tabus
whereψ = 1−θ .Fromnowonweassumethatthesymbolswithinthesequence
M−1
areindependent—thoughthisisnottrueforrealDNAsequences,itdoespro-
vide us a simple model, which turns out to be very efficient for our purposes.
Theprobabilitymodel(1)extendedtoblocksofnsymbolsisthen
(cid:2) (cid:2)
P(yn|xn;θ)=θ in=−01χ(yi=xi)ψ in=−01χ(yi(cid:4)=xi) =θnmψn−nm, (2)
whereχ(condition)is1ifconditionistrue;0otherwise.Thenumberofmatching
bases in blocks xn and yn is indicated by n . Consequently, the maximum
m
likelihood estimate of the parameter θ is θˆ(yn,xn) = nm, and the maximized
n
likelihoodis
P(yn|xn;θˆ(yn,xn))=(cid:3)nm(cid:4)nm(cid:5) n−nm (cid:6)n−nm. (3)
n n(M −1)
Itturnsoutthattherearemanypossibleblocksthataredissimilartotheextent
thatthereisnogaininusingxnforthemasregressor.Collectingonlythe“good”
blocksforwhichNMLisprofitableintothesetY ,wenormalizethemaximized
xn
likelihoodtoobtaintheuniversalmodel
P(yn|xn;θˆ(yn,xn))
Pˆ(yn|xn) = (cid:2)
P(zn|xn;θˆ(zn,xn))
zn∈Yxn (cid:7)nm(cid:8)nm(cid:3) n−nm (cid:4)n−nm
= (cid:2) (M −n1)n−m(cid:5)n(nM(cid:6)−1(cid:7))m(cid:8)m(cid:3) n−m (cid:4)n−m, (4)
m∈(cid:5)n m n n(M−1)
where the second expression makes use of the fact that the “good” blocks are
clearly those with high probability in (3), and since this probability depends
onlyonn ,constrainingY isstraightforwardlydonebyconstrainingn toa
m xn m
set(cid:5) ⊆{0,...,n}.
n
2.1.2 TheNMLCodeLength. From(4)itfollowsthatthecodelengththat
isspentonencodingtheblock yn knowingtheregressorblock xn is
L (n ,(cid:5) )=−log Pˆ(yn|xn). (5)
NML m n 2
(Theargumentsofthecodelength,n and(cid:5) indicatethatthislengthdepends
m n
onlyonthenumberofmatchingbasesandthesettowhichn isconstrained.)
m
Theblockxnismostconvenientlyspecifiedbyapointertoitsstart,andonebit
indicatingifdirectorpalindromematchingisused.Sincethe(cid:1)thblock ynstarts
atposition(cid:1)n,andweexcluderegressorsequencesoverlappingthisblock,the
pointer to xn can be encoded with log ((cid:1)n−n+1) bits. However, to prevent
2
explodingrun-timerequirements,inpracticeitisoftenimportanttointroduce
awindowW ={(cid:1)n−n−w+1,...,(cid:1)n−n}ofsizewinwhichthesearchforthe
starting position of the regressor is carried out. This way the pointer can also
be sent in log w bits, so the number of bits needed to completely specify the
2
block xn iseither1+log ((cid:1)n−n+1)or1+log w,whicheverissmaller.This
2 2
meansthattheoverallnumberofbitstheNML-algorithmspendsonencoding
yn is
L (n ,(cid:5) )=L (n ,(cid:5) )+log min{((cid:1)−1)n+1,w}+1. (6)
1 m n NML m n 2
ACMTransactionsonInformationSystems,Vol.23,No.1,January2005.