Mercurial > repos > iuc > raceid_clustering

comparison scripts/trajectoryinspect.R @ 0:4ea021bd7513 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/raceid3 commit f880060c478d42202df5b78a81329f8af56b1138
author iuc
date Thu, 22 Nov 2018 04:43:57 -0500
parents
children a4b734cd253b
comparison
equal deleted inserted replaced
-1:000000000000 0:4ea021bd7513
1 #!/usr/bin/env R
2 VERSION = "0.2"
3
4 args = commandArgs(trailingOnly = T)
5
6 if (length(args) != 1){
7 message(paste("VERSION:", VERSION))
8 stop("Please provide the config file")
9 }
10
11 suppressWarnings(suppressPackageStartupMessages(require(RaceID)))
12 suppressWarnings(suppressPackageStartupMessages(require(FateID)))
13 source(args[1])
14
15 test <- list()
16 test$side = 3
17 test$line = 2.5
18 second <- test
19 second$cex = 0.5
20 second$line = 2.5
21
22 do.trajectoryinspection.stemID <- function(ltr){
23 makeBranchLink <- function(i,j,k){
24 ingoing <- paste(sort(c(i,j)), collapse=".")
25 outgoing <- paste(sort(c(j,k)), collapse=".")
26 messed <- sort(c(ingoing,outgoing))
27 return(list(messed[[1]], messed[[2]]))
28 }
29
30 zzz <- do.call(getproj, c(ltr, trjsid.getproj))
31 bra <- branchcells(
32 ltr,
33 do.call("makeBranchLink", as.list(trjsid.branchcells.ijk))
34 )
35 write.table(
36 head(bra$diffgenes$z, trjsid.numdiffgenes),
37 file=out.diffgenes)
38
39 par(mfrow = c(2,2), cex=0.5)
40 print(do.call(plotmap, c(bra$scl, final=FALSE, fr=FALSE)))
41 print(do.call(mtext, c("Initial Clusters (tSNE)", test)))
42 print(do.call(plotmap, c(bra$scl, final=TRUE, fr=FALSE)))
43 print(do.call(mtext, c("Final Clusters (tSNE)", test)))
44 print(do.call(plotmap, c(bra$scl, final=FALSE, fr=TRUE)))
45 print(do.call(mtext, c("Initial Clusters (F-R)", test)))
46 print(do.call(plotmap, c(bra$scl, final=TRUE, fr=TRUE)))
47 print(do.call(mtext, c("Final Clusters (F-R)", test)))
48 }
49
50 do.trajectoryinspection.fateID <- function(ltr){
51 n <- do.call(cellsfromtree, c(ltr, trjfid.cellsfrom))
52 x <- getfdata(ltr@sc)
53
54 trjfid.filterset$x = x
55 trjfid.filterset$n = n$f
56 fs <- do.call(filterset, c(trjfid.filterset))
57 trjfid.getsom$x = fs
58 s1d <- do.call(getsom, c(trjfid.getsom))
59 trjfid.procsom$s1d = s1d
60 ps <- do.call(procsom, c(trjfid.procsom))
61
62 y <- ltr@sc@cpart[n$f]
63 fcol <- ltr@sc@fcol
64
65 trjfid.plotheat$xpart = y
66 trjfid.plotheat$xcol = fcol
67
68 ##Plot average z-score for all modules derived from the SOM:
69 trjfid.plotheat$x = ps$nodes.z
70 trjfid.plotheat$ypart = unique(ps$nodes)
71 print(do.call(plotheatmap, c(trjfid.plotheat)))
72 print(do.call(mtext, c("Average z-score for all modules derived from SOM", test)))
73 ##Plot z-score profile of each gene ordered by SOM modules:
74 trjfid.plotheat$x = ps$all.z
75 trjfid.plotheat$ypart = ps$nodes
76 print(do.call(plotheatmap, c(trjfid.plotheat)))
77 print(do.call(mtext, c("z-score profile of each gene ordered by SOM modules", test)))
78 ##Plot normalized expression profile of each gene ordered by SOM modules:
79 trjfid.plotheat$x = ps$all.e
80 trjfid.plotheat$ypart = ps$nodes
81 print(do.call(plotheatmap, c(trjfid.plotheat)))
82 print(do.call(mtext, c("Normalized expression profile of each gene ordered by SOM modules", test)))
83 ##Plot binarized expression profile of each gene (z-score < -1, -1 < z-score < 1, z-score > 1):
84 trjfid.plotheat$x = ps$all.b
85 trjfid.plotheat$ypart = ps$nodes
86 print(do.call(plotheatmap, c(trjfid.plotheat)))
87 print(do.call(mtext, c("Binarized expression profile of each gene", test)))
88 ## This should be written out, and passed back into the tool
89 ## to perform sominspect
90 return(list(fs=fs,ps=ps,y=y,fcol=fcol,nf=n$f))
91 }
92
93 do.trajectoryinspection.fateID.sominspect <- function(domo){
94 g <- trjfidsomi.use.genes
95 if (class(g) == "numeric"){
96 g <- names(ps$nodes)[ps$nodes %in% g]
97 }
98
99 typ = NULL
100 if (!is.null(trjfidsomi.use.types)){
101 typ = sub(trjfidsomi.use.types,"", domo$nf)
102 }
103
104 trjfidsomi$x = domo$fs
105 trjfidsomi$y = domo$y
106 trjfidsomi$g = g
107 trjfidsomi$n = domo$nf
108 trjfidsomi$col = domo$fcol
109 trjfidsomi$types = typ
110
111 ## The average pseudo-temporal expression profile of this group
112 ## can be plotted by the function plotexpression:
113 par(mfrow = c(1,1))
114 test$cex = 1
115 second$line = 1.5
116 if (trjfidsomi$name == "Title") trjfidsomi$name = ""
117 print(do.call(plotexpression, c(trjfidsomi)))
118 mess2 <- paste(c(trjfidsomi.use.genes), collapse=", ")
119 mess1 <- "Average pseudo-temporal expression profile"
120 print(do.call(mtext, c(mess1, test)))
121 print(do.call(mtext, c(mess2, second)))
122 }
123
124 ltr <- in.rdat
125
126 pdf(out.pdf)
127 if (perform.stemID) do.trajectoryinspection.stemID(ltr)
128 if (perform.fateID) {
129 domo <- do.trajectoryinspection.fateID(ltr)
130 if (perform.fateID.sominspect) {
131 do.trajectoryinspection.fateID.sominspect(domo)
132 }
133 }
134 dev.off()