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

BiocManager::install("GEOquery") # 下载GEOquery包
library(GEOquery)
library(tidyverse)
setwd('C:/Users/30521/Desktop/R自学/GEO-analysis') # 设置工作路径；注意只能是“/”
gset<-getGEO("GSE299988",destdir=".",AnnotGPL=F,getGPL=F)
# getGEO(GEO,GSEMatrix = TRUE,AnnotGPL = FALSE,destdir = ".",getGPL = TRUE)，这是getGEO()函数的使用格式。
# GSEMatrix指是否下载整理好的 Series Matrix 表达矩阵
# AnnotGPL指是否同时下载平台注释表（GPL）
# destdir指数据保存到本地的路径（默认为"."）
# getGPL指是否自动下载对应的 GPL 平台文件。

class(gset) # 读取的为list类型
gset[[1]] # 查看信息，GSE只有一个平台就是1，可能有列表中其他平台，要查看就是更改成gset[[2]]等等
exp<-exprs(gset[[1]])# 提取基因表达矩阵数据，由于没有载入注释文件，矩阵的行名还是探针名

# 1、标准化
boxplot(exp,outline=F,las=2) # 绘制表达量分布箱线图，快速检查数据的整体质量和标准化效果。在同一条直线上则是已经标准化了
# boxplot(x, outline, las) x为表达矩阵，outline为是否显示异常值点，las为控制x轴标签方向：1 = 水平，2 = 垂直
library(limma)
exp<-normalizeBetweenArrays(exp) # 手动标准化，再去查看箱图

# 2、归一化
range(exp) # 查看range范围，数值过大需要归一化，log2处理
exp<-log2(exp+1)

# 3、基因注释，GPL导入
gpl<-read.table("GPL21185-21174.txt",header=T,fill=T,sep="\t",comment.char="#",stringsAsFactors=F,quote='')
# header=T列名保留
# fill=T文件中某些行的列数不一致时，自动用空值填充缺失的列
# sep="\t"指定列之间的分隔符是制表符（Tab）
# comment.char="#"忽略文件中以 # 开头的注释行
# stringsAsFactors=F把字符型数据（如探针 ID、基因名）保留为字符串，而不是自动转换成因子
# quote=''告诉 R 文件中没有引号包裹的字符串，避免把基因名或 ID 中的特殊字符误判为引号
colnames(gpl) # 查看列名
ids<-gpl[,c("ID","GENE_SYMBOL")] # 选取查看的两列
colnames(ids)<-c("probe_id","symbol")
library(stringr)
ids$symbol<-trimws(str_split(ids$symbol,"//",simplify=T)[,3])# 如果是第3个是基因，就选择第3列
# stringr包是tidyverse生态里专门处理字符串的工具包
# str_split()是按指定分隔符拆分字符串
# simplify=T默认返回列表，设为TRUE后会直接返回矩阵
# trimws()去除提取出的第 3 个基因名两端的空白字符，确保基因名格式干净
ids<-ids[ids$symbol!='',] # 删除了空白探针
ids<-ids[ids$symbol!='---',]
ids<-ids[ids$symbol!='--- ',] # 具体情况具体格式删除
ids<-ids[ids$probe_id %in% rownames(exp),] # 保留ids中与exp匹配的行
exp2<-exp[ids$probe_id,] # 选取exp中与ids相互匹配的行
table(rownames(exp2)==ids$probe_id)
exp3<-cbind(ids,exp2) # 合并两者为数据框

# 4、基因去重（一个基因会有多个探针）
exp3<-aggregate( . ~symbol,exp3,max) # 合并同一基因的多探针表达量
# . ~ symbol去重的核心逻辑，“.”代表除了symbol列之外的所有列，~symbol代表按照symbol列进行分组
# max指取每个基因对应探针的最大表达量作为代表值，也可以选择mean平均
rownames(exp3)<-exp3$symbol
exp3<-exp3[,-c(1,2)] # 再删去合并的前两列

# 5、获取临床数据
library(Biobase) # pData是Biobase里面的函数
rt1<-pData(gset[[1]]) # 先前exprs函数是获取的gset中基因表达矩阵
# pData函数用来从 ExpressionSet 对象中提取 phenotype data（表型数据），也就是每个样本的临床信息，比如分组、性别、年龄、疾病分期等。
write.csv(rt1,file='GSE299988_cli.csv',row.names=T)
#（1）疾病分组
colnames(rt1)
group_list<-ifelse(rt1$`tissue types:ch1`=='Normal','Normal','Tumor')
# 将Normal与Tumor细胞分开
# 更多组别的分组：group_list<-ifelse(str_detect(rt1$`title`,'normal_adjace'),'Normal',ifelse(str_detect(rt1$`title`,'positiveLN_RAIrefractive'),'P_tumor','N_tumor'))
# str_detect可以模糊字符选择
rt1$group<-group_list
rt2<-rt1 %>% dplyr::select("group")
# %>%管道操作符，把左边的rt1作为右边函数的输入
# dplyr::select()用于从数据框中选择特定的列
# （2）分组排列
rt3<-rt2 %>% arrange(rt2) # 升序排列，Normal在前，Tumor在后
# rt3<-rt2 %>% arrange(desc(group)) Tumor在前，Normal在后

# 6、数据保存
# exp3与rt3的列行一一对应
exp3<-exp3[,rownames(rt3)]
identical(rownames(rt3),colnames(exp3)) # 返回为T确认一一对应
save(rt3,exp3,file="GSE299988.rda") # 保存数据
# 之后可以直接load("GSE299988.rda")进行使用

 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

library(GEOquery)
library(tidyverse)
setwd('C:/Users/30521/Desktop/R自学/GEO-analysis/GSE200539')
# 中间部分与上述GEO数据处理部分一致
# 基因注释，soft导入（没有Download full table…的情况下）
GPL<-getGEO(filename="GPL23159_family.soft.gz",destdir=".")
gpl<-GPL@dataTable@table

