/BrownianInsertor.R

http://github.com/sbotond/phylosim · R · 660 lines · 230 code · 62 blank · 368 comment · 33 complexity · 652347048cfdbbb71074e00555530cbf MD5 · raw file

  1. ## Copyright 2009 Botond Sipos
  2. ## See the package description for licensing information.
  3. ##
  4. ##########################################################################/**
  5. #
  6. # @RdocClass BrownianInsertor
  7. #
  8. # @title "The BrownianInsertor class"
  9. #
  10. # \description{
  11. # The \code{BrownianInsertor} class inherits from the \code{DiscreteInsertor}
  12. # or \code{ContinuousInsertor} class depending on the \code{type} constructor argument
  13. # ("discrete" or "continuous").
  14. #
  15. # This process generates the insert sequence based on the sites flanking the insertions as follows:
  16. # \itemize{
  17. # \item An insert length is sampled by calling the function stored in the \code{proposeBy} virtual field.
  18. # \item A sequence object is constructed.
  19. # \item The processes attached to both flanking sites are attached to the insert sequence. If there are no common processes, the processes from a randomly chosen site will be attached to the insert.
  20. # \item The site-process specific parameters are sampled from Brownian paths with linear trends having the values from the flanking sites as endpoints.
  21. # }
  22. #
  23. # The "noisiness" of the Brownian path can be controlled through the \code{scale} virtual field/constructor parameter.
  24. #
  25. # @classhierarchy
  26. # }
  27. #
  28. # @synopsis
  29. #
  30. # \arguments{
  31. # \item{name}{Object name.}
  32. # \item{type}{Process type (see above).}
  33. # \item{scale}{Brownian path scale parameter.}
  34. # \item{...}{Additional arguments.}
  35. # }
  36. #
  37. # \section{Fields and Methods}{
  38. # @allmethods
  39. # }
  40. #
  41. # \examples{
  42. # # create a BrownianInsertor process, discrete type
  43. # p<-BrownianInsertor(
  44. # type="discrete",
  45. # scale=0.05,
  46. # sizes=1:4,
  47. # probs=c(3/6,1/6,1/6,1/6),
  48. # rate=0.05
  49. # )
  50. # # get object summary
  51. # summary(p)
  52. # # plot insert length distribution
  53. # plot(p)
  54. # # create a nucleotide sequence, attach processes
  55. # s<-NucleotideSequence(string="AAAAAAAAAAA",processes=list(list(p,JC69())))
  56. # # create simulation object
  57. # sim<-PhyloSim(root.seq=s, phylo=rcoal(2))
  58. # # simulate and show alignment
  59. # Simulate(sim)
  60. # sim$alignment
  61. # # check the rate multipliers and insertion tolerances in one of the sequences
  62. # res<-sim$sequences[[2]]
  63. # getRateMultipliers(res,p)
  64. # getInsertionTolerance(res,p)
  65. # }
  66. #
  67. # @author
  68. #
  69. # \seealso{
  70. # DiscreteInsertor ContinuousInsertor GeneralInsertor GeneralInDel
  71. # }
  72. #
  73. #*/###########################################################################
  74. setConstructorS3(
  75. "BrownianInsertor",
  76. function(
  77. name="Anonymous",
  78. type="discrete",
  79. scale=0.001,
  80. ...
  81. ) {
  82. if(type == "continuous"){
  83. this<-ContinuousInsertor(
  84. ...
  85. );
  86. }
  87. else if (type == "discrete") {
  88. this<-DiscreteInsertor(
  89. ...
  90. );
  91. }
  92. else {
  93. throw("Invalid insertor process type!\n");
  94. }
  95. this<-extend(
  96. this,
  97. "BrownianInsertor",
  98. .scale = NA,
  99. .type = type
  100. );
  101. # Using virtual field to clear Id cache:
  102. this$name<-name;
  103. # setting scale
  104. this$scale<-scale;
  105. this$generateBy<-function(process=NA,length=NA,target.seq=NA,event.pos=NA,insert.pos=NA){
  106. if(is.na(target.seq)){
  107. return(NA);
  108. }
  109. if(is.na(length) | (length(length) == 0) | length == 0){
  110. throw("Invalid insert length!\n");
  111. }
  112. # The start and end of the Brownian path:
  113. start;
  114. end;
  115. proc<-list();
  116. if( (event.pos == 1) || (event.pos == target.seq$.length) ){
  117. start<-clone(target.seq$.sites[[event.pos]]);
  118. start$.state=NA;
  119. end<-clone(start);
  120. proc<-getProcesses(start);
  121. } else {
  122. start<-clone(target.seq$.sites[[insert.pos]]);
  123. start$.state=NA;
  124. end<-clone(target.seq$.sites[[insert.pos + 1]]);
  125. end$.state=NA;
  126. proc.start<-getProcesses(start);
  127. proc.end<-getProcesses(end);
  128. # Calculate the intersection of process list:
  129. proc<-PSRoot$intersect.list(proc.start,proc.end);
  130. # No common processes:
  131. if(length(proc) == 0){
  132. coin.flip<-sample(c(0,1),1);
  133. if (coin.flip) {
  134. proc <- getProcesses(start)
  135. end <- clone(start)
  136. } else {
  137. proc <- getProcesses(end)
  138. start <- clone(end)
  139. }
  140. }
  141. }
  142. # Create the insert sequence:
  143. class.seq<-class(target.seq)[[1]];
  144. insert<-do.call(class.seq,list(length=length));
  145. setProcesses(this=insert,value=list(proc));
  146. # For every process...
  147. for (p in proc){
  148. # ... and site specific parameter:
  149. for(param in getSiteSpecificParamIds(p)){
  150. start.value<-getParameterAtSite(p,site=start,id=param)$value;
  151. end.value<-getParameterAtSite(p,site=end,id=param)$value;
  152. path<-seq(from=start.value,to=end.value,length.out=(insert$.length + 2));
  153. path<-path[2:(length(path)-1)];
  154. brownian.path<-abs(BrownianInsertor$BrownianPath(p=path, a=this$.scale));
  155. # Tolerance values are probabilities, do not alow them to wander beyond 1:
  156. if(param == "insertion.tolerance" || param == "deletion.tolerance"){
  157. brownian.path[which(brownian.path > 1)]<-1;
  158. }
  159. setParameterAtSites(
  160. insert,
  161. process = p,
  162. id = param,
  163. value = brownian.path
  164. );
  165. }
  166. }
  167. return(insert);
  168. }
  169. return(this);
  170. },
  171. enforceRCC=TRUE
  172. );
  173. ##
  174. ## Method: checkConsistency
  175. ##
  176. ###########################################################################/**
  177. #
  178. # @RdocMethod checkConsistency
  179. #
  180. # @title "Check object consistency"
  181. #
  182. # \description{
  183. # @get "title".
  184. # }
  185. #
  186. # @synopsis
  187. #
  188. # \arguments{
  189. # \item{this}{An object.}
  190. # \item{...}{Not used.}
  191. # }
  192. #
  193. #
  194. # \value{
  195. # Returns an invisible TRUE if no inconsistencies found in the object, throws
  196. # an error otherwise.
  197. # }
  198. #
  199. # @author
  200. #
  201. # \seealso{
  202. # @seeclass
  203. # }
  204. #
  205. #*/###########################################################################
  206. setMethodS3(
  207. "checkConsistency",
  208. class="BrownianInsertor",
  209. function(
  210. this,
  211. ...
  212. ){
  213. wp<-this$writeProtected;
  214. if (wp) {
  215. this$writeProtected<-FALSE;
  216. }
  217. may.fail<-function(this) {
  218. this$scale<-this$scale;
  219. if( (this$.type != "discrete") && (this$.type != "continuous") ){
  220. throw("BrownianInsertor type is invalid!\n");
  221. }
  222. }
  223. tryCatch(may.fail(this),finally=this$writeProtected<-wp);
  224. NextMethod();
  225. },
  226. private=FALSE,
  227. protected=FALSE,
  228. overwrite=FALSE,
  229. conflict="warning",
  230. validators=getOption("R.methodsS3:validators:setMethodS3")
  231. );
  232. ##
  233. ## Method: BrownianPath
  234. ##
  235. ###########################################################################/**
  236. #
  237. # @RdocMethod BrownianPath
  238. #
  239. # @title "Generate a Brownian path"
  240. #
  241. # \description{
  242. # @get "title".
  243. #
  244. # This method generates a Brownian path given the scale parameter a (determining "noisiness")
  245. # and the vector p describing the trends. More useful as a static method.
  246. # }
  247. #
  248. # @synopsis
  249. #
  250. # \arguments{
  251. # \item{this}{A BrownianInsertor object.}
  252. # \item{p}{Path parameter (a numeric vector).}
  253. # \item{a}{Scale paramater (a numeric vector of length one).}
  254. # \item{...}{Not used.}
  255. # }
  256. #
  257. # \value{
  258. # A numeric vector.
  259. # }
  260. #
  261. # \examples{
  262. # path<-BrownianInsertor$BrownianPath(a=2,p=1:10);
  263. # path
  264. # plot(path)
  265. # }
  266. #
  267. # @author
  268. #
  269. # \seealso{
  270. # @seeclass
  271. # }
  272. #
  273. #*/###########################################################################
  274. setMethodS3(
  275. "BrownianPath",
  276. class="BrownianInsertor",
  277. function(
  278. this,
  279. p=NA,
  280. a=NA,
  281. ...
  282. ){
  283. generate_brownian<-function(length, a){
  284. cumsum(rnorm(length,0,sd=a));
  285. }
  286. generate_bridge <- function (length,a){
  287. b <- generate_brownian(length,a)
  288. b - (1:length)/(length) * b[length]
  289. }
  290. generate_path <- function (p,a){
  291. n <- length(p);
  292. b <- generate_bridge (n+1,a);
  293. p + b[1:n];
  294. }
  295. return(generate_path(p,a));
  296. },
  297. private=FALSE,
  298. protected=FALSE,
  299. overwrite=FALSE,
  300. conflict="warning"
  301. );
  302. ##
  303. ## Method: getType
  304. ##
  305. ###########################################################################/**
  306. #
  307. # @RdocMethod getType
  308. #
  309. # @title "Get the type of the BrownianInsertor process"
  310. #
  311. # \description{
  312. # @get "title".
  313. #
  314. # If type is \code{discrete}, than the process inherits from \code{DiscreteInsertor}.
  315. # If the type is \code{continuous}, then the object inherits from \code{ContinuousInsertor}.
  316. # }
  317. #
  318. # @synopsis
  319. #
  320. # \arguments{
  321. # \item{this}{A BrownianInsertor object.}
  322. # \item{...}{Not used.}
  323. # }
  324. #
  325. # \value{
  326. # A character vector of length one.
  327. # }
  328. #
  329. # \examples{
  330. # p<-BrownianInsertor(type="discrete")
  331. # # get type
  332. # getType(p)
  333. # # get upstream classes
  334. # class(p)
  335. # p<-BrownianInsertor(type="continuous")
  336. # # get type via virtual field
  337. # p$type
  338. # # get upstream classes
  339. # class(p)
  340. # }
  341. #
  342. # @author
  343. #
  344. # \seealso{
  345. # @seeclass
  346. # }
  347. #
  348. #*/###########################################################################
  349. setMethodS3(
  350. "getType",
  351. class="BrownianInsertor",
  352. function(
  353. this,
  354. ...
  355. ){
  356. this$.type;
  357. },
  358. private=FALSE,
  359. protected=FALSE,
  360. overwrite=FALSE,
  361. conflict="warning",
  362. validators=getOption("R.methodsS3:validators:setMethodS3")
  363. );
  364. ##
  365. ## Method: setType
  366. ##
  367. ###########################################################################/**
  368. #
  369. # @RdocMethod setType
  370. #
  371. # @title "Forbidden action: setting the type of a BrownianInsertor object"
  372. #
  373. # \description{
  374. # @get "title".
  375. #
  376. # The type can be set only from the \code{type} constructor argument and cannot be changed later.
  377. # }
  378. #
  379. # @synopsis
  380. #
  381. # \arguments{
  382. # \item{this}{An object.}
  383. # \item{value}{Not used.}
  384. # \item{...}{Not used.}
  385. # }
  386. #
  387. # \value{
  388. # Throws an error.
  389. # }
  390. #
  391. # @author
  392. #
  393. # \seealso{
  394. # @seeclass
  395. # }
  396. #
  397. #*/###########################################################################
  398. setMethodS3(
  399. "setType",
  400. class="BrownianInsertor",
  401. function(
  402. this,
  403. value,
  404. ...
  405. ){
  406. throw("The type of the BrownianInsertor objects can be set only from the constructor!\n");
  407. },
  408. private=FALSE,
  409. protected=FALSE,
  410. overwrite=FALSE,
  411. conflict="warning",
  412. validators=getOption("R.methodsS3:validators:setMethodS3")
  413. );
  414. ##
  415. ## Method: getScale
  416. ##
  417. ###########################################################################/**
  418. #
  419. # @RdocMethod getScale
  420. #
  421. # @title "Get scale parameter"
  422. #
  423. # \description{
  424. # @get "title".
  425. # }
  426. #
  427. # @synopsis
  428. #
  429. # \arguments{
  430. # \item{this}{A BrownianInsertor object.}
  431. # \item{...}{Not used.}
  432. # }
  433. #
  434. # \value{
  435. # A numeric vector of length one.
  436. # }
  437. #
  438. # \examples{
  439. # # create a BrownianInsertor object
  440. # p<-BrownianInsertor(scale=0.002)
  441. # # set/get scale parameter
  442. # setScale(p,0.1)
  443. # getScale(p)
  444. # # set/get scale parameter via virtual field
  445. # p$scale<-0.1
  446. # p$scale
  447. # }
  448. #
  449. # @author
  450. #
  451. # \seealso{
  452. # @seeclass
  453. # }
  454. #
  455. #*/###########################################################################
  456. setMethodS3(
  457. "getScale",
  458. class="BrownianInsertor",
  459. function(
  460. this,
  461. ...
  462. ){
  463. this$.scale;
  464. },
  465. private=FALSE,
  466. protected=FALSE,
  467. overwrite=FALSE,
  468. conflict="warning",
  469. validators=getOption("R.methodsS3:validators:setMethodS3")
  470. );
  471. ##
  472. ## Method: setScale
  473. ##
  474. ###########################################################################/**
  475. #
  476. # @RdocMethod setScale
  477. #
  478. # @title "Set scale parameter"
  479. #
  480. # \description{
  481. # @get "title".
  482. # }
  483. #
  484. # @synopsis
  485. #
  486. # \arguments{
  487. # \item{this}{A BrownianInsertor object.}
  488. # \item{value}{A numeric vector of length one.}
  489. # \item{...}{Not used.}
  490. # }
  491. #
  492. # \value{
  493. # value (invisible).
  494. # }
  495. #
  496. # \examples{
  497. # # create a BrownianInsertor object
  498. # p<-BrownianInsertor(scale=0.002)
  499. # # set/get scale parameter
  500. # setScale(p,0.1)
  501. # getScale(p)
  502. # # set/get scale parameter via virtual field
  503. # p$scale<-0.1
  504. # p$scale
  505. # }
  506. #
  507. # @author
  508. #
  509. # \seealso{
  510. # @seeclass
  511. # }
  512. #
  513. #*/###########################################################################
  514. setMethodS3(
  515. "setScale",
  516. class="BrownianInsertor",
  517. function(
  518. this,
  519. value,
  520. ...
  521. ){
  522. .checkWriteProtection(this);
  523. if(!is.numeric(value) || (length(value) != 1)){
  524. throw("The value of the scale paramter must be a numeric vector of length 1!\n");
  525. }
  526. this$.scale<-value;
  527. return(invisible(this$.scale));
  528. },
  529. private=FALSE,
  530. protected=FALSE,
  531. overwrite=FALSE,
  532. conflict="warning",
  533. validators=getOption("R.methodsS3:validators:setMethodS3")
  534. );
  535. ##
  536. ## Method: summary
  537. ##
  538. ###########################################################################/**
  539. #
  540. # @RdocMethod summary
  541. #
  542. # @title "Summarize the properties of an object"
  543. #
  544. # \description{
  545. # @get "title".
  546. # }
  547. #
  548. # @synopsis
  549. #
  550. # \arguments{
  551. # \item{object}{An object}
  552. # \item{...}{Not used.}
  553. # }
  554. #
  555. # \value{
  556. # Returns a PSRootSummary object.
  557. # }
  558. #
  559. # \examples{
  560. #
  561. # # create an object
  562. # p<-BrownianInsertor(
  563. # type="discrete",
  564. # scale=0.05,
  565. # sizes=1:4,
  566. # probs=c(3/6,1/6,1/6,1/6),
  567. # rate=0.05
  568. # )
  569. # # get a summary
  570. # summary(p)
  571. # }
  572. #
  573. # @author
  574. #
  575. # \seealso{
  576. # @seeclass
  577. # }
  578. #
  579. #*/###########################################################################
  580. setMethodS3(
  581. "summary",
  582. class="BrownianInsertor",
  583. function(
  584. object,
  585. ...
  586. ){
  587. this<-object;
  588. .addSummaryNameId(this);
  589. this$.summary$"Type"<-this$.type;
  590. this$.summary$"Brownian path scale parameter"<-this$.scale;
  591. NextMethod();
  592. },
  593. private=FALSE,
  594. protected=FALSE,
  595. overwrite=FALSE,
  596. conflict="warning",
  597. validators=getOption("R.methodsS3:validators:setMethodS3")
  598. );