DynOmics to identify delays and co-expression patterns across time course experiments

12 Pages·2016·0.47 MB·English

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
*[e-mail]
+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