# 由于这个注释文件没有Gene symbol，有NM_序列，因此提取这一列
ids<-gpl[,c("ID","SPOT_ID.1")] # 选取查看的两列
colnames(ids)<-c("probe_id","symbol") # 重命名
library(stringr)
ids$symbol<-trimws(str_split(ids$symbol,"//",simplify=T)[,1]) # 选取//分割的第一列NM_序列
# 下载两个包：BiocManager::install(c("org.Hs.eg.db", "AnnotationDbi"))
# org.Hs.eg.db包是人类基因注释数据库包，Hs 代表 Homo sapiens（人），eg 代表 Entrez Gene。包含了人类基因的各种标识符（如 Entrez Gene ID、Gene Symbol、Ensembl ID 等）之间的映射关系
# AnnotationDbi包是所有注释数据库包的核心接口包
library(org.Hs.eg.db)
library(AnnotationDbi)
symbol_ids<-mapIds(x=org.Hs.eg.db,keys=ids$symbol,keytype="REFSEQ",column="SYMBOL",multivals="first")
# org.Hs.eg.db指定要使用的注释数据库，这里是人类基因注释库
# keytype告诉输入的类型
# column指定想要转换后的目标 ID 类型
# multivals="first" 处理 “一个输入 ID 对应多个输出 ID” 的情况。当一个 REFSEQ ID 对应多个基因符号时，只取第一个结果，避免返回列表，保证输出是一个向量
symbol_ids
ids$symbol<-symbol_ids # NM_序列转换为基因名
ids<-na.omit(ids) # 删除NA

# 注意：exp行名默认是字符，如果ids$symbol是纯数字，应当转变为字符格式ids$probe_id<-as.character(ids$probe_id)

exp2<-exp[ids$probe_id,] # 选取exp中与ids相互匹配的行
table(rownames(exp2)==ids$probe_id)
exp3<-cbind(ids,exp2) # 合并两者为数据框

# 之后与之前GEO数据处理部分一致

 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

setwd('C:/Users/30521/Desktop/R自学/GEO-analysis/GSE200539/Affymetrix原始数据处理')
library(tidyverse)
# 下载：BiocManager::install("oligo")
# Bioconductor 生态里专门用来处理 Affymetrix 芯片原始数据的 R 包
library(oligo)
# 1、数据导入
FileName<-list.files(path='./GSE200539_RAW') # 注意提前解压
FileName
sample_id<-str_sub(FileName,1,10) # 重命名，此处选取前10个字符，即命名为如GSM6037488的样式
rt<-data.frame(FileName=FileName,row.names=sample_id,group=rep(c('CTL','IDMF','DRF','MMF'),each=3))
# 分组数据参考网站上的Sample的分组，共4组，每组3个样本
cel_files<-list.celfiles('./GSE200539_RAW',full.names=T,listGzipped=T)
# list.celfiles()是oligo包的核心函数，用于搜索指定目录下的所有 .CEL文件
# './GSE200539_RAW'当前目录下的GSE200539_RAW文件夹
# full.names=T返回完整的文件路径（包含文件夹名）
# listGzipped=T支持识别并列出被gzip压缩的 .CEL.gz文件
cel_files
cel_data<-read.celfiles(cel_files) # 通过cel_files获取的所有文件路径导入并读取

# 2、质控处理
# 下载：BiocManager::install('arrayQualityMetrics')
# arrayQualityMetrics包是Bioconductor中专门用于芯片数据质量控制的工具包，能自动生成详细的质控报告，包含多种图表来评估数据质量。
library(arrayQualityMetrics)
# 版本问题手动下载吧……然后根据指示下一堆依赖包
library(GEOquery)
expr<-exprs(cel_data) #读取数据为原始表达量的矩阵
colnames(expr)<-sample_id # 列名重命名
eset<-ExpressionSet(assayData=expr,phenoData=AnnotatedDataFrame(rt))
# ExpressionSet()是Biobase包的核心函数，用来整合表达数据、样本信息和注释信息
# AnnotatedDataFrame()把rt数据框转换成 ExpressionSet 要求的注释数据框格式
# 此处是转换为ExpressionSet用于质控分析
arrayQualityMetrics(eset,outdir='GSE200539_QulityMetrics',force=T,do.logtransform=T,intgroup='group')
# force=T如果输出文件夹已经存在，强制覆盖它，避免报错
# do.logtransform=T对表达数据进行对数转换，这是芯片数据分析的常规预处理步骤
# intgroup='group'指定用于分组的表型信息列（比如 rt 数据框里的 group 列）

#rma标准化,对原始 CEL 文件数据进行背景校正、归一化和汇总
eset_rma<-oligo::rma(cel_data)
class(eset_rma) # "ExpressionSet"
exp<-exprs(eset_rma)
colnames(exp)<-str_sub(colnames(exp),1,10) # 保留列名前面10个字符
boxplot(exp,outline=F,notch=T,las=2) # 箱图看结果

# 这里的exp与之前章节得到的直接读取的exp相同
# 之后部分是基因注释，与先前相同

 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

#清空环境
rm(list=ls())
# BiocManager::install("GEOquery")
library(GEOquery)
library(tidyverse)
library(limma)
setwd('E:/R_do/阿伦生信教学代码/GEO/ilumina原始数据处理/') #设置工作路径
gset <- getGEO("GSE68004", destdir = ".", AnnotGPL = F, getGPL = F)#手动下载，自动读取
class(gset) 
gset[[1]] #查看信息
#提取基因表达数据  
exp <- exprs(gset[[1]])
#1.标准化####
boxplot(exp, outline=F,las=2)
library(limma)
exp <- normalizeBetweenArrays(exp)
boxplot(exp, outline=F,las=2)
#2.归一化####
range(exp)#查看range范围，如果数值过大，需要归一化，log2处理
range(exp,na.rm = T)
# [1]   -22.82553 17174.43372


#1.观察数据结构####
dat1 <- data.table::fread('GSE68004_non-normalized_ProbeRowSignalDetectionDataset_5383_143523_1.txt')
dat2 <- data.table::fread('GSE68004_non-normalized_ProbeRowSignalDetectionDataset_5383_143529_1.txt')
#2.读取数据####
x1 <- read.ilmn(files = 'GSE68004_non-normalized_ProbeRowSignalDetectionDataset_5383_143523_1.txt',
                expr = '750',#数据列名 相同的内容
                other.columns='Detection', #p值
                probeid='ID_REF') #芯片列名


