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
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
|
getDEGs <- function(exprset,genelist,degtype = "deseq2"){
library(tidyverse)
library(tinyarray)
if(!(degtype %in% c("deseq2","limma","dger"))){
stop("请输入正确的差异分析方法")
}
exprset[exprset<0] <- NA
exprset <- na.omit(exprset)
group_list <- tinyarray::make_tcga_group(exprset)
pca.plot <- draw_pca(exprset,group_list)
pca.plot <- pca.plot + theme_classic()
Cairo::CairoPDF(file = "pcaplot.pdf",width = 10,height = 10)
print(pca.plot)
dev.off()
if(max(exprset)<50){exprset <- ceiling(2^exprset - 1)}
valid_gene <- apply(exprset,1, function(x){sum(x==0) <= 0.25*ncol(exprset)})
exprset <- exprset[valid_gene,]
colData <- data.frame(row.names = colnames(exprset),condition = group_list)
valid_genes <- dplyr::intersect(genelist,row.names(exprset))
if(!length(valid_genes)){stop("你的基因在表达矩阵里全都不存在")}
if(degtype=="deseq2"){
library(DESeq2)
dds <- DESeqDataSetFromMatrix(
countData = exprset,
colData = colData,
design = ~condition)
dds <- DESeq(dds)
res <- results(dds, contrast = c("condition",rev(levels(group_list))))
resOrdered <- res[order(res$pvalue),]
DEG <- as.data.frame(resOrdered)
DEG <- na.omit(DEG)
detach(package:DESeq2)
logFC_cutoff <- 2
DEG$change = as.factor(
ifelse(DEG$pvalue < 0.05 & abs(DEG$log2FoldChange) > logFC_cutoff,
ifelse(DEG$log2FoldChange > logFC_cutoff ,'UP','DOWN'),'NOT'))
}
genelist_up <- DEG %>% filter(.,change == "UP") %>% filter(.,rownames(.) %in% all_of(genelist))
genelist_down <- DEG %>% filter(.,change == "DOWN") %>% filter(.,rownames(.) %in% all_of(genelist))
print(paste(("总上调表达基因有"),nrow(DEG[DEG$change=="UP",]),sep = ""))
print(paste(("总下调表达基因有"),nrow(DEG[DEG$change =="DOWN",]),sep = ""))
print(paste(("上调表达的重要基因有"),nrow(genelist_up),sep = ""))
print(paste(("下调表达的重要基因有"),nrow(genelist_down),sep = ""))
write.csv(DEG,"DESeq2_DEG.csv")
write.csv(group_list,"group_list.csv")
write.csv(colData,"colData.csv")
write.csv(genelist_up,"genelist_up.csv")
write.csv(genelist_down,"genelist_down.csv")
library(EnhancedVolcano)
E1 <- EnhancedVolcano(DEG,
lab = row.names(DEG),
selectLab = genelist,
x = 'log2FoldChange',
xlab = bquote(~Log[2]~ 'Fold Change'),
y = 'padj',
title = 'Tumor versus Normal',
pCutoff = 10e-2,
FCcutoff = 2,
pointSize = 2.0,
labSize = 2.0,
legendPosition = 'right',
legendLabSize = 12,
drawConnectors = TRUE,
widthConnectors = 0.5,
col = c('grey', 'green', 'darkgreen', 'red'),
colAlpha = 0.5,
colConnectors = 'black')
Cairo::CairoPDF(file = "DEGVolcano.pdf",width = 10,height = 10)
print(E1)
dev.off()
detach(package:EnhancedVolcano)
colData <- colData %>%
mutate(colData,id = row.names(colData)) %>%
arrange(colData,desc(condition))
dat <- exprset %>% select(.,all_of(colData$id),everything())
dat_group_list <- colData$condition
cg1 <- rownames(DEG)[DEG$change !="NOT"]
cg1 <- dplyr::intersect(cg1,genelist)
if(length(cg1)==0){
print("你的基因列表里没有差异表达基因，但还是给你画个热图。")
h1 = draw_heatmap(dat[genelist,],
dat_group_list,
cluster_cols = FALSE,
legend = TRUE,
show_rownames = TRUE,
split_column = FALSE,
annotation_legend = TRUE,
n_cutoff = 2)
Cairo::CairoPDF(file = "DEGHeatmap.pdf",width = 10,height = 10)
print(h1)
dev.off()
}
if(length(cg1)>10){
cg1 <- cg1[1:10]
h1 = draw_heatmap(dat[cg1,],
dat_group_list,
cluster_cols = FALSE,
legend = TRUE,
show_rownames = TRUE,
split_column = FALSE,
annotation_legend = TRUE,
n_cutoff = 2)
Cairo::CairoPDF(file = "DEGHeatmap.pdf",width = 10,height = 10)
print(h1)
dev.off()
}else{
h1 = draw_heatmap(dat[cg1,],
dat_group_list,
cluster_cols = FALSE,
legend = TRUE,
show_rownames = TRUE,
split_column = FALSE,
annotation_legend = TRUE,
n_cutoff = 2)
Cairo::CairoPDF(file = "DEGHeatmap.pdf",width = 10,height = 10)
print(h1)
dev.off()
}
TNcolData <- colData
return(TNcolData)
print("File Saved. Unnecessary Packages detached. Proceed To Next Step.")
}
|