Mercurial > repos > iuc > raceid_inspectclusters

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