x2 <- read.ilmn(files = 'GSE68004_non-normalized_ProbeRowSignalDetectionDataset_5383_143529_1.txt',
                expr = '586',#数据列名 相同的开头
                other.columns='Detection', #p值
                probeid='ID_REF') #芯片列名
#合并
pro_id <- intersect(rownames(x1),rownames(x2))
x1 <- x1[pro_id,]
x2 <- x2[pro_id,]
x <- cbind(x1, x2)

#3.校准化-log处理####
y <- neqc(x) 
#4.提取表达数据####
exp <- y$E
boxplot(exp,las=3,outline=F)
range(exp)
#[1]  4.161602 14.037369
#去批次
library(sva)
batch <- c(rep(1, ncol(x1)), rep(2, ncol(x2))) #输入批次信息
exp_combat <- ComBat(dat = exp, batch = batch, mod = NULL, par.prior = TRUE)
boxplot(exp_combat,las=3,outline=F)
range(exp_combat)
#[1]  3.807521 14.370650
#5.整理临床数据
rt <- pData(gset[[1]])
colnames(rt)
rt <- dplyr::select(rt,'title','source_name_ch1')
sample_id <- substring(rt$title,nchar(rt$title)-8)#提取字符-8就是最后8个
sample_id
rownames(rt) <- sample_id

rt$group <- trimws(str_split(rt$source_name_ch1,',',simplify = T)[,2])
rt <- dplyr::select(rt,group)

#保持行名列名一致
exp_combat <- exp_combat[,rownames(rt)]
identical(rownames(rt),colnames(exp_combat))
#[1] TRUE
save(exp_combat,rt,file = 'exp_combat.rad')

#质控检查
rt1 <- rt[colnames(x$E), , drop = FALSE]
library(arrayQualityMetrics)
colnames(x$E) <- rownames(rt1)
eset <- ExpressionSet(assayData = x$E,phenoData = AnnotatedDataFrame(data = rt1))
arrayQualityMetrics(eset,
                    outdir='GSE68004_quality',
                    force=T,
                    intgroup='group',
                    do.logtransform=T)


rm(list=ls())
load('exp_combat.rad')

 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

# control type: 0是对照，1是基因
# glsaboveBG:信号强度是否高于背景，结果是1则说明值是可靠的

# 1、原始数据读取
setwd("C:/Users/30521/Desktop/R自学/GEO-analysis/GSE211806")
library(limma)
library(tidyverse)
FileName<-list.files(path="./GSE211806_RAW")
FileName
sample_id<-str_sub(FileName,1,10)
tissue_type<-rep(c("PLC/PRF/5","SNU-449"),each=12)
group<-rep(c("control","TGFb","Galunisertib","TGFb + Galunisertib"),each=3)

Ids<-data.frame(FileName=FileName,
                sample_id=sample_id,
                tissue_type=tissue_type,
                group=group,
                row.names=sample_id,
                stringsAsFactors=F)
table(Ids$group)
write.table(Ids,file="GSE211806_Ids.txt",quote=F,sep="\t",row.names=F)
save(Ids,file="GSE211806_Ids.rda")

GSE211806_raw<-read.maimages(files=Ids$FileName,
                             source="agilent", # 代表经过哪种程序得到的，有些Agilent芯片是通过genepix处理的
                             path="GSE211806_RAW",
                             names=Ids$sample_id,
                             other.columns="gIsWellAboveBG", # 读取进去是为了判断是否高于背景值
                             green.only=T) # 代表是个单色芯片(Cy3)，双色芯片F(Cy3,Cy5)，见官方样本信息页的Label
# GSE211806_raw的E代表expression，Eb代表background
# 添加标签
GSE211806_raw$Ids<-Ids

# probe type summary
table(GSE211806_raw$genes$ControlType)
# 0说明不是control，需要检测这60901个基因，1是阳性对照，-1是阴性对照
head(GSE211806_raw$E)[,1:5]
dim(GSE211806_raw)

# 保存数据
save(GSE211806_raw,file="GSE211806_raw_Agilent.rda")
boxplot(log2(GSE211806_raw$E),ylab=expression(log[2](intensity)),las=2,
        col=factor(Ids$group),outline=F)
source('Function/PCA.R')
PCA(log2(GSE211806_raw$E),
    ntop=nrow(GSE211806_raw$E),
    group=Ids$group,
    show_name=F)# 或者T

library(arrayQualityMetrics)

# 函数不接受ElistRaw对象，只接受ExpressionSet
library(GEOquery)
rownames(Ids)<-Ids$sample_id
eset<-ExpressionSet(assayData=GSE211806_raw$E,
                    phenoData=AnnotatedDataFrame(data=Ids))

arrayQualityMetrics(eset,
                    outdir="GSE211806_QulitMetrics",
                    do.logtransform=T,
                    force=T,
                    intgroup="group")

# 个体差异较大
# 背景校正及标准化
bgc<-backgroundCorrect(RG=GSE211806_raw,
                       method="normexp",
                       offset=50,
                       normexp.method="mle")

normdata<-normalizeBetweenArrays(bgc,method="quantile")
# normalizeBetweenArrays适用于单色芯片，normalizeinArrays适用于双色芯片
sum(is.na(normdata$E))
boxplot(normdata$E,ylab=expression(log[2](intensity)),las=2,
        col=factor(Ids$group),outline=F)
PCA(normdata$E,
    ntop=nrow(normdata$E),
    group=Ids$group,
    show_name=F)# 或者T

# 可再做一步质控（选）
eset2<-ExpressionSet(assayData=normdata$E,
                    phenoData=AnnotatedDataFrame(data=Ids))

arrayQualityMetrics(eset2,
                    outdir="GSE211806_QulitMetrics_norm",
                    do.logtransform=T,
                    force=T,
                    intgroup="group")
exp<-exprs(eset2)
# 行名问题
rownames(normdata$E)<-normdata$genes$ProbeName
# 选取探针id的行，舍去前三行
exp<-normdata$E[4:nrow(normdata$E),]
rt<-dplyr::select(Ids,"tissue_type","group")
save(exp,normdata,Ids,rt,file='GSE211806_Agilent_processed.rda')

 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

