bioRxiv preprint first posted online Sep. 20, 2016; doi: http://dx.doi.org/10.1101/076257. The copyright holder for this preprint (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under a CC-BY-NC-ND 4.0 International license. DynOmics to identify delays and co-expression patterns across time course experiments Jasmin Straube1, Bevan Emma Huang2+, and Kim-Anh Leˆ Cao3*+ 1QFAB,InstituteforMolecularBiosciences,TheUniversityofQueensland,QueenslandBiosciencePrecinct,St Lucia,Australia, 2JanssenResearch&Development,LLC,DiscoverySciences,MenloPark,USA 3TheUniversityofQueenslandDiamantinaInstitute,TranslationalResearchInstitute,Woolloongabba,Australia *[email protected] +theseauthorscontributedequallytothiswork ABSTRACT Dynamicchangesinbiologicalsystemscanbecapturedbymeasuringmolecularexpressionfromdifferentlevels(e.g.,genes andproteins)acrosstime. Integrationofsuchdataaimstoidentifymoleculesthatshowsimilarexpressionchangesovertime; suchmoleculesmaybeco-regulatedandthusinvolvedinsimilarbiologicalprocesses. Combiningdatasourcespresents a systematic approach to study molecular behaviour. It can compensate for missing data in one source, and can reduce falsepositiveswhenmultiplesourceshighlightthesamepathways. However,integrativeapproachesmustaccommodate thechallengesinherentin‘omics’data,includinghigh-dimensionality,noise,andtimingdifferencesinexpression. Ascurrent methodsforidentificationofco-expressioncannotcopewiththislevelofcomplexity,wedevelopedanovelalgorithmcalled DynOmics. DynOmics is based on the fast Fourier transform, from which the difference in expression initiation between trajectoriescanbeestimated. Thisdelaycanthenbeusedtorealignthetrajectoriesandidentifythosewhichshowahigh degreeofcorrelation. Throughextensivesimulations,wedemonstratethatDynOmicsisefficientandaccuratecomparedto existingapproaches. Weconsidertwocasestudieshighlightingitsapplication,identifyingregulatoryrelationshipsacross ‘omics’datawithinanorganismandforcomparativegeneexpressionanalysisacrossorganisms. Introduction High-throughput‘omics’platformssuchastranscriptomics,proteomics,andmetabolomicsenablethesimultaneousmonitoring ofthousandsofbiologicalmolecules(transcripts,proteins,andmetabolites),typicallythroughasinglestaticexperiment.1 The recentdecreaseincostofsuchtechnologicalplatformshasmadepossiblethestudyofdynamicbiologicalprocessesbyinstead quantifymoleculesatseveraltimepoints. Thisallowsdeeperinsightintothebehaviourofthemoleculesinsituationsranging fromdevelopmentalprocessestodrugresponse. Thesetimecourse‘omics’experimentsenabletheidentificationofregulators, andmaygiveabetterunderstandingofthestructureanddynamicsofbiologicalsystems. The statistical analysis of dynamic ‘omics’ experiments is difficult. Applying traditional statistical methods for static experimentsislimited,sinceeachtimepointwillbetreatedasindependent,ignoringpotentiallyimportantcorrelationsbetween samplingtimes. Indeed,realisingthepotentialpowerofferedintimecoursestudiestoinvestigateawidevarietyofchanges isnontrivial. Analyticalchallengesarefurthercomplicatedbynoise, smallsamplesizespertimepoint, andfewsampled timepoints. Inthepastdecade,severalmethodshavebeenproposedtoanalysetimecourse‘omics’data,withaparticular focus on microarray and RNA-Seq data. These methods perform differential expression analysis using spline fitting,2,3 Bayesianmethods,4–7 Gaussianprocesses,8–10 andatwo-stepregressionapproach(maSigPro11). Othermethodsfocuson clusteringexpressionprofilestoidentifyco-expressedtrajectories,e.g.,asubsetofmoleculesforwhichexpressionchanges occursimultaneouslyacrosstime.3,7,12–16 Targetedco-expressionanalysiscanalsobeperformedusingvariousmodel-based applications to retrieve data sets from databases given specific query data.17–19 Finally, a third category of methods was proposedbasedonbiologicalpathwayanalysis20,21seethedetailedreviewofSpiesetal.22 Co-expressionanalysiscanprovide valuableinsightintotheroleofmoleculesduringbiologicalprocesses,23–25 butfacessignificantchallengesindealingwith differenttypesof‘omics’andtheirvariationinmolecularresponsetimes. Thesetimingdifferencesordelaysintheinitiationor suppressionofmoleculeexpressionareacommonphenomenoninbiologyandoccuracrossbothdifferentmolecularlevels andorganisms. Forexample,thestudyofregulatoryprocessesafterenvironmentalchangeshasrevealedthatthereisoftena measurabledelayfromthetimeofsignalintroductiontomolecularresponse.23,24,26–28 Thiscanresultfromdifferencesinthe reactionkineticsbetweenanenzymeanditssubstrate,presenceofaninhibitor,oralteredbindingaffinitiesoftranscription factors. SuchprocessescanbestudiedthroughtimecoursemiRNAandmRNAdata,sincemiRNAsplayanimportantrole bioRxiv preprint first posted online Sep. 20, 2016; doi: http://dx.doi.org/10.1101/076257. The copyright holder for this preprint (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under a CC-BY-NC-ND 4.0 International license. ingenetranslationregulationinmanyorganisms,througheithermRNAtranslationinhibitionormRNAdegradation.29 This abilityofmiRNAstofinetunegeneexpressionandtranslationinabroadrangeofimportantbiologicalprocessesisofbroad interestinmedicine.30–32 WhilecorrelationanalysisisfrequentlyusedtoanalysemiRNA-mRNAtimecoursedata,33,34 it mayhavelimitedpowerinsituationswheredelayeddynamicexpressionchangesofmiRNAsrelativetomRNAhavebeen observed.30,35 Delayscanalsohindergeneexpressioncomparisonsacrossorganisms,sinceevenhighlyconservedprocessesmayvary in timing. The pre-implantation embryonic development (PED) is a highly conserved process across mammals, reflected throughtheprogressionofthesamemorphologicstages.36 Nevertheless,attemptstocomparePEDinmammalsbasedongene expressiondatahavefacedchallengesduetodifferencesintimingofgenomeactivationandregulatoryprocesses.37 Hence ignoringthesedelaysinco-expressionanalysiscanmasktrueassociations;thefirststepshouldinsteadbetodetectandquantify thetimedelaybetweenmolecules. Thiswillenableidentificationoffunctionallyrelatedmoleculesregardlessofthedifferences inthetimingofexpressionchanges,aswellasallowingquantificationofsimilaritiesanddifferencesbetweentheobserved responsesinmoredetail. Todate,veryfewmethodsfortimecourse‘omics’dataaccountfortimedelaybetweenmoleculeexpressionlevels. Aijo et al.10 recently proposed DynB, a set of methods based on Gaussian processes, to quantify RNA-Seq gene expression dynamics. This allows rescaling of time profiles, but only between replicates (i.e., at the sample level) rather than at the moleculeexpressionlevel. ThemostcommonlyusedapproachformoleculesistoconsiderthePearsoncorrelation,34,38despite its obvious limitations for detecting co-expressed molecules when their expression change occurs at different time points. LaggedPearsoncorrelation,a.k.a. Pearsoncross-correlationforlaggedtimeseries,circumventsthislimitationbyintroducing artificialdelaysorlagsinthetimeexpressionprofilesforeverypossibletimeshift. Themethodeventuallyappliesthedelay thatmaximisesthecorrelationwiththeoriginalprofile,butcanbepronetooverestimationofdelay. More sophisticated approaches for time course ‘omics’ data come at the expense of computational cost. Shi et al.39 proposedanprobabilisticmodelbasedonmultipledatasetstabularcombinationstoidentifypairwisetranscriptionfactorand gene(TF-G)pairsunderdifferentexperimentalconditions. Thisapproachhasshowntoreducefalsepositivepredictionsbut requiresatimeconsuminglearningsteponexistingandknownTF-Gpairdata.39 DynamicTimeWarping(DTW)40–42 is analgorithmthatalignsthetimepointsoftwotrajectoriestominimisethedistancebetweenthem. Itcanthereforeidentify similaritiesbetweentrajectorieswhichmayvaryinphaseandspeed. Onevariation,DTW4Omics24 identifiesco-expressed moleculeswithapermutationtest,butthiscanbecomputationallyexpensive. Analternateapproach25 utilisesacombined statisticbasedonHiddenMarkovModels(HMM)andPearsoncorrelation. HMMsaretrainedonasetoftrajectorieswherea distributionofvaluesisconsideredforeachtimepoint. Thisgeneratesaprobabilitytoobserveatrajectoryunderthetrained modelthatcantoleratesmalldelays. Whilepromising,thisapproachcannotdetectlargedelays. Additionally,bothitand DTW4Omicscanonlyidentifypositivelycorrelatedtrajectories,requiringheaviercomputationalcoststoexhaustivelyexplore potentialassociations. While integrating time course experiments from different ‘omics’ functional levels is the key to identifying dynamic molecularinteractions,itschallengingnaturehasthusfarpreventedmuchmethodologicaldevelopment. Difficultieslienot onlyinthecomputationrequiredbycomplexalgorithms,butalsoinvariationintypesofcorrelation,levelsofnoiseinthe expressionprofiles,andthedelaysthemselves. WepresentDynOmics,anovelalgorithmtodetect,estimate,andaccountfordelaysbetween‘omics’timeexpression profiles. ThealgorithmisbasedonthefastFouriertransform(FFT),43whichhasalreadybeenshowntosuccessfullydetect periodicallyexpressedgenesintranscriptomicsexperiments.44–46 BycombiningtheFFTangulardifferencebetweenreference andquerytrajectorieswithlaggedPearsoncorrelation,weareabletocharacterisethedirectionandmagnitudeofdelay,whether thereferenceandqueryarepositivelyornegativelycorrelated. Afteraccountingfortheestimateddelays,similarprofilescan beclusteredforfurtherinsight. SimulationresultsshowthatDynOmicsoutperformscurrentmethodstodetecttimeshift,both intermsofsensitivityandspecificity. Weapplyittotwobiologicalcasestudies: onefocusingontheintegrationofmiRNAs andmRNAsinmouselungdevelopment,andoneontheconservationofgenemolecularprocessesacrossmultipleorganisms (mouse,bovine,human)duringPED.Inbothcases,DynOmicsisabletounraveltimingdifferencesbetween‘omics’functional levels,demonstratingitswideapplicability. DynOmicsisavailableasanRpackageonbitbucket.47 Material and Methods Theexpressionchangesofmoleculesmonitoredintimecourseexperimentsoftenformsimpletemporary,sustainableorcyclic patternsthatcanbemodelledasmixturesofoscillating/cyclicpatternsusingthediscreteFouriertransform(FT).48Weintroduce DynOmics,anovelmethodthatfirstconvertstrajectoriestothefrequencydomainusingtheFFT,fromwhichitextractsthe frequencyofthemaincyclicpattern. Condensingthetrajectorytoinformationonthemainfrequencyisthenusedtoidentify whethertwotrajectoriesarerelatedorassociated,whileignoringthenoiseineachtimeexpressionprofile. 2/12 bioRxiv preprint first posted online Sep. 20, 2016; doi: http://dx.doi.org/10.1101/076257. The copyright holder for this preprint (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under a CC-BY-NC-ND 4.0 International license. FourierTransform Foragiventimeseriesx=(x ,...,x,...x ),measuredattimepointst=1,...T ,theFTdecomposesxintocircularcomponents 1 t T orcyclicpatternsforeachfrequencyk=1,...,T−1as: Xk= 1 T∑−1xte−i2πkTt . (1) T t=0 Astheamplitudeatfrequencyk=0simplydescribesthey-axisoffset(i.e.,theglobaldifferencesofexpressionlevels),this frequencyisnotincludedinouranalysiscontext. Equation(1)canbewrittenwithpolarcoordinateswithrealpartaand imaginarypartbasX =a +b i. k k k (cid:113) Foreachfrequencyk=1,...,T−1wecancalculatetheamplituder ofthecomponentasr = a2+b2. Theamplitude k k k k reflectsthecontributionofthekthcyclicpatterntotheoveralltrajectory,andthepatternwithmaximumamplituder˜ describes k themainshapeofthetimeseries. TheargumentArg(X )istheoffsetofthecyclicpattern,definedas: k  (cid:113) 2arctan a2k+bbk2k−ak, ifak>0orbk(cid:54)=0, Arg(X =a +b i)= 0, ifak>0andbk=0, k k k uπn,defined, iiffaak=<00aannddbbk==00., k k Wecantransformtheargumenttothephaseangle(delay)φ indegreesby: k 180∗Arg(X ) k φ = . k π Together,theamplitudeandphaseangledescribeeachfrequencycomponent,andthesetofthesequantitiesisknownasthe frequencydomainrepresentation. DynOmics WedescribeDynOmics, anovelmethodtoestimatedelaysbetweenareferencexandqueryygiveninfrequencydomain representation. First,weidentifyK asthefrequencyofthepatternwithmaximumamplitudeforx,i.e.,themainreference patternfrequency.Then,forbothxandy,weextractphaseanglesatthisfrequencyφx=φ ,φy=φ anddefine∆ =φx−φy xK yK xy asthedifferencebetweenthephaseangles. InFFTliterature,∆ isoftenexpressedintherangeof[−180,180]. Tosimplify xy representationinDynOmics,when∆ <0,weadd360sothat∆ isintherangeof0to359. ∆ indicatesboththesignof xy xy xy thecorrelationbetweenxandyandthesignofthedelay,asseeninFigure1. Thetrajectoriesxandycanbeeitherpositively (Figure1abf)ornegativelycorrelated(Figure1cde),withadelaythatwerefertoasnegative,i.e.,thereferencexispriorto thequeryy(Figure1be)orpositive,i.e.,thereferencexisdelayedwithrespecttothequeryy(Figure1cf). Specificangular differencecasesincludewhen∆ =0(positivecorrelation,butnodelay,Figure1a)andwhen∆ =180(negativecorrelation, xy xy nodelay,Figure1d). WecanestimatethedelaybetweentwotrajectoriesbasedontheFTfrequency,thelengthofthetime seriesand∆ as xy ∆ xy δ = , xy (360/T) K whereδ rangesfrom0to T. Inordertokeepdelayestimatesandhencesignalshiftsassmallaspossible,wecollapsethese xy K valuestotherangeof[− T , T ]bysettingδ =δ −T when270≤∆ ≤359,andδ =δ − T when90≤∆ ≤270. 4K 4K xy xy K xy xy xy 2K xy Wenotethatforqueryprofileseitherpositivelyornegativelycorrelatedwiththereference,thismeansthatδ <0represents xy positivedelay,andδ >0representsnegativedelay. xy UsinglaggedPearsoncorrelationtoincreaseaccuracyindelayestimation. Thedelayestimatebasedontheangular differencepresentedaboveisbasedonapproximatingboththereferenceandquerybythepatternatthemainfrequencyforthe reference. Thisapproximationworkswellwhenthesignalsfrombothqueryandreferencearedominatedbythatmainpattern andrelatively‘noise-free’. However,whenmultiplefrequencieshavesubstantialcontributionstotheoverallsignal,wecan improvetheestimatebymaximisingthelaggedPearsoncorrelationcoefficientoversmallperturbationsinthedelay. Specifically,letδ =(cid:98)δ (cid:101)denoteourinitialdelayestimate,roundedtotheclosestinteger. LetL beasetoflags{l} 0 xy representing perturbations to this initial delay estimate. For each lag, we construct trajectories x and y by shifting the l l 3/12 bioRxiv preprint first posted online Sep. 20, 2016; doi: http://dx.doi.org/10.1101/076257. The copyright holder for this preprint (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under a CC-BY-NC-ND 4.0 International license. a b c d e f 1.0 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● nsity 00..05 ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● ● In●puxt: Reference e Int−0.5 ● y: Query −1.0 2.5 5.0 7.5 2.5 5.0 7.5 2.5 5.0 7.5 2.5 5.0 7.5 2.5 5.0 7.5 2.5 5.0 7.5 Time Figure1. Relationshipbetweenangulardifferences,correlationanddelayforareferencetrajectoryx(reddots)anda querytrajectoryy(greenline). Thetrajectoriesarea)positivelycorrelatedwithnodelay(∆ =0);b)positivelycorrelated xy withnegativedelay(0<∆ ≤90);c)negativelycorrelatedwithpositivedelay(90<∆ <180);d)negativelycorrelatedwith xy xy nodelay(∆ =180);e)negativelycorrelatedwithnegativedelay(180<∆ <270);f)positivelycorrelatedwithpositive xy xy delay(270≤∆ <360). xy original trajectories so that if l <0, x =x ,...,x and y =y ,...,y ; if l >0, conversely, x =x ,...,x and l 1+|l| T l 1 T−|l| l 1 T−|l| y =y ,...,y . The(lagged)Pearsoncorrelationcoefficientbetweenthetwotrajectoriesx andy isdefinedas: l 1+|l| T l l cor(x,y)= T1∑tT=−1|l|(xlt−x¯l)(ylt−y¯l) , (2) l l (cid:113) (cid:113) T1∑tT=−1|l|(xlt−x¯l)2 T1∑tT=−1|l|(ylt−y¯l)2 wherex¯ (y¯)isthesamplemeanacrosstimepointsforeachtrajectory. Wedeterminetheoptimaldelayforagivensetoflags l l asthatforwhichthePearsoncorrelationcoefficientismaximised: δ∗=argmax|cor(x,y)|, (3) l l l∈L withc∗=cor(xδ∗,yδ∗)thenusedtoassessthestrengthanddirectionoftheassociation. WeconsidertwosetsL ,L oflagswhichrepresentperturbationsofδ : 1 2 0 L ={δ ,−δ } 1 0 0 L ={δ −1,δ ,δ +1}, 2 1 1 1 whereδ istheresultoftheoptimizationoverl∈L . Thesetwooptimisationsthusallowustocomparetheinitialestimate 1 1 withthatintheoppositedirection,andthenwithdelaysinalocalneighbourhood. Whiletheoptimisationsdoincreasethe computationrequired,ourrestrictiontolocalperturbationsminimisestheadditionalcomputationwhileimprovingtheestimate inthepresenceofnoise. SensitivityandspecificityofDynOmicscomparedtootherstudies WecomparedDynOmicsperformancewithcurrent availablemethods(DTW4Omics,24PearsonandlaggedPearsoncorrelation)usingmeasuresofsensitivityandspecificitywhile identifyingassociationsinsimulateddata. Thesimulateddataweregeneratedbasedonsimilarscenariosto25withdifferent parameters. Wesimulateddifferentexpressionpatternswithdifferentdelays(−2,1,0,1,2)aswellasdifferentnoiselevels N (0,σ2);σ =0.1,0.2,0.3,0.5. Wealsosimulateddifferentnumberoftimepoints(7and14timespoints). Anumberof7 timespointscharacterisesbestconventional‘omics’timecourseexperiments. Weobservedthatforthesimulatedscenariowith7timepoints,DynOmicsincreasedsensitivitycomparedtocommonlyused methodsbyatleast8%,whilestillremaininghighlyspecific(>0.9). With14timepoints,allmethodsperformedsimilarly, including DynOmics. Pearson correlation which does not take time delays into account performed the worst in terms of sensitivityinallscenarios,demonstratingthatordinarycorrelationmeasuresarenotsufficienttodetectassociationswhen trajectoriesaredelayed. AdetaileddescriptionofthesimulationstudyandtheresultsisprovidedintheSupportingMaterial SectionA. Implementation and computation time DynOmics is implemented in R and uses the FFT implemented in the function fft()fromthestatsRpackage49forthedecompositionofthetimeseries.DynOmicsutilisestheRpackageparallelto performcalculationonCPUsinparallelwherepossible.DynOmics’computationtimewastestedandcomparedtoDTW4Omics onsimulateddatasetswithseventimepointsandtheLungOrganogenesisstudydescribedbelowwith14timepoints. On simulateddatawithonereferenceand100queriesDynOmicsrequiredtwoseconds,whileDTW4Omicsrequired30seconds. OntheLungOrganogenesisstudytheassociationof50referencesand50queriestookDynOmicsfoursecondscomparedto 600secondsforDTW4Omics. 4/12 bioRxiv preprint first posted online Sep. 20, 2016; doi: http://dx.doi.org/10.1101/076257. The copyright holder for this preprint (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under a CC-BY-NC-ND 4.0 International license. Case studies LungOrganogenesis Description of the study The study of Dong et al.34 investigated the dynamic regulation of miRNAs in mouse lung organogenesisbymeasuringtheexpressionof516miRNAsand45,105mRNAsontwobiologicalreplicatesatseventime points(embryoday12,14,16,18;postnatalday2,10,30)inlungs(GSE21053,AffymetrixMouseGenome4302.0Array). Thedataweanalysedwerepre-processedintheoriginalstudy. Subsequently,alinearmixedeffectmodelsplines(LMMS) modellingapproachdevelopedpreviously3wasusedtoobtainrepresentativetrajectoriesover14equallyspacedtimepoints betweenembryoday12andpostnatalday30. Inadditiontoallowinginterpolationtoevenoutspacingbetweentimepoints, LMMScanhandleunbalanceddesigns-whenthenumberofobservationspertimepointisunequal,oriftherearemissingdata. WefurtherfilteredthedatatoretainonlymiRNAsandmRNAsdeclaredasdifferentiallyexpressedovertimeusinglmmsDE3 (FDR≤0.05). ThefinaldatasetanalysedwithDynOmicsincluded105miRNAsand11,326mRNAs. Analysisstrategy. WecomparedassociationsdetectedbetweenmiRNAandmRNApairsforbothrawandLMMSmodelled trajectories,usingeitherclassicalPearsoncorrelation(onrawandLMMSmodelleddata)orDynOmics(onLMMSmodelled data). MiRNAsareknowntobeabletotargettranscriptionregulatorsandthereforeleadtotheindirectexpressionofmany mRNAsdownstream.30 Inthisstudy,however,wefocusedondirecttargetsofmiRNA,andthereforesoughttoidentifynegative correlations between miRNAs and mRNAs, i.e., increased miRNA expression levels associated to a decreased (inhibited) mRNAexpressionlevels,orviceversa. AssociationsweredeclaredforallmiRNA-mRNApairswhosePearsoncorrelation coefficientwas<−0.9.ThemRNAsassociatedtoagivenmiRNAwerecomparedwithmiRNAtargetspredictedfromsequence similarityfrommicroRNA.org(GoodmirSVRscore,ConservedmiRNA,releaseAugust2010),50TargetScan(release6.2)51 andmiRDB(Version5).52 Weonlycompareddatabaseentrieswithanexactidentifiermatchtotheanalysed105miRNAs, leaving14miRNAsformiRDB,33forTargetScanand86formicroRNA.org. Pathwayenrichmentanalysiswasperformed (cid:114) (cid:114) usingQIAGEN’sIngenuity PathwayAnalysis(IPA ,QIAGENRedwoodCity,www.qiagen.com/ingenuity). MammalianEmbryonicDevelopment Descriptionofthestudy. Xieetal.36 investigatedthedynamicexpressionofhuman,mouseandbovinetranscriptsduring PED.Theexpressionlevelsinhuman(30,283mRNAs),mouse(19,607mRNAs)andbovine(13,898mRNAs)weremonitored duringsixtoeightcomparablecellstages(oozygote(onlybovine),zygote,two-,four-,eight-,16-cell(bovineonly),morula andblastocyte)intwo(human,mouse)andthree(bovine)embryoreplicates(GSE18290,Affymetrixmircoarrays: Mouse Expression430AArray,HumanGenomeU133Plus2.0,BovineGenomeArray). Analysisstrategy. Wefirstconvertedthecellstages(zygotetoblastocyte)intoquantitativetimepoints(onetoseven)for inputintomodelling. Foreachorganism,expressiontrajectoriesweremodelledusingLMMSwith14regularlyspacedtime points. Humantranscriptsweretakenasreferences,withreference-querypairsrestrictedtoorthologoussequenceswithmouse andbovineasspecifiedintheAffymetrixfileHG-U133 Plus 2.na26.ortholog.csv. Sevenhumantranscriptsdidnotmatchany identifierintheorthologyfileandwereremoved. Atotalof81,966orthologoustranscriptpairswereanalysed(48,566mouse, 33,400bovine),wherereferencesand/orqueriesmayhavebeenincludedinmultiplepairs. WeappliedDynOmicstoevery orthologoustranscriptpairtoassessdelaysinexpressionlevelsbetweenorganismsanddeclaredassociationwhentheabsolute correlationexceeded0.9. PathwayenrichmentanalysiswasperformedusingIPA. Results LungOrganogenesis Firstly,focusingonthemiRNAsasreferencetrajectories,wecomparedtheperformanceofPearsoncorrelationontherawand LMMSmodelleddata. Wedefinedtheaverageagreementasthenumberofassociationsidentifiedincommonbetweenthetwo methodsdividedbythenumberofassociationsobservedbyonemethodaveragedoverallmiRNAs(SupportingTableS3). WefoundthatmodellingrepresentativetrajectoriesusingLMMSsubstantiallyincreasedthenumberofassociations,byover 80%comparedtorawdata. Thisislikelyduetotheremovalofnoisewhenmodellingthetrajectories.3 Wenextcomparedthe performanceofPearsoncorrelationwithDynOmicsfortheLMMSmodelleddata. DynOmicsidentifiedonaverage18%more associations,indicatingthatthesimplecorrelationanalysiswasnotsufficienttodetectalldelaysinexpressionbetweenmiRNA andmRNA. Secondly, we analysed the overlap of these putative miRNA targets with the miRNA targets predicted from sequence similarity. SupportingTablesS4-S7summariseforeachmiRNAandeachmethodthenumberofputativetargetsandtheoverlap withthepredictedtargetsfromTargetScan,microRNA.organdmiRDB.Fortherawdata,weobservedlowoverlapbetween predictedandputativetargets(rangingfrom0to0.4%miRDB,1.8%microRNA.org,and4.8%TargetScan). Thenumberof overlapsincreasedfortheLMMSmodelleddatawiththemajorityofmiRNA-mRNApairschangingexpressionsimultaneously 5/12 bioRxiv preprint first posted online Sep. 20, 2016; doi: http://dx.doi.org/10.1101/076257. The copyright holder for this preprint (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under a CC-BY-NC-ND 4.0 International license. (i.e.,adelayof0). However,thepercentageofoverlapwasstillsmall(rangingfrom0to3%miRDB,3.5%microRNA.org,and 4.8%TargetScan;SupportingFigureS5). Finally,weinvestigatedwhethertheputativedelayswereofbiologicalrelevanceformiRNA-mRNApairs. ThreemiRNAs in particular, mu-miR-429, mmu-let-7g, and mmu-miR-134, were associated with a large number of negatively delayed mRNAs,representedinFigure21-3. AnalysisofthesedelayedmRNAsusingIPAidentifiedformmu-miR-429enrichment ofthe‘PhospholipaseCSignaling’pathway(P=1.21×10−14),formmu-let-7gthe‘AxonalGuidanceSignaling’pathway (P=4.0×10−11),andformmu-miR-134the‘MitoticRolesofPolo-Like-Kinase’pathway(P=1.29×10−8). Thesepathways havebeendescribedas(cid:19)beinginvolvedineitherembryonicorlungdevelopment. PhospholipaseCwasassociatedwithfetal lungcellproliferationinrats53 andplaysanimportantroleinorganogenesisandembryonicdevelopment.54 Someaxonal guidancemoleculeslikenetrinshavebeensuspectedtoplayaroleinlungbranching,55 whileEphrinB2orsemaphorin3C werefoundtobeinvolvedinalveolargrowthanddevelopment.56,57 Finally,polo-like-kinases(PLKs)arehighlyconservedin mammalsandareimportantforearlyembryonicdevelopment.58 PLKsareknowntoregulatecellcycleprogressionbutlittleis knownabouttheirroleinlungdevelopment. However,over-expressionofPLKshasbeenassociatedwithmalignancyandpoor prognosisinlungcancer,andPLKsarethereforeatargetfortherapy.59 a1 a1 a3 mmu−miR−429 before mmu−let−7g before mmu−miR−134 before miRNA 2 2 2 original n n n o o o si si si inverted es 0 es 0 es 0 pr (cid:237)(cid:21) pr pr x x x E E E −2 −2 −2 (cid:39)(cid:72)(cid:79)(cid:68)(cid:92)(cid:86)(cid:3)(cid:82)(cid:73)(cid:3) 5 10 5 10 5 10 (cid:80)(cid:53)(cid:49)(cid:36)(cid:3) Time Time Time b1 b2 b3 (cid:237)(cid:20) mmu−miR−429 after mmu−let−7g after mmu−miR−134 after (cid:92) (cid:237)(cid:21) 2 2 2 (cid:237)(cid:22) n n n o o o si si si (cid:237)(cid:23) es 0 es 0 es 0 xpr xpr xpr (cid:237)(cid:24) E E E −2 −2 −2 (cid:237)(cid:25) (cid:237)(cid:23) 0 5 10 0 5 10 0 5 10 Time Time Time Figure2. MiRNAandmRNAexpressionassociationsinLungOrganogenesisstudy. ScaledLMMSmodelled expressionlevels(y-axis)aredepictedovertimein14equallyspacedtimeunitsfromembryoday12topostnatalday30 (x-axis)forthemiRNAsmmu-miR-429,mmu-let-7g,andmmu-miR-134(redlines). Solidlinesdepictactualscaledexpression levels,whiledashedlinesdepictinvertedscaledexpressionlevelstoaccountforthenegativecorrelationwithmRNA.Modelled expressionlevelsofthemRNAsidentifiedasassociatedwitheachmiRNAusingDynOmicsaredisplayed(DynOmics correlation<−0.9,delay<0)a)beforeandb)aftershiftingthetrajectoriesusingtheDynOmicsestimateddelay. Theblue colorgradientreflectstheamountofdelay. MammalianPre-implantationEmbryonicDevelopment WeappliedDynOmicstoidentifydelaysinorthologoustranscriptexpressionofmouseandbovinerelativetohumanduring (cid:237)(cid:25) PED.Foranabsolutecorrelationthresholdof0.9,weidentified32,329(67%)orthologouspairsasbeingassociatedbetween humanandmouse,and26,769(80%)betweenhumanandbovine,summarisedinTable1withrespecttothedifferenttypesof delay. Ofthetranscriptsdisplayingassociati(cid:23)on,weobservedtha(cid:27)tthemajorityofth(cid:20)e(cid:21)mouse(56%)and(cid:20)b(cid:25)ovine(67%)transcripts werenotdelayedcomparedtotheorthologoushumantranscripts(cid:91). Interestingly,20%ofmousetranscripts(comparedto10% inbovine)changedexpressionpriortothehumanorthologoustranscript. Thiscouldreflecttimingdifferencesinthezygote genomeactivationofmousePEDatthegeneexpressionlevel.36,37 PathwayanalysisusingIPAwasperformedonthehumanorthologsforthethreetypesofdelay(negative,nodelayand 6/12 bioRxiv preprint first posted online Sep. 20, 2016; doi: http://dx.doi.org/10.1101/076257. The copyright holder for this preprint (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under a CC-BY-NC-ND 4.0 International license. Table1. OrthologoustranscriptsidentifiedasassociatedbyDynOmics. Number(percentage)ofmouseandbovine transcriptsidentifiedasassociatedwithorthologoushumantranscriptsatanabsolutecorrelationthresholdof0.9. Thenumber ofassociationsaredividedaccordingtodifferenttypesofdelay,indicatingwhetherchangesinexpressionlevelsofthemouse andbovinetranscriptsoccurredpriorto(delay>0),simultaneouslyto(delay=0),orafter(delay<0)expressionchangesofthe orthologoushumantranscript. Delay MousevsHuman(%) BovinevsHuman(%) >0 6,582(20) 2,766(10) 0 18,065(56) 17,906(67) <0 7,682(24) 6,097(23) Total 32,329 26,769 positive)relativetomouseorbovineorthologs. Table2liststhetopthreecanonicalpathwaysidentifiedasenrichedforeach type of delay and organism. The majority of trajectories whose expression levels changed in mouse prior to human were involved in EIF2 Signaling (P=7.94×10−18), mTOR Signaling (P=5.64×10−12) and regulation of eIF4 and p70S6K Signaling(P=5.72×10−11). EIF2SignalingandeIF4andp70S6KSignalingplayanimportantroleintranslationregulation andmTORSignalingisanimportantpathwayinembryonicdevelopment.60 Thesesamepathwayswerealsohighlightedina recentstudyusingRNA-Sequencingtechnologies61onhumanduringearlyembryonicdevelopment(4-cell,8-cell,morula,and blastocytestages). EIF2Signaling(P=1.75×10−25)andtheregulationofeIF4(P=3.48×10−0.9)werealsofoundtobe enrichedinbovine;however,thegenesinvolvedinthesepathwayschangedexpressionafterhumanexpressionchanges. AsanillustrativeexamplewedisplaythetrajectoriesoftheorthologoustranscriptsinvolvedinEIF2Signalinginhuman andmousewithrespecttothetypeofdelay(Figure3). Table2. IPAenrichmentanalysisofhumanorthologsforthreetypesofdelayrelativetomouse/bovinetranscripts. ThetopthreeIPAenrichedpathwaysarelisted. Associatedtranscriptswereanalysedseparatelywithrespecttothedelay: positive(negative)delayindicatesthatthemouseorbovineortholog’sexpressionchangesoccurredpriorto(after)thehuman expressionchanges. Nodelayindicatesthatallexpressionchangesoccurredsimultaneously. Delay compared Organism Pathway Pvalue tohuman EIF2Signaling 7.94×10−18 Mouse mTORSignaling 5.64×10−12 RegulationofeIF4andp70S6KSignaling 5.72×10−11 >0 ProteinUbiquitinationPathway 6.28×10−09 Bovine AmyloidProcessing 4.36×10−08 GlucocorticoidReceptorSignaling 1.03×10−06 RoleofMacrophages,FibroblastsandEndothelialCellsinRheuma- 1.84×10−31 Mouse toidArthritis RoleofOsteoblasts,OsteoclastsandChondrocytesinRheumatoid 6.88×10−27 Arthritis 0 AxonalGuidanceSignaling 1.29×10−22 ProteinKinaseASignaling 7.11×10−16 Bovine ThrombinSignaling 1.44×10−15 AcutePhaseResponseSignaling 4.51×10−15 EphrinReceptorSignaling 8.39×10−13 Mouse MolecularMechanismofCancer 5.26×10−12 BCellReceptorSignaling 1.54×10−10 <0 EIF2Signaling 1.75×10−25 Bovine RegulationofeIF4 3.48×10−09 ProteinUbiquitinationPathway 5.11×10−09 Mouse, EIF2Signaling 1.59×10−17 >0,0,<0 Bovine RegulationofeIF4 5.25×10−10 Acetyl-CoABiosynthesisI(PyruvateDehydrogenaseComplex) 8.71×10−06 7/12 bioRxiv preprint first posted online Sep. 20, 2016; doi: http://dx.doi.org/10.1101/076257. The copyright holder for this preprint (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under a CC-BY-NC-ND 4.0 International license. a Human 1 2 3 2 n o si Reference pres 0 Ref x E −2 5 10 5 10 5 10 Time b Mouse 1 2 3 Delays of orthologs 2 1 n o si 2 es 0 3 pr x 4 E −2 5 6 5 10 5 10 5 10 Time Figure3. EIF2Signaling. Modelledtranscriptsexpressionlevels(scaledforeachtimepointforvisualpurposes,y-axis) withrespecttotime(x-axis)involvedinEIF2Signalinga)inhumanwithb)theirorthologsinmouse(DynOmicscorrelation >0.9,delay>0). Hierarchicalclusteringwasperformedonthehumantranscriptstoextractthreemainexpressionpatternsin EIF2Signaling(a;1-3). Thethreemainpatternsofexpressioninhumansa)werevisualisedinseparateplots(1-3). Themouse expressionprofilesinb)wereseparatedbytheclassificationoftheirhumanorthologs(1-3)andwerecolouredaccordingtothe DynOmicsestimatesofdelay. Wealsoperformedenrichmentanalysesforhumanorthologsforalltranscriptsidentifiedasassociated,acrossallthreetypes ofdelay. WehighlighttheconservedprocessofAcetyl-CoABiosynthesisI,sinceithasnotoccurredintheenrichmentlooking atthedelayedorthologsindividually. Acetyl-CoAlevelswerefoundtoplayaroleintheacetylationofproteinsandmay playaroleinregulationofembryogenesis.62 UsingDynOmics,weidentifieddifferentresponsedynamicsacrossorganisms forfouroutofsixtranscripts(dihydrolipoamidebranchedchaintransacylase(DBT),dihydrolipoamides-acetyltransferase (DLAT),dihydrolipoamidedehydrogenase(DLD)andpyruvatedehydrogenase(Lipoamide)beta(PDHB))thatareconserved andinvolvedinthisprocess(Table3). Discussion Todate,veryfewmethodshavebeendevelopedtointegratetimecourse‘omics’datathatarerobusttodelaysinexpression between co-expressed molecules. The integration task is particularly challenging as the data are often characterised by a highlevelofnoiseandmeasuredonasmallnumberoftimepoints. OuralgorithmDynOmicsaddressesthesechallengesby modellingtimecoursetrajectories,identifyingdelaysandre-aligningtrajectoriestodeterminethedegreeofmutualdependency betweenreferenceandquerytrajectories. Modellingtimecoursetrajectoriesisanimportantstepinthisprocess,asmostmethodsdevelopedtointegratetimecourse data,suchasDTW4Omics24andHMMs25requireasinputonlyasinglevaluepertimepoint.Inthisstudyweusedadata-driven modellingapproachbasedonlinearmixedmodelsplines3tosummarisethetimecoursedataappropriately,toreducenoise, andtointerpolateadditionaltimepointswithinthetimecourse. Wefoundthatwhilethemodellingstepmayremovesome associationsbetweenreferenceandquerytrajectories,e.g.,intheLungOrganogenesiscasestudy,itultimatelyincreasedthe numberoffindingsbyconsiderablyreducingtheamountofnoiseinthedata. Inaddition,modellingeachtrajectoryasanoisy functionoftimeallowsintegrationofdatasetswithdifferenttimeintervalsornumbersoftimepoints,aswedemonstratedin 8/12 bioRxiv preprint first posted online Sep. 20, 2016; doi: http://dx.doi.org/10.1101/076257. The copyright holder for this preprint (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under a CC-BY-NC-ND 4.0 International license. Table3. Acetyl-CoABiosynthesisIorthologoustranscripts. OrthologoustranscriptsidentifiedasassociatedbyDynOmics andinvolvedintheAcetyl-CoABiosynthesisIpathway. Genenames,transcriptIDsinhuman,bovineandmouseareindicated, aswellastheestimatedDynOmicsdelayandthePearsoncorrelationbetweenthereferencetrajectoryandthequerytrajectory aftershiftingbasedontheDynOmicsdelayestimate. Genename TranscriptIDHuman TranscriptIDOrganism Organism DynOmicsDelay PearsonCorrelation DBT 205369 x at BT.18489.1.A1 AT Bovine -2 0.99 DBT 205369 x at 1449118 AT Mouse -5 0.98 DLAT 211150 s at 1426264 AT Mouse 3 0.92 DLAT 211150 s at 1426265 X AT Mouse 3 0.91 DLD 230426 at BT.27889.1.S1 AT Bovine 4 0.99 DLD 230426 at 1423159 AT Mouse 4 0.9 PDHB 208911 s at BT.2973.2.S1 A AT Bovine -2 0.98 PDHB 208911 s at BT.2973.3.A1 AT Bovine 3 0.97 PDHB 208911 s at 1416090 AT Mouse 3 0.97 themammalianembryonicdevelopmentcasestudy. Theselectionofanappropriatethresholdtodefineassociationsbetweenco-expressiontrajectoriesisnottrivial,anddepends onthecharacteristicsofthedatathemselves. Forouranalyses,wespecifiedacorrelationthresholdof0.9,aswewereonly interestedinhighlyconcordantexpressiontrajectories. TheroleofmiRNAsasgeneexpressionregulatorsisanexcitingnewsubjectofstudy,asitisestimatedthattheycontrol one-thirdoftheexpressionofthehumangenome.63 Moreover,sincemiRNAsappeartobethemasterswitchinbiological processes,theyarethetargetoffuturetherapeuticdevelopment.31,32 IntheLungOrganogenesisstudy,themiRNA-mRNA associations that we identified with DynOmics largely did not agree with the predictions from the databases TargetScan, microRNA.organdmiRDB.Onepossiblereasonforthegenerallackofagreementcouldbethatthepredictedtargetswerenot expressedintheexperiment,orwerenottargetedatall. ThosemRNAthatdidagreerepresentedsubtledelaysbetweenmiRNA andmRNAtrajectoriesthatmayindicatehighsequenceaffinities. Theotherassociations,whichincludedlargerdelaysthat werenotidentifiedbystandardcorrelationanalysis,34maynotbeassimilarinsequenceandhencewerenotpredictedasmiRNA targetsinthedatabases.50–52 Indeed,ourresultssuggestthatsequenceinformationalonemaynotsufficetodeterminewhether miRNAsareexpressedandregulatespecificmRNAundercertainconditions. Alternately,thelargenumberofmiRNA-mRNA associationsidentifiedbyDynOmicsmayrepresentmRNAswhichareindirecttargetsofmiRNAs. Determiningwhetherthese mRNAsaretrulydirecttargetsofmiRNAswillrequirefurtherexperimentalvalidation,buttheenrichmentanalysisshowedthat themRNAswereinvolvedinmeaningfulbiologicalprocessesrelatedtoLungOrganogenesis,e.g.,lungcellproliferation,lung branchingandalveolardevelopment. Thus,inthiscontext,DynOmicshasthepotentialtoidentifynoveltargetsofmiRNAsto aidintherapeuticdevelopment. OurstudyemphasisestheimportanceofjointlystudyingmiRNAandmRNAexpressionto understandthemechanismsofmiRNAregulation. Modelorganismspresentasimplerandmoreconvenientalternativetodirectlystudydiseaseinhumans. Inthemammalian pre-implantation embryonic development study, we showed that DynOmics could identify delayed conserved expression betweendifferentorganisms. Thisisachallengingtask,astimingdifferencesofexpressionchangescanoccurbothinmetabolic processesandacrossorgansfordifferentorganisms.37 Bycorrectingforthesetimingdifferences,DynOmicscanthereforehelp toinfergenefunctionsacrossorganisms,andtherebyintegrateinformationinwholebiologicalprocesses. Suchintegrationmay inturnidentifywhichorganismsprovidesuitablemodelsforhumandiseaseanddrugdiscoveryduetotheconservationin processes.64 Currently, DynOmics has been used to identify associations between datasets of moderate size (∼100 references and ∼10,000queries). Thecomputationaltimewouldincreaseforlargedatasets(∼10,000referencesandqueries). Onesolution couldbetoclusterprofilespriortoapplyingDynOmics,toidentifyspecificpatternsofinterestovertimeasqueriesand/or references. Asthealgorithmisbasedonindependentpairwisecomparisons,parallelcomputingcouldalsobeusedtodecrease thecomputationalburden. Alternatively,asshownintheLungOrganogenesisstudy,theDynOmicsanalysiscanbeperformed onasmallernumberofqueriesselectedbasedonpriorknowledgeorbiologicalassumptions. Conclusion Delays in molecular expression are an acknowledged and important phenomenon in many areas of biology. Here we demonstrated the need for and value of methods that are robust to delays, by showcasing the benefit of accurate delay estimatestointerpretresponsedynamicsandidentifyconservedmolecularmechanisms. DynOmicsovercomesthechallengeof 9/12 bioRxiv preprint first posted online Sep. 20, 2016; doi: http://dx.doi.org/10.1101/076257. The copyright holder for this preprint (which was not peer-reviewed) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is made available under a CC-BY-NC-ND 4.0 International license. integratingdatawithtimingdifferencesofexpressionchangesandthereforepresentsaneffectivetooltostudytime-sensitive molecularexpression. Theintegrationofmultipletimecourse‘omics’dataisbecomingnecessaryinordertounderstand abiologicalsystem’sformation,actionsandregulationwithhighconfidence. OuralgorithmDynOmicsprovidesaunique opportunity to study molecular interactions between multiple functional levels of a single system or multiple organisms. DynOmicsisimplementedintheopensourceprogramminglanguageRandisfreelyavailableviabitbucket. Acknowledgements This work was supported by the Wound Management Innovation CRC (established and supported under the Australian Government’s Cooperative Research Centres Program) [JS], the Australian Cancer Research Foundation (ACRF) for the DiamantinaIndividualisedOncologyCareCentreandNationalHealthandMedicalResearchCouncilCareerDevelopment (NHMRC)fellowship[APP1087415toKALC]andtheAustralianResearchCouncil[DE120101127toBEH]. Contributions JSdevelopedthemethodologies,implementedtheapproaches,performedthestatisticalanalysesandwrotethemanuscript. KALCandBEHhelpedwritingthemanuscript. Allauthorsparticipatedinthedesignofthestudyandreviewedthemanuscript. Competing financial interests Theauthorsdeclarenocompetingfinancialinterests. References 1. Ritchie,M.D.,Holzinger,E.R.,Li,R.,Pendergrass,S.A.&Kim,D. Methodsofintegratingdatatouncovergenotype- phenotypeinteractions. Nat.Rev.Genetics16,85–97(2015). 2. Storey, J.D., Xiao, W., Leek, J.T., Tompkins, R.G.&Davis, R.W. Significanceanalysisoftimecoursemicroarray experiments. PNAS102,12837–42(2005). 3. Straube,J.,Gorse,A.-D.,Huang,B.E.&LeˆCao,K.-A. Alinearmixedmodelsplineframeworkforanalyzingtimecourse ‘omics’data. PLOSONE10,e0134540(2015b). 4. Tai,Y.C.,Speed,T.P.etal. Amultivariateempiricalbayesstatisticforreplicatedmicroarraytimecoursedata. TheAnnals ofStatistics34,2387–2412(2006). 5. Aryee,M.J.,Gutie´rrez-Pabello,J.A.,Kramnik,I.,Maiti,T.&Quackenbush,J. Animprovedempiricalbayesapproachto estimatingdifferentialgeneexpressioninmicroarraytime-coursedata: Betr(bayesianestimationoftemporalregulation). BMCbioinformatics10,409(2009). 6. Stegle,O.etal. Arobustbayesiantwo-sampletestfordetectingintervalsofdifferentialgeneexpressioninmicroarray timeseries. J.Comp.Biol17,355–367(2010). 7. Leng,N.etal. Ebseq-hmm: abayesianapproachforidentifyinggene-expressionchangesinorderedrna-seqexperiments. Bioinformaticsbtv193(2015). 8. Kalaitzis,A.A.&Lawrence,N.D. Asimpleapproachtorankingdifferentiallyexpressedgeneexpressiontimecourses throughgaussianprocessregression. BMCbioinformatics12,1(2011). 9. Heinonen,M.etal. Detectingtimeperiodsofdifferentialgeneexpressionusinggaussianprocesses: anapplicationto endothelialcellsexposedtoradiotherapydosefraction. Bioinformaticsbtu699(2014). 10. A¨ijo¨, T. et al. Methods for time series analysis of rna-seq data with application to human th17 cell differentiation. Bioinformatics30,i113–i120(2014). 11. Conesa,A.,Nueda,M.J.,Ferrer,A.&Talo´n,M. masigpro: amethodtoidentifysignificantlydifferentialexpression profilesintime-coursemicroarrayexperiments. Bioinformatics22,1096–1102(2006). 12. De´jean,S.,Martin,P.G.,Baccini,A.&Besse,P. Clusteringtime-seriesgeneexpressiondatausingsmoothingspline derivatives. EURASIPJBioinformSystBiol2007,1–10(2007). 13. Luan,Y.&Li,H.Clusteringoftime-coursegeneexpressiondatausingamixed-effectsmodelwithb-splines.Bioinformatics 19,474–482(2003). 14. Ernst,J.,Nau,G.J.&Bar-Joseph,Z. Clusteringshorttimeseriesgeneexpressiondata. Bioinformatics21,i159–i168 (2005). 10/12