-
Notifications
You must be signed in to change notification settings - Fork 0
/
FigS8_CIBERSORTx
150 lines (126 loc) · 6.08 KB
/
FigS8_CIBERSORTx
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
#STEP1 Prepare Reference data with snRNA-sea
library(Seurat)
library(ggplot2)
library(dplyr)
library(biomaRt)
library(Matrix)
set.seed(1234)
rnaAggr <- readRDS("~/ATAC_IRI_project/IRI_RNA.rds")
Idents(rnaAggr) <- "name"
PTCrna <- subset(rnaAggr,idents = c("PTS1","PTS2","PTS3","NewPT1","NewPT2"))
PTCrna <- subset(PTCrna,downsample = 3816) #down size the data (3816 cells for each of PTS1/2/3,NewPT1/2 = a total of 19080 cells)
cluster.ids <- c("NormPTC","NormPTC","NormPTC","InjuredPTC","FRPTC") #"PTS1/2/3 = "NormPTC", NewPT1= "InjuredPTC", NewPT2 = "FRPTC
names(cluster.ids) <- levels(PTCrna)
PTCrna <- RenameIdents(PTCrna, cluster.ids)
PTCrna@meta.data[["celltype"]] <- PTCrna@active.ident
Idents(PTCrna) <- "celltype"
mat <- PTCrna@assays[["RNA"]]@counts
mat <- as.data.frame(mat) #27133x19080
#convert mouse genes to human genes
human = useMart("ensembl", dataset = "hsapiens_gene_ensembl", host = "https://dec2021.archive.ensembl.org/")
mouse = useMart("ensembl", dataset = "mmusculus_gene_ensembl", host = "https://dec2021.archive.ensembl.org/")
mouse_gene <- rownames(PTCrna)
converted_genelist = getLDS(attributes = "mgi_symbol", filters = "mgi_symbol", values = mouse_gene ,
mart = mouse, attributesL = "hgnc_symbol", martL = human, uniqueRows=T)
converted_genelist <- converted_genelist %>% distinct(HGNC.symbol,.keep_all=TRUE)
mat$MGI.symbol <- rownames(mat)
mat <- merge(mat,converted_genelist,by = "MGI.symbol")
rownames(mat)<- mat$HGNC.symbol
mat <- mat[,2:19081]
meta <- PTCrna@meta.data
meta <- t(meta)
meta <- meta[20:21,]
mat <- rbind(meta,mat)
mat <- mat[-1,]
rownames(mat)[1] <- "GeneSymbol"
write.table(mat,"IRI_mat2.txt",quote = F,col.names = F,sep="\t") #Reference file
#STEP2
#RUN on Linux(#docker run)
#cibersortx/fractions --username xxxxx --token xxxxxxxxxxxxxx --single_cell TRUE --mixture RPTEC_RNAseq.txt --refsample IRI_mat2.txt --rmbatchSmode TRUE --perm 100
Deconv <- read.table("path/to/out/CIBERSORTx_Adjusted.txt",header = T)
library(reshape2)
lm = melt(Deconv[,1:4], id = "Mixture")
ggplotColours <- function(n = 6, h = c(0, 360) + 15){
if ((diff(h) %% 360) < 1) h[2] <- h[2] - 360/n
hcl(h = (seq(h[1], h[2], length = n)), c = 100, l = 65)
}
colours <- ggplotColours(n=3)
FigS8B <- ggplot(lm, aes(x = Mixture, fill = variable, y = value)) +
geom_bar(stat = "identity", colour = "black") +
scale_y_continuous(limits = c(0, 1),expand = c(0, 0)) +theme(
panel.background = element_rect(fill='transparent')) +
scale_fill_manual(values = colours[c(3,2,1)])
lm_cont <- lm[c(13:15,28:30,43:45),] #only controls
mean <- c(mean(lm_cont$value[1:3]),mean(lm_cont$value[4:6]),mean(lm_cont$value[7:9]))
celltype <- c(lm_cont$variable[c(1,4,7)])
test <- data.frame(celltype,mean)
FigS8A <- ggpie(test, "mean", label = "celltype",lab.pos = "in",
fill = "celltype", color = "black",
palette = colours[c(3,2,1)]) #610x470
#normal
summary(lm_cont$value[1:3]) #0.194
sd(lm_cont$value[1:3]) #0.013
#FRPTC
summary(lm_cont$value[4:6]) #0.304
sd(lm_cont$value[4:6]) #0.0090
#Injured
summary(lm_cont$value[7:9]) #0.502
sd(lm_cont$value[7:9]) #0.023
lm2 <- lm[c(43:45,34:36,31:33,40:42,37:39),]
lm2$Mixture2 <- c("siCont","siCont","siCont",
"siCREB5","siCREB5","siCREB5",
"siHIVEP2","siHIVEP2","siHIVEP2",
"siNFAT5","siNFAT5","siNFAT5",
"siPPARA","siPPARA","siPPARA")
FigS8D <- ggbarplot(lm2, x = "Mixture2", y = "value",
add = c("mean_sd", "jitter"),color = "Mixture2",
position = position_dodge(0.8),ylim = c(0.25, 0.75))+scale_y_continuous(expand = c(0, 0)) #400x440
lm3 <- lm[c(28:30,19:21,16:18,25:27,22:24),]
lm3$Mixture2 <- c("siCont","siCont","siCont",
"siCREB5","siCREB5","siCREB5",
"siHIVEP2","siHIVEP2","siHIVEP2",
"siNFAT5","siNFAT5","siNFAT5",
"siPPARA","siPPARA","siPPARA")
FigS8C <- ggbarplot(lm3, x = "Mixture2", y = "value",
add = c("mean_sd", "jitter"),color = "Mixture2",
position = position_dodge(0.8),ylim = c(0, 0.5))+scale_y_continuous(expand = c(0, 0))
#400x440
library(multcomp)
lm2$Mixture2 = as.factor(lm2$Mixture2)
res.aov <- aov(value ~ Mixture2, data =lm2)
summary(res.aov)
#Df Sum Sq Mean Sq F value Pr(>F)
#Mixture 4 0.03490 0.008725 11.86 0.00082 ***
# Summary of the analysis
Injured_turkey <- TukeyHSD(res.aov)
#diff lwr upr p adj
#siCREB5-siCont -0.079288389 -0.152163435 -0.006413342 0.0318915
#siHIVEP2-siCont -0.034037284 -0.106912330 0.038837763 0.5637304
#siNFAT5-siCont 0.066867401 -0.006007646 0.139742447 0.0760774
#siPPARA-siCont 0.006289212 -0.066585835 0.079164258 0.9983208
#siHIVEP2-siCREB5 0.045251105 -0.027623941 0.118126151 0.3135079
#siNFAT5-siCREB5 0.146155790 0.073280743 0.219030836 0.0004512
#siPPARA-siCREB5 0.085577600 0.012702554 0.158452647 0.0205625
#siNFAT5-siHIVEP2 0.100904685 0.028029638 0.173779731 0.0072274
#siPPARA-siHIVEP2 0.040326495 -0.032548551 0.113201542 0.4134346
#siPPARA-siNFAT5 -0.060578189 -0.133453236 0.012296857 0.1173002
library(multcomp)
lm3$Mixture2 = as.factor(lm3$Mixture2)
res.aov <- aov(value ~ Mixture2, data =lm3)
summary(res.aov)
#Df Sum Sq Mean Sq F value Pr(>F)
#Mixture 4 0.027539 0.006885 31.05 1.29e-05 ***
# Summary of the analysis
TukeyHSD(res.aov)
#$Mixture
#diff lwr upr p adj
#siCREB5-siCont 0.018232217 -0.02178207 0.058246501 0.5848284
#siHIVEP2-siCont 0.022666666 -0.01734762 0.062680950 0.3926589
#siNFAT5-siCont -0.095008203 -0.13502249 -0.054993919 0.0001096
#siPPARA-siCont -0.020503073 -0.06051736 0.019511211 0.4822554
#siHIVEP2-siCREB5 0.004434449 -0.03557983 0.044448733 0.9955843
#siNFAT5-siCREB5 -0.113240420 -0.15325470 -0.073226136 0.0000234
#siPPARA-siCREB5 -0.038735290 -0.07874957 0.001278994 0.0588599
#siNFAT5-siHIVEP2 -0.117674869 -0.15768915 -0.077660585 0.0000165
#siPPARA-siHIVEP2 -0.043169740 -0.08318402 -0.003155456 0.0334159
#siPPARA-siNFAT5 0.074505130 0.03449085 0.114519414 0.0008193