setwd('C:/Users/30521/Desktop/R自学/GEO-analysis/GSE309139')
library(tidyverse)
exp<-data.table::fread("GSE309139_raw_counts.txt",header=T)
# 用data.table::fread()高效读取RNA-seq原始计数文件
# header=T 确保文件第一行（样本ID）作为列名
table(is.na(exp)) # 查看是否有缺失值
exp[is.na(exp)]<-0
table(is.na(exp)) # 查看是否有缺失值
exp<-column_to_rownames(exp,var='Gene')
# 把基因名那一列挪到最前面成为行名

# 1、设置分组信息
sample_id<-colnames(exp)
sample_id
group<-rep(c('L','H','LL','HL'),each=6) # 具体分组看具体数据
group
rt<-data.frame(row.names=sample_id,group=group)
identical(rownames(rt),colnames(exp))

# 2、DESeq2分析
# BiocManager::install("DESeq2")
library(DESeq2)
rt$group<-factor(rt$group) # 变成因子的格式
rt$group
# 构建DESeq2对象
dds<-DESeqDataSetFromMatrix(countData=exp,colData=rt,design=~group)
# design=~group根据group分组来比较差异
# dds把原始计数数据、样本分组信息、实验设计逻辑打包在一起，后续所有分析（标准化、差异检验）都基于这个对象

# 低表达基因过滤（自定义，可选）
filter<-rowSums(counts(dds))>1 # 筛选每行加起来大于1的基因
dds<-dds[filter,]
# 差异表达分析
dds<-DESeq(dds)
# 格式为c("因子名"，“实验组”，“对照组”)
res<-results(dds,contrast=c("group","H","L"))
# 提取差异表达结果
res1<-as.data.frame(res)
# log2FoldChange差异倍数的对数值：H 组表达量 / L 组表达量 的 log₂ 结果
# lfcSE是log2FC 的标准误，代表差异倍数的统计波动程度（值越小，log2FC 越可靠）
# pvalue越小越拒绝原假设H0，越显著性差异
# padj是校准后的P值，控制多检验的假阳性
# 结果筛选padj<0.05 & abs(log2FoldChange)>=1，选出差异表达
degs<-filter(res1,padj<0.05 & abs(log2FoldChange)>=1)
view(degs)

# 绘制离散图
plotDispEsts(dds)
# 离散图横轴是基因经过标准化后的平均表达量，纵轴离散度，衡量基因表达量的变异程度
# P值直方图，看P值显著基因多少
hist(res$padj)

# 3、count数据vst转化，这里得到exp1的数据就不是整数了，可以理解为之前章节的exp，之后就是一样的标准化操作等等
vst<-vst(dds,blind=F)
# vst()将RNA-seq的原始计数数据转换为适合可视化和下游分析的数值
# blind=F变换时考虑实验设计信息，不会忽略分组差异，这样能更真实地保留生物学差异
exp1<-assay(vst)
range(exp1)
exp1<-as.data.frame(exp1)
save(exp,exp1,rt,file='GSE309139.rda')

# 注意，这里得到exp1之后记得看看需要标准化，归一化那一系列操作，之后再进行热图绘制。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16

setwd('C:/Users/30521/Desktop/R自学/GEO-analysis/limma差异分析')
library(limma)
library(pheatmap)
library(tidyverse)
library(ggpubr)

# 数据读取，直接使用之前第一部分处理好的数据
load('C:/Users/30521/Desktop/R自学/GEO-analysis/GSE299988/GSE299988.rda')
rt<-rt3
exp<-exp3
str(exp) # 数据全是字符型：$ GSM9051410: chr  "8.561352698" "5.812654031" "2.853925376" "13.89898258" ...
exp<-exp %>% mutate(across(everything(),as.numeric))
# 把数据框 exp 里所有列的数据类型都转换成数值型
str(exp) # 把数据转换成了数字型：$ GSM9051410: num  8.56 5.81 2.85 13.9 8.65 ...
identical(rownames(rt),colnames(exp))
# 确认行名列名一致

 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

group_list<-factor(rt$group,levels=c("Normal","Tumor"))
# 把临床数据的分组弄成因子形式
design<-model.matrix(~group_list)
# model.matrix()：把因子型的分组信息转化为统计模型可识别的数值矩阵
# 基于上面的 group_list，生成的设计矩阵会以 Normal 作为截距（基线），用一列数值来表示 Tumor 组与基线的差异
design
# 线性模型拟合
fit<-lmFit(exp,design)
# 贝叶斯检验
fit2<-eBayes(fit)
deg<-topTable(fit2,coef=2,number=Inf)
# 是差异表达基因结果数据框，包含了所有基因的差异分析结果
DEG=na.omit(deg) # 去除可能的NA值
write.table(DEG,file='DEG1.txt',sep="\t",row.names=T,col.names=NA,quote=F)
# quote=F：不在字符串两端添加引号，避免后续读取时出现多余的引号干扰

# 通常logFC设置为1，adj.P.Val<0.05，下述代码进行显著水平与正负相关性的区分
logFC_cutoff<-1
type1=(DEG$adj.P.Val<0.05)&(DEG$logFC < -logFC_cutoff)
type2=(DEG$adj.P.Val<0.05)&(DEG$logFC > logFC_cutoff)
DEG$change=ifelse(type1,"DOWN",ifelse(type2,"UP","NOT"))
write.table(DEG,file='DEG2.txt',sep="\t",row.names=T,col.names=NA,quote=F)
table(DEG$change) # 统计上述3中情况："DOWN","UP","NOT"的个数

