-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy path05-lincs.Rmd
150 lines (121 loc) · 4.65 KB
/
05-lincs.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
# LINCS
This is LINCS
## Signatures
Putative drug therapies were identified with signature-based connectivity analysis utilizing the Library of Integrated Network-based Signatures (LINCS) database.
```{r, results='asis'}
#Input: Annotated DrugFindr investigate signature results
generate_moa_report <- function(X) {
X %>%
select(integratedMoas, Target, GeneTargets) %>%
filter(integratedMoas != "" & !is.na(integratedMoas)) %>%
separate_rows(integratedMoas, sep = "\\|") %>%
separate_rows(GeneTargets, sep = "\\|") %>%
mutate(across(c(integratedMoas, GeneTargets), str_trim)) %>%
group_by(integratedMoas) %>%
summarise(
Target = paste(unique(Target), collapse = "|"),
GeneTargets = paste(unique(GeneTargets), collapse = "|"),
N = n()
) %>%
arrange(desc(N))
}
#Input: Annotated DrugFindr investigate signature results | DEG data
generate_gene_report <- function(X, data) {
X %>%
select(Symbol = GeneTargets) %>%
filter(Symbol != "" & !is.na(Symbol)) %>%
separate_rows(Symbol, sep = "\\|") %>%
mutate(Symbol = str_trim(Symbol)) %>%
count(Symbol) %>%
arrange(desc(n)) %>%
inner_join(global_state$hgnc, by = "Symbol") %>%
inner_join(data, by = "Symbol")
}
#Input: DEG signature
#Output: Named list of relevant LINCS results
get_ilincs_results <- function(X) {
if(global_state$species != "human") { #convert mouse or rat orthologs to human genes for lincs
orthologs <- orthologs(X$Symbol, human = F, species = global_state$species) %>%
select(Symbol = human_symbol, Symbol_old = symbol)
X <- X %>%
rename(Symbol_old = Symbol) %>%
inner_join(orthologs, by = "Symbol_old") %>%
select(-Symbol_old)
}
signatures <- investigateSignature(
X,
filterThreshold = 0,
outputLib = "CP",
geneColumn = "Symbol",
logfcColumn = "log2FoldChange",
pvalColumn = "pvalue"
)
concordant <- signatures %>%
filter(Similarity >= 0) %>% inner_join(get_ilincs_metadata(.$TargetSignature), by = "TargetSignature") %>%
left_join(global_state$lincs_fda, by = c("Target" = "sm_name"))
discordant <- signatures %>%
filter(Similarity < 0) %>% inner_join(get_ilincs_metadata(.$TargetSignature), by = "TargetSignature") %>%
left_join(global_state$lincs_fda, by = c("Target" = "sm_name"))
concordant_moa_report <- generate_moa_report(concordant)
discordant_moa_report <- generate_moa_report(discordant)
concordant_gene_report <- generate_gene_report(concordant, X)
discordant_gene_report <- generate_gene_report(discordant, X)
discordant_pathways <- do_enrichr(discordant_gene_report$Symbol)
system("sleep 3")
concordant_pathways <- do_enrichr(concordant_gene_report$Symbol)
lst(concordant,
discordant,
concordant_moa_report,
discordant_moa_report,
concordant_gene_report,
discordant_gene_report,
discordant_pathways,
concordant_pathways)
}
global_state$data %<>%
map(~ update_list(., results = c(.$results, lincs = list(get_ilincs_results(.$data)))))
global_state$data %>%
map(~ knit_child("Rmd/lincs_signatures.Rmd", envir = environment(), quiet = TRUE)) %>%
list_c() %>%
cat(sep = "\n")
```
## MOAS and Genes
These are LINCS Perturbagens MoAs and Genetargets
```{r, results='asis'}
global_state$data %>%
map(~ knit_child("Rmd/lincs_metadata.Rmd", envir = environment(), quiet = TRUE)) %>%
list_c() %>%
cat(sep = "\n")
```
## Pathways
These are LINCS Pathways
```{r, results='asis'}
global_state$data %>%
map(~ knit_child("Rmd/lincs_pathways.Rmd", envir = environment(), quiet = TRUE)) %>%
list_c() %>%
cat(sep = "\n")
```
## PAVER
These are PAVER plots of the pathways from all comparisons
```{r, results='asis'}
input = global_state$data %>%
map(~ bind_rows(.$results$lincs$concordant_pathways, .$results$lincs$discordant_pathways %>% mutate(Combined.Score = Combined.Score * -1))) %>%
bind_rows(.id = "Group") %>%
select(GOID, CS = Combined.Score, Group) %>%
mutate(CS = sign(CS) * log2(abs(CS)) + 1) %>%
mutate(Type = ifelse(CS > 0, "C", "D")) %>%
pivot_wider(names_from = c(Group, Type), values_from = CS, names_sep = "-")
minClusterSize = 5
maxCoreScatter = 0.33
minGap = (1 - maxCoreScatter) * 3 / 4
LINCS_PAVER_result <- quiet(PAVER::prepare_data(input, global_state$embeddings, global_state$term2name) %>%
PAVER::generate_themes(maxCoreScatter = maxCoreScatter,
minGap = minGap,
minClusterSize = minClusterSize))
global_state$results <- c(global_state$results, lst(LINCS_PAVER_result))
list(LINCS_PAVER_result, "logCS") %>%
knit_child(text = readLines("Rmd/paver_report.Rmd"),
envir = environment(),
quiet = TRUE) %>%
cat(sep = "\n")
```