library(ggplot2)
library(cowplot)
library(ggrepel)
DEG<-mutate(DEG,Gene_symbol=row.names(DEG)) # 把Gene_symbol弄到最后一列
UPDEG<-subset(DEG,change=="UP")
UPDEG_5<-top_n(x=UPDEG,n=-5,wt=P.Value)
# n = -5：- 表示按升序取最小的 5 个，+ 表示按降序取最大的 5 个
# wt = P.Value：按 P.Value 列来排序，P 值越小代表差异越显著
# 由此选取最小的5组就是最小的5个P，最显著的5个
DOWNDEG<-subset(DEG,change=="DOWN")
DOWNDEG_5<-top_n(x=DOWNDEG,n=-5,wt=P.Value)
p <- ggplot(data = DEG,
            aes(x = logFC,
                y = -log10(adj.P.Val))) +
  geom_point(alpha = 0.5, size = 4.5,
             aes(color = change)) +
  ylab("-log10(adj.P.Val)") +
  xlab("log2(Fold Change)") +
  scale_color_manual(values = c("#1F77B4", "grey", "#FF7F0E")) +
  geom_vline(xintercept = c(-1, 1), lty = 5, col = "black", linewidth = 0.8) +
  geom_hline(yintercept = -log10(0.05), lty = 4, col = "black", linewidth = 0.8) +
  theme_half_open() +
  # 标记最显著的5个上调基因
  geom_text_repel(data = UPDEG_5, 
                  aes(label = Gene_symbol), 
                  vjust = 1.5,
                  box.padding = unit(0.35, "lines"),
                  point.padding = unit(0.35, "lines"),
                  segment.color = 'grey50',
                  nudge_y = 0.1,
                  direction = "y") +
  # 标记最显著的5个下调基因
  geom_text_repel(data = DOWNDEG_5, 
                  aes(label = Gene_symbol), 
                  vjust = 1.5,
                  box.padding = unit(0.35, "lines"),
                  point.padding = unit(0.35, "lines"),
                  segment.color = 'grey50',
                  nudge_y = 0.1,
                  direction = "y")
p # 查看火山图
# 保存为pdf
pdf(file="volcano.pdf",width=10,height=6)
print(p)
dev.off()

 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

library(pheatmap)
identical(rownames(rt),colnames(exp))
DEG<-read.table("DEG2.txt",sep="\t",check.names=F,stringsAsFactor=F,header=T,row.names=1)
table(DEG$change) 
diff<-rownames(DEG)[DEG$change != "NOT"] # 选出仅有显著上调或下调的基因
exp_diff<-exp[diff,]
range(exp_diff)

# 补：从DEG中选出目标基因来做热图
DEG_target <- DEG[grepl(pattern = "CASP1|ZDHHC|NLRP3|IL1B|IL18|GSD", x = DEG$name, ignore.case = TRUE), ]
target<-exp[rownames(DEG_target),]
# 作图
# 绘制差异基因表达热图
map<-pheatmap(
  exp_diff,
  annotation_col = rt, # 列注释（样本分组信息）
  color = colorRampPalette(c("blue", "white", "red"))(100), # 颜色渐变（100个梯度）
  cluster_cols = F, # 不对样本（列）进行聚类
  cluster_rows = T,  # 对基因（行）进行聚类
  show_colnames = F, # 不显示样本名（避免拥挤）
  show_rownames = T, # 显示基因名
  border_color = NA, # 去掉单元格边框
  scale = 'row', # 按行标准化（让每个基因的表达在不同样本间可比较）
  treeheight_row = 20, # 行聚类树的高度（可选）
  fontsize_row = 8, # 行名字体大小（可选）
  fontsize_col = 8 # 列名字体大小（可选）
)
pdf(file='heatmap.pdf',width=6,height=20)
print(map)
dev.off()

# 特定展示（可选）
library(tidyverse)
diff=DEG[DEG$change != "NOT",]
up<-diff %>% top_n(25,logFC)
dw<-diff %>% top_n(-25,logFC)
all<-c(rownames(up),rownames(dw))
exp_diff2<-exp[all,]
map2<-pheatmap(
  exp_diff2,
  annotation_col = rt, # 列注释（样本分组信息）
  color = colorRampPalette(c("blue", "white", "red"))(100), # 颜色渐变（100个梯度）
  cluster_cols = F, # 不对样本（列）进行聚类
  cluster_rows = T,  # 对基因（行）进行聚类
  show_colnames = F, # 不显示样本名（避免拥挤）
  show_rownames = T, # 显示基因名
  border_color = NA, # 去掉单元格边框
  scale = 'row', # 按行标准化（让每个基因的表达在不同样本间可比较）
  treeheight_row = 20, # 行聚类树的高度（可选）
  fontsize_row = 8, # 行名字体大小（可选）
  fontsize_col = 8 # 列名字体大小（可选）
)
pdf(file='heatmap2.pdf',width=6,height=8)
print(map2)
dev.off()

 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

setwd('C:/Users/30521/Desktop/R自学/GEO-analysis/GSE309139')
library(tidyverse)
exp<-data.table::fread("GSE309139_raw_counts.txt",header=T)
# 用data.table::fread()高效读取RNA-seq原始计数文件
# header=T 确保文件第一行（样本ID）作为列名
table(is.na(exp)) # 查看是否有缺失值
exp[is.na(exp)]<-0
table(is.na(exp)) # 查看是否有缺失值
exp<-column_to_rownames(exp,var='Gene')
# 把基因名那一列挪到最前面成为行名

# 1、设置分组信息
sample_id<-colnames(exp)
sample_id
group<-rep(c('L','H','LL','HL'),each=6) # 具体分组看具体数据
group
rt<-data.frame(row.names=sample_id,group=group)
identical(rownames(rt),colnames(exp))

# 2、DESeq2分析
# BiocManager::install("DESeq2")
library(DESeq2)
rt$group<-factor(rt$group) # 变成因子的格式
rt$group
# 构建DESeq2对象
dds<-DESeqDataSetFromMatrix(countData=exp,colData=rt,design=~group)
# design=~group根据group分组来比较差异
# dds把原始计数数据、样本分组信息、实验设计逻辑打包在一起，后续所有分析（标准化、差异检验）都基于这个对象

# 低表达基因过滤（自定义，可选）
filter<-rowSums(counts(dds))>1 # 筛选每行加起来大于1的基因
dds<-dds[filter,]
# 差异表达分析
dds<-DESeq(dds)
# 格式为c("因子名"，“实验组”，“对照组”)
res<-results(dds,contrast=c("group","H","L"))
# 提取差异表达结果
res1<-as.data.frame(res)
# log2FoldChange差异倍数的对数值：H 组表达量 / L 组表达量 的 log₂ 结果
# lfcSE是log2FC 的标准误，代表差异倍数的统计波动程度（值越小，log2FC 越可靠）
# pvalue越小越拒绝原假设H0，越显著性差异
# padj是校准后的P值，控制多检验的假阳性
# 结果筛选padj<0.05 & abs(log2FoldChange)>=1，选出差异表达
write.table(res1,file='DEG1.txt',sep='\t',row.names=T,col.names=NA,quote=F)
degs<-filter(res1,padj<0.05 & abs(log2FoldChange)>=1)
view(degs)

 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

library(tidyverse)
library(clusterProfiler)
library(org.Hs.eg.db)
library(enrichplot)
library(ggplot2)
setwd('C:/Users/30521/Desktop/R自学/GEO-analysis/富集分析')
DEG<-read.table("DEG2.txt",sep="\t",check.names=F,stringsAsFactors=F,header=T,row.names=1)
table(DEG$change)
diff<-rownames(DEG)[DEG$change!='NOT']

# 将基因ID从Symbol转换为Entrez ID（富集分析主要用Entrez ID 或 Ensembl ID来索引基因，因此要进行转换）
gene_entrez <- bitr(diff, fromType ="SYMBOL", toType ="ENTREZID", OrgDb = org.Hs.eg.db)
# OrgDb = org.Hs.eg.db：使用人类（Homo sapiens）的注释数据库
cat("所有基因ID转换", nrow(gene_entrez),"/", length(diff)) # 看转换成功多少个
untrans <- setdiff(diff, gene_entrez$SYMBOL) # 两部分的差，得出哪些转换成功
untrans 
isValid <- "C9orf92" %in% keys(org.Hs.eg.db, keytype = "SYMBOL")
isValid # 查看T还是F，判断没有被转换的基因是否被收录在org.Hs.eg.db里面

# 使用mapIds函数转换
library(AnnotationDbi)
entrez_ids <- mapIds(
  x = org.Hs.eg.db,
  keys = untrans,       # 输入的基因别名列表
  keytype = "ALIAS",    # 明确指定输入类型为别名（ALIAS）
  column = "ENTREZID",  # 转换为ENTREZID作为中间步骤
  multiVals = "first"   # 多个匹配时取第一个
)
entrez_df <- data.frame(ENTREZID = entrez_ids, check.names = FALSE)
entrez_df <- na.omit(entrez_df)
entrez_df <- entrez_df %>% rownames_to_column(var = 'SYMBOL')
#合并
gene_entrez2 <- rbind(gene_entrez,entrez_df)

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21

GO <- enrichGO(gene = gene_entrez2$ENTREZID,
               OrgDb = org.Hs.eg.db,
               keyType ="ENTREZID",
               ont ="ALL", # "BP", "MF", "CC"或"ALL"
               pAdjustMethod ="BH",
               pvalueCutoff =0.05,
               qvalueCutoff =0.05,readable = T)
# GO 分为三个分支：BP（Biological Process，生物学过程）；MF（Molecular Function，分子功能）；CC（Cellular Component，细胞组分）；ALL 表示同时分析三个分支
# pAdjustMethod = "BH"：多重检验校正方法，BH（Benjamini-Hochberg）是最常用的方法，用于控制假阳性率
# readable = T：将结果中的Entrez ID自动转换为更易读的基因Symbol，方便解读
GO_res <- GO@result
write.table(GO_res,file="GO_res.txt",sep="\t",quote=F,row.names = F)
#简单作图
barplot(GO, showCategory = 20)
dotplot(GO, showCategory = 20)
#分类展示，按 GO 的三个分支（BP/MF/CC）分别展示最显著的富集条目
barplot(GO, showCategory = 5,,split="ONTOLOGY") +
  facet_grid(ONTOLOGY~., scale='free')
# showCategory = 5：每个分类下只展示最显著的5个GO条目
# split="ONTOLOGY"：按ONTOLOGY（即BP/MF/CC三个分支）对结果进行分组
# facet_grid(ONTOLOGY~.)：将分组后的结果分面展示，生成三个独立的子图

 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

KEGG <- enrichKEGG(gene = gene_entrez2$ENTREZID,
                   organism     = 'hsa', # 物种：hsa = Homo sapiens（人类）
                   pvalueCutoff = 0.05, 
                   qvalueCutoff =0.05) 
KEGG_res <- KEGG@result
write.table(KEGG_res,file="KEGG_res.txt",sep="\t",quote=F,row.names = F)
#简单作图
barplot(KEGG, showCategory = 20)
dotplot(KEGG, showCategory = 20)
#通路图可视化
#BiocManager::install('pathview')
library(pathview)
DEG <- rownames_to_column(DEG,var = 'SYMBOL')
pgenes <- inner_join(DEG,gene_entrez2,by = 'SYMBOL')
pgenes <- dplyr::select(pgenes,ENTREZID,P.Value)

pvalue <- -log10(pgenes$P.Value)
gene_names <- as.character(pgenes$ENTREZID)

gene_data <- pvalue
names(gene_data) <- gene_names
#选择兴趣的通路
hsa04978 <- pathview(gene.data  = gene_data,
                     pathway.id = "hsa04978", #通路名字
                     species    = "hsa", #代表人类
                     out.suffix = "pathview",  
                     limit      = list(gene= 10, cpd= 1)) #对基因数量限制，cpd化合物数量

 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

library(tidyverse)
library(limma)
setwd('C:/Users/30521/Desktop/R自学/GEO-analysis/富集分析/GSEA分析')#设置工作路径
load('GSE299988.rda')
rt <- rt3
exp <- exp3
str(exp) 
exp <- exp %>% mutate(across(everything(), as.numeric))
str(exp)
#确认行名列名一致#确认行名列名一致
identical(rownames(rt),colnames(exp))
#组间差异分析
group_list <- factor(rt$group,levels = c("Normal","Tumor"))
group_list
design <- model.matrix(~group_list)
#比较矩阵命名
design
#线性模型拟合
fit <- lmFit(exp,design)
#贝叶斯检验
fit2 <- eBayes(fit)
deg <- topTable(fit2, coef = 2, number = Inf)
DEG = na.omit(deg) 
DEG <- rownames_to_column(DEG,var = 'Gene')

# 至此其实和差异分析的处理相类似

#整理输入数据
library(enrichplot)
library(clusterProfiler)
library(org.Hs.eg.db)
#转换ENTREZID（如果在GSEA中下载的是ENTREZID则需要转换，下载的Symbol则不需要）
# genelist <- bitr(DEG$Gene, fromType="SYMBOL",
#                  toType="ENTREZID", OrgDb='org.Hs.eg.db')
# DEG <- inner_join(DEG,genelist,by=c("Gene"="SYMBOL"))
geneList <- DEG[,2]
names(geneList) <- as.character(DEG[,'Gene'])
#排序
geneList <- sort(geneList, decreasing = TRUE) # sort默认是从小到大排列
geneList

#读取gmt文件
gmt <- read.gmt('h.all.v2026.1.Hs.symbols.gmt')
GSEA <- GSEA(geneList,
             TERM2GENE = gmt, 
             pvalueCutoff = 0.05,
             pAdjustMethod = "BH",
             minGSSize = 20,
             maxGSSize = 500) 

res <- GSEA@result
write.csv(res,file = 'H_GSEA.csv')
# 选择一条通路进行展示
p <- gseaplot2(GSEA,geneSetID = 'HALLMARK_ALLOGRAFT_REJECTION',title = 'HALLMARK_ALLOGRAFT_REJECTION',color = 'red')
p
# 选择多条通路进行展示
p1 <- gseaplot2(GSEA,geneSetID = 1:3,pvalue_table = T, subplots = 1:2)
p1  

 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

library(tidyverse)
library(limma)
setwd('C:/Users/30521/Desktop/R自学/GEO-analysis/富集分析/GSEA分析')#设置工作路径
load('GSE299988.rda')
rt <- rt3
exp <- exp3
str(exp) 
exp <- exp %>% mutate(across(everything(), as.numeric))
str(exp)
#确认行名列名一致#确认行名列名一致
identical(rownames(rt),colnames(exp))
#提取 自定义基因,整理高低表达组
rt1 <- as.data.frame(t(exp['PTPRC',]))
# 选取的基因表达如果是大于这一组中的中位数，则认为是上调，完成单个基因的高表达或低表达区分
rt1$group <- ifelse(rt1$PTPRC > median(rt1$PTPRC),'up','down')
# 依照group进行排列，中间的'PTPRC'没用
rt1 <- arrange(rt1,'PTPRC',group)
identical(rownames(rt1),colnames(exp))
# 让exp的列名与rt1的行名一致
exp <- exp[,rownames(rt1)]
identical(rownames(rt1),colnames(exp))

# 组间差异分析
group_list <- factor(rt1$group,levels = c("down","up"))
group_list
design <- model.matrix(~group_list)
# 比较矩阵命名
design
# 线性模型拟合
fit <- lmFit(exp,design)
# 贝叶斯检验
fit2 <- eBayes(fit)
deg <- topTable(fit2, coef = 2, number = Inf)
DEG = na.omit(deg) 
DEG <- rownames_to_column(DEG,var = 'Gene')
library(enrichplot)
library(clusterProfiler)
library(org.Hs.eg.db)
# 转换ENTREZID
# genelist <- bitr(DEG$Gene, fromType="SYMBOL",
#                  toType="ENTREZID", OrgDb='org.Hs.eg.db')
# DEG <- inner_join(DEG,genelist,by=c("Gene"="SYMBOL"))
geneList <- DEG[,2]
names(geneList) <- as.character(DEG[,'Gene'])
# 排序
geneList <- sort(geneList, decreasing = TRUE)
geneList

# 读取gmt文件，去官网下载的C2基因集
gmt <- read.gmt('c2.cp.kegg_medicus.v2026.1.Hs.symbols.gmt')
GSEA <- GSEA(geneList,
             TERM2GENE = gmt,
             pvalueCutoff = 0.05, #改成1获得全部结果
             pAdjustMethod = "BH",
             minGSSize = 20,
             maxGSSize = 500) 
res <- GSEA@result
write.csv(res,file = 'KEGG_GSEA.csv')
p <- gseaplot2(GSEA,geneSetID = 'KEGG_MEDICUS_REFERENCE_COPI_VESICLE_FORMATION',
               title = 'KEGG_MEDICUS_REFERENCE_COPI_VESICLE_FORMATION',color = 'red')
p
p1 <- gseaplot2(GSEA,geneSetID = 1:3)
p1  

p2 <- gseaplot2(GSEA,
                geneSetID = c('KEGG_MEDICUS_REFERENCE_WNT_SIGNALING_PATHWAY',
                              'KEGG_MEDICUS_REFERENCE_ITGA_B_FAK_CDC42_SIGNALING_PATHWAY'))
p2

#3.是不是可以自创基因集？更改文件内容即可
library(enrichplot)
library(clusterProfiler)
library(org.Hs.eg.db)
gmt <- read.gmt('GOBP_MITOPHAGY.v2025.1.Hs.gmt')
# GO 生物学过程（GOBP）下的「线粒体自噬」功能模块

 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

rm(list = ls())
library(tidyverse)
setwd('C:/Users/30521/Desktop/R自学/GEO-analysis/VN图')

# 1.数据读取
DEG <- read.table("DEG2.txt",sep = "\t",check.names = F,stringsAsFactors = F,header = T,row.names = 1)
table(DEG$change)
diff <- rownames(DEG)[DEG$change!='NOT']

dat1 <- read.csv('ferroptosis.csv')
dat1 <- dat1$Gene.Symbol

dat2 <- read.csv('mitophage.csv')
dat2 <- dat2$Gene.Symbol

a <- list(DEGs=diff,
          ferroptosis=dat1,
          mitophage=dat2)

#演示数据，sample随机抽样（此处知道怎么选取颜色即可,pass）
library(RColorBrewer)
display.brewer.pal(9,'Set1')
brewer.pal(9, "Set1")[3]
"#4DAF4A"
display.brewer.pal(8,'Accent')

'#E41A1C'（红）  '#377EB8'（蓝）  '#4DAF4A'（绿）   '#984EA3'（紫）
'#FF7F00'（橙）  '#FFFF33'（黄）  '#A65628'（棕）   '#F781BF'（粉）
'#1B9E77'（深绿）'#D95F02'（深橙）'#7570B3'（深紫） '#E7298A'（深粉）
'#66A61E'（亮绿）'#E6AB02'（亮黄）'#A6761D'（棕黄） '#666666'（深灰）
'#A6CEE3'（浅蓝）'#1F78B4'（中蓝）'#B2DF8A'（浅绿） '#33A02C'（中绿）
'#FB9A99'（浅红）'#E31A1C'（中红）'#FDBF6F'（浅橙） '#FF7F00'（中橙）

set.seed(123) #随机种子（pass）
a <- list(a=sample(1:100, 35),
          b=sample(1:100, 45),
          c=sample(1:100, 55),
          d=sample(1:100, 80),
          e=sample(1:100, 90),
          f=sample(1:100, 100))

# 2.VN作图####
# 2.1 ggvenn包 最多4组
library(ggvenn)
ggvenn(a)

# 不太理解这一段美化为啥有问题，可能得看看具体ggvenn的用法
ggvenn(
  a,
  columns = NULL,
  show_elements = FALSE,
  show_percentage = TRUE,
  digits = 1,
  fill_color = c("blue", "yellow", "green", "red"),
  fill_alpha = 0.3,
  stroke_color = "black",
  stroke_alpha = 1,
  stroke_size = F,
  stroke_linetype = "solid",
  set_name_color = "black",
  set_name_size = 6,
  text_color = "black",
  text_size = 3,
  label_sep = ",",
  count_column = NULL,
  show_outside = c("auto", "none", "always"),
  auto_scale = FALSE
)

#2.2 venn包 最多7组  
# BiocManager::install('venn')
library(venn) 
venn(a)

{
venn(
  a,
  ilabels = "counts",  # 显示数字
  zcolor = c("#FF9999", "#66B2FF", "#99FF99"),  # 自定义柔和花瓣色（粉、蓝、绿）
  opacity = 0.6,       # 透明度
  box = FALSE,         # 是否除外边框
  ilcs = 1,            # 圈内数字大小
  sncs = 1.2,          # 组名大小
  border = "white",    # 边缘颜色
  lwd = 1.1            # 边框线条粗细
)
  }

# 3.输出重叠结果####
res <- Reduce(intersect, a)
res
write.csv(res,file = 'venn.csv',row.names = F)

 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

rm(list=ls())
setwd('C:/Users/30521/Desktop/R自学/GEO-analysis/箱线图') #设置工作路径
library(tidyverse)
library(ggplot2)
library(ggpubr)
load('GSE299988.rda')

# 数据准备
exp <- exp3
str(exp)
exp <- exp %>% mutate(across(everything(), as.numeric))
str(exp)
# 分组信息
rt <- rt3

# 差异表达可视化验证（箱线图）
# 将exp转置为矩阵，再弄成数据框，选取列名为上述韦恩图得到的交集基因
exps <- exp %>% t() %>% as.data.frame() %>% dplyr::select("PPARGC1A",'SNCA')
identical(rownames(exps),rownames(rt))
# [1] TRUE 

# 合并
rts <- cbind(exps,rt)

# 转换数据为长格式，方便绘图
rts_long <- rts %>% pivot_longer(cols = c("PPARGC1A",'SNCA'),
                                 names_to = "genes",values_to = "Expression")

#箱线图 gglot2
library(ggplot2)
p <- ggplot(data = rts_long,aes(x = genes,y = Expression,
                                fill = group))+
  geom_boxplot() +stat_compare_means(method ='wilcox.test') +theme_classic()+
  scale_fill_manual(values = c("#e67e22", "#3498db"))
p


# 小提琴图 ggpubr
library(ggplot2)
library(ggpubr)
p <- ggviolin(rts_long,x = "genes",y = "Expression",
              fill = "group",
              palette = c("#e67e22", "#3498db"),
              add.params = list(fill = "white")) +
  stat_compare_means(aes(group = group),
                     method = "wilcox.test", #method='wilcox.test' ;method = "t.test"
                     label = "p.signif",symnum.args = list(cutpoints = c(0, 0.001, 0.01, 0.05, 1),
                                                           symbols = c("***","**","*","ns")))
p


'#E41A1C'（红）  '#377EB8'（蓝）  '#4DAF4A'（绿）   '#984EA3'（紫）
'#FF7F00'（橙）  '#FFFF33'（黄）  '#A65628'（棕）   '#F781BF'（粉）
'#1B9E77'（深绿）'#D95F02'（深橙）'#7570B3'（深紫） '#E7298A'（深粉）
'#66A61E'（亮绿）'#E6AB02'（亮黄）'#A6761D'（棕黄） '#666666'（深灰）
'#A6CEE3'（浅蓝）'#1F78B4'（中蓝）'#B2DF8A'（浅绿） '#33A02C'（中绿）
'#FB9A99'（浅红）'#E31A1C'（中红）'#FDBF6F'（浅橙） '#FF7F00'（中橙）

mycol <- c('#984EA3','#B2DF8A')

mycol <- c("#e67e22", "#3498db")

{
p1 <- ggviolin(rts_long, x = "genes", y = "Expression", fill = "group", 
               palette = mycol, legend = "top", alpha = 0.5,
               add = "boxplot", 
               add.params = list(width = 0.3, #箱图的宽度
                                 alpha = 0.8, outlier.shape = NA),  # 调整箱线图宽度和透明度
               outlier.shape = NA, color = alpha("black", 0.5)) +             
  stat_compare_means(aes(group = group),
                     method = "wilcox.test",label = "p.signif",
                     symnum.args = list(cutpoints = c(0, 0.001, 0.01, 0.05, 1),
                                        symbols = c("***","**","*","ns"))) + 
  geom_jitter(aes(color = group),
              position = position_jitterdodge(
                jitter.width = 0.25,
                dodge.width = 0.8),   #散点的位置
              size = 1.2,
              alpha = 0.3) +
  scale_color_manual(values = mycol) +
  theme(panel.grid = element_line(color = NA),
        panel.grid.minor = element_line(color = NA),
        panel.border = element_rect(fill = NA, color = NA),
        axis.text.x = element_text(size = 10, face = "plain", hjust = 0.5),
        axis.title.x = element_blank(),
        axis.text.y = element_text(size = 10, face = "plain", hjust = 0.5),
        axis.title.y = element_text(vjust = 0.5, size = 10, face = "bold"))
}

p1