3 RNAseq差异分析

此版为基础版流程，如果未接触过这个版本，可以考虑直接跳过这个版本，直接学习新版流程。详细见RNAseq新版简化流程。

3.1 差异分析

3.1.1 加载包

rm(list = ls());gc()
library(LZ)
library(YZ)
library(tibble)
library(data.table)
library(DESeq2)
library(parallel)
library(BiocParallel)
library(ggplot2)
cat(" 您电脑线程为:", detectCores())
# 如果是12代以后的interCPU,建议最高不超过6或8。服务器可加大设置，
#  但不能大于线程总数。
# 此处如果电脑性能一般，建议直接使用n=4或者2。
n = 4
# register(MulticoreParam(n)) 苹果和linux电脑使用这句替代下句
register(SnowParam(n))

3.1.2 数据准备

RNAseq下游分析必须准备两个文件：表达矩阵表格文件，分组表格文件
将gene_count.csv，group.csv放在工作目录下

  • gene_count.csv 矩阵数据格式(数值型，整型)
ID sample1 sample2 sample3 sample4
gene1 34 23 56 23
gene2 35 23 12 23
gene3 12 78 78 78
  • group.csv 分组数据格式：需要组的行名 包含于 表达谱的列名 rownames(group) %in% colname(eset)
Sample Type BATCH
sample1 tumor 1
sample2 tumor 1
sample3 normal 2
sample4 normal 2
  • 注意：
    表达文件中的基因名是SYMBOL还是ESembleID。如为EsembleID,要注意是有小数点的ID还是没有小数点的。 有小数点的形式为这样：ESEM00000123.34，没有点的是ESEM00000123。
    还有记住基因名这列的列名，建议统一设置为ID。

3.1.3 文件夹准备

本文中有时也将文件夹称为目录，这两者等价。建议每个项目新建一个文件夹，例如本项目新建了一个名为 LZexample的文件夹，然后再在这个文件夹下建了一个data文件夹，以后data目录专门用来存原始文件，例 如RNAseq分析所需的eset.csv, group.csv或者更加原始的文件。

目录结构建议：本项目的目录初始结构，建议每个项目按着这个形式来。项目文件夹LZexample这个文 件夹名建议改成有意义的名称，一眼便能看出这个项目是什么数据或者什么目的，而data文件夹名不建议更改。 图片

3.1.4 差异分析预设置

# 设置工作目录，即之后所有的操作如果不指定文件夹，都将会在这个文件夹下进行
setwd("C:/data/LZexample")  # 按需更改成你的项目文件夹
getwd() # 检查是否更改工作目录成功了？
# Windows系统下默认的文件夹路径是 "C:/data/LZexample" 这中斜杆分隔文件夹，
# 如果是直接从win复制而来的，请将斜杠\改成反斜杠/, 就如下面设置的这样。（改成\\也行）

# 设置此次分析的标记
mark <- "T_C"  # 此次分析的标记(记录谁比谁或和筛选阈值，建议设置的有意义)
# 设置结果输出的文件夹，按需设定，可保持默认。
#  第一次分析可以不用改，但如是第二分析，必须至少修改mark,outdir其中一个，
#  否则会覆盖第一次的结果。
outdir <- "result"  
outdirsub <- paste0(outdir, "/",mark)
outdirsub.gsea <- paste0(outdirsub, "/gsea")
outdirsub.rich <- paste0(outdirsub, "/rich")
# 按以上设置，结果将会保存在当前工作目录下的result/T_C文件夹下

# 差异分析阈值设定
ffdr <- 0.1
fpval <- 0.05
flogfc <- 1

3.1.5 差异分析

# 读取并整理数据（如果都是按照上面的要求来的，不需要改这里的参数）
glist <- DEG_prepareData(eset_file = "data/eset.csv", #表达数据的相对路径
                         group_file = "data/group.csv", #分组文件的相对路径
                         id_dot = F,  # ESEM是否有点，有点设为T
                         col.by = "ID",  # 基因名列的列名
                         annot_trans = F, # 是否要注释，如果是EsembleID就需要设置为T
                         f_mark = mark)

# 差异分析 deseq2三部曲
dds <- DEG_DESeq2.dds(exprset.group=glist, batch = F)
DEG_DESeq2.pca(dds, outdir = outdirsub) # 此处有warning信息，不用管。
dds_list <- DEG_DESeq2.ana(dds)

# 差异后分析
# 构建GSEA官网软件分析所需格式文件 #此处会有warning,不用管
DEGres_ToGSEA(diffan.obj = dds_list, outdir = outdirsub.gsea) 
# all_father中记录了
#             差异分析的总表，默认阈值的差异基因表，
#             上调基因列表，下调基因列表
#             以及R-GSEA分析所需要的所有mRNA的表达排序列表。
#  是我们后续各种分析的万恶之源（因此命名all_father）
#  上述数据同时保存于当前工作目录/outdirsub.rich，
#  文件是一个多sheet的xlsx表
all_father <- DEGres_ToRICH(diffan.obj = dds_list, p=fpval, q=ffdr, 
                            f=flogfc, mark=mark, outdir = outdirsub.rich)
save.image(file = paste0(outdirsub, "/1.diff.img.RDATA")) # 保存中间数据

3.2 火山图

需要修改的是以下几个值：
df_valcano：文件读取时的文件路径
label_gene: 展示基因列表
不清楚的话按默认设置即可。仅修改label_gene: 展示基因列表即可。

# 本流程中不需要运行，后续想再次分析时可从此步开始
# load("./result/T_C/1.diff.img.RDATA") 
# library(LZ)
library(ggpubr);library(ggrepel);library(ggsci);library(scales)
library(tidyverse);library(dplyr);library(pheatmap);library(RColorBrewer)
# 导入火山图需要的数据，即差异分析后的未筛选表格(我们也称这个对象为resdf,
#  resdf文件涵盖差异分析的所有结果信息，可以做后续所有基于差异分析或者基因
#  列表的所有分析，如果后续分析时使用其它数据，请按这个resdf的格式改数据，
#  主要就是把数据的列名改成和resdf的列名相同，即可用此包的函数分析画图)
#  即df_valcano <- readxl::read_xlsx("xxx.xlsx", sheet = 1)
df_valcano <- all_father$DIFF.ALL
names(df_valcano) # 对应的列名必须为Gene, log2FC, PValue, FDR

# 模式设定 
# filterc参数设置
# "both": pvalue, padj均考虑模式
# "fdr":仅考虑fdr值模式, 
# "pvalue": 仅考虑p值模式
# 如果火山图出现横线位置“不对”时，
#   意味着fdr的考虑导致大量横线上基因被判定为非差异基因
#   请将filterc 设置为pvalue可以解决这个现象

# 火山图 无标记 （图片会自动保存）
DEGp_Volcano2(df_valcano, filterc = "both", 
              pvalue=0.05, fdr=0.1, logfc.p=1, logfc.n=-1,
              outdir="result", filename.base = "DEG_xx")

# 火山图 有标记 （图片会自动保存）
# 设定需要标记的marker gene
label_gene <- c('TFRC', 'ACSL1', 'LPCAT3', 'PCBP1', 'FTH1', 'SLC11A2',
                'SLC39A8', 'SAT1', 'FTL', 'GSS')
# 查看想展示的基因在不在差异分析总表中
# label_gene %in% df_valcano$Gene %>% all
# pic_data %>% filter(Row.names %in% label_gene)
DEGp_Volcano2(df_valcano, filterc = "both", 
              pvalue=0.05, fdr=0.1, logfc.p=1, logfc.n=-1,
              outdir="result", filename.base = "DEG_xx",
              label_geneset = NULL)

3.3 热图

3.3.1 热图表格数据预处理

library(LZ)
library(RColorBrewer)
library(pheatmap)
library(ggplot2)

resdf <- readxl::read_xlsx("D:/xx/DIFF.an_p0.05q0.2fc1.xlsx", 
                           sheet = 1)
sig.df <- resdf %>% 
  dplyr::filter(abs(log2FC) > 1, PValue < 0.05, FDR < 0.1) %>% 
  dplyr::arrange(log2FC) %>% data.frame(check.names = F)
rownames(sig.df) <- sig.df$Gene

# 画图的表达数据
sig.expr <- sig.df[, 8:(ncol(sig.df)-1)]
sig.expr %>% range()
# 画图的显示基因（这里表示最大差异前10和后10,
#   可以将括号里的两个10改为100，就是显示最大差异前100后100）
show.gene <- c(head(rownames(sig.df), 10), tail(rownames(sig.df), 10))

3.3.2 热图分组数据预处理

R中为一个行名为样本名的一列dataframe, excel中表现为两列，第1列为样本名，第2列为分组。
新手建议excel里操作然后保存为csv格式，然后使用groupdf <- read.csv("your csv", row.names = 1)读入。

rowname Tpye
S1-11A 11A
S2-11A 11A
S3-11A 11A
R1-01A 01A
R2-01A 01A
R3-01A 01A 
# 方法一，新手不建议用这种方法
groupdf <- data.frame(
  row.names = colnames(sig.expr),
  # 要一一对应才行
  Type = c(rep("11A", 3), rep("01A", 3)),
  check.names = F)
groupdf[,1] <- as.factor(groupdf[,1])
# 方法二，新手推荐这种方法，先在excel中处理好数据
groupdf <- read.csv("your csv", row.names = 1)
groupdf[,1] <- as.factor(groupdf[,1])

3.3.3 选定一种算法作图

# 检查groupdf的行名是否全部与sig.expr的列名一一对应
# 这句如果输出为False，不可继续进行画图操作，要重新调整数据。
all(rownames(groupdf) == colnames(sig.expr))
# 作图，默认算法
DEGp_Pheat(
  sig.expr,
  gene = show.gene, 
  pheno = groupdf,
  outdir = "D:/xx/yy",
  f_mark = "heatmap",
  nlevel = 2,
  pic_w = 4.5,
  pic_h = 4.5,
  c_w = 16,
  c_h = 9.8, #7.4 # 9.9
  c_colors = c("blue", "#FFFFCC", "red"),
  color.num = 100,
  fontsize = 10,
  fontsize.col = 12,
  angle = "90",
  cluster_rows = T,
  cluster_cols = T
)

3.3.4 循环各种聚类算法作图

m <- c("ward.D", "ward.D2", "single", "complete", 
"average", "mcquitty", 
"median", "centroid")
d <- c("correlation", "euclidean", "maximum", "manhattan",
       "canberra", "binary", "minkowski")
for (i in m) {
  for (i.d in d) {
    cat(i, i.d, "\n")
    DEGp_Pheat(
      sig.expr,
      gene = show.gene, 
      pheno = groupdf,
      outdir = "D:/xx/yy",
      f_mark = paste0("heatmap.m.", i, "_", i.d),
      nlevel = 2,
      pic_w = 5,
      pic_h = 8,
      c_w = 18,
      c_h = 10,
      f_z = 10,
      f_z.col = 14,
      angle = "90",
      cluster_rows = T,
      cluster_cols = T,
      clustering_method = i,
      clustering_distance_rows=i.d
    )  
  }
}

3.4 富集分析

3.4.1 GO & KEGG 富集分析

3.4.1.1 一键脚本(批量处理)

这是一个一键脚本，请新建一个单独的文件写这段脚本，然后按这个脚本的顶部注释 修改resdf outd fc.list处即可，运行即可批量出不同FC的富集分析结果。

# 此脚本为GO、KEGG分析（需要一个输入文件即可，为差异分析流程后的resdf文件）
# 即为第一步(或1脚本)的结果的一个结果文件（DIFF_an_***.xlsx）
# 即为resdf文件，此文件是差异分析后的总表
# 注意如果采用了其他的分析方法得到差异分析后表，运行这个脚本时可能需要更改
# 列名即我们的resdf对象的列名为Gene, log2FC,PValue,FDR，需要与这些个列名保
# 持一致。
# 此脚本中的需要修改的位于 /// *** /// 行中，另外还有一个LZ::setproxy()行，
#   如果没有代理工具，或者代理工具不支持http代理，或者端口不通，请不要运行。
rm(list = ls());gc() # 清空所有对象，慎用，必要时用
suppressMessages({ suppressWarnings({
  library(LZ);library(YZ)
  library(tidyverse);library(data.table)
  library(clusterProfiler);library(enrichplot)
  library(topGO);library(Rgraphviz)
  library(RColorBrewer);library(ggsci);library(pheatmap)
  library(readxl)
}) })
# 若无代理工具，切勿运行 
# LZ::setproxy() # 高危！！！新手不要运行此行，会使当前窗口断网！！！
# Sys.getenv('http_proxy') Sys.setenv('http_proxy'='') Sys.setenv('https_proxy'='')

# 如果是自己提供的表格，要按需修改列名为标准的resdf格式的列名：
# 即: 表格必须含列名Gene, log2FC, PValue, FDR这四列，列名必须为这四个，
# 提前在xlsx中修改好，然后取消下面这句的注释符，运行
# resdf <- readxl::read_xlsx("result/rnaseqOR-NC/rich/DIFF.an_OR-NC.xlsx",
#                  sheet = 1) %>% as.data.frame()

# 按流程跑下来是运行这句，注意这句和上面的注释掉的是二选一，不要重复运行
resdf <- all_father$DIFF.ALL
# 输出目录
outd = "result/xx/rich" 
# logFC阈值, 多个阈值的话，
#   写成fc.list <- list('1.2' = log2(1.2), '2' = log2(2))
# 注意！！！！！！：括号里log2(2)的2，和引号里'2'的2都要需同步要改。！！！
# 否则可能会覆盖结果 
# logFC阈值, 多个阈值
fc.list <- list('1.5' = log2(1.5), '2' = log2(2), '4' = log2(4))
# logFC阈值, 单个阈值运行这句，也是二选一，不要重复运行
# fc.list <- list('2'=log2(2))

# 预处理数据符合GOKEGG分析的要求
# # 不同fc条件下的GOgenelist list(ALL, UP, DOWN)
gogenelist <- lapply(fc.list, function(x) {
  DEG_prepareGOglist(resdf, logfc = x, p = 0.05, fdr = 0.1) })
# gogenelist %>% length()
# 对logFC迭代，每个FC新建一个目录，用来存upgene, downgene, allgene的GO结果
enrich <- DEG_runENRICH(genelist = gogenelist, outdir = outd, 
                        glist.save = T, rungo = T, runkegg = T, rapid = T)

3.4.1.2 简易GO,KEGG一次分析

如果已经得到了差异基因列表，且无需批量分析，可以进行这个简易分析。
数据格式: head(genelist.lh)
[1] “AARS1” “AATF” “ABCB7” “ABCE1” “ABHD11” “ABHD12”

# 简易GO,KEGG一次分析(即：已经得到了差异基因列表)
# LZ::setproxy() # 代理设置，新手别碰，会断网
# 差异基因列表
genelist.lh <- pic.list$sig.data$Gene 
# 转换ID
gene_df <- bitr(genelist.lh, fromType = "SYMBOL", toType = c("ENTREZID", "UNIPROT"), 
                OrgDb = 'org.Hs.eg.db')
# GO分析
go.lh <- DEG_GO(gene_df, orgdb = "org.Hs.eg.db", sigNodes = 20, 
                resultdir="./result/proteinOR-NC", filemark = "p1.5_g_2")
go.lhdf <- sapply(go.lh, function(x) x@result, simplify = T)
write_xlsx(go.lhdf, path = "./result/xx/lh_go.all.xlsx")
# KEGG分析
kegg.lh <- DEG_KEGG(gene_df)
write_xlsx(kegg.lh$pSigDF, path = "./result/xx/lh_kegg.all.xlsx")

3.4.1.3 GO、KEGG分析结果可视化

# dotplot go
# 读取go分析保存的表格
#dotData <- go$GODF$"倍数"$变化趋势(BP)
# 自己提供表格读取
#dotData <- readxl::read_xlsx("kegg.xlsx", sheet = 1) 
dotData <- enrich$GODF$"2"$all
# 筛选数据（按需配合其他筛选）
pic.dot <- DEGp_Dotplot2(dotData, head = 20, title = 'TOP of GO', 
                         resultdir = "./result/proteinOR-NC", 
                         filemark = 'GO_top', pcompress = 2,
                         pic.save = T)
# dotplot kegg
# 读取kegg分析保存的表格，格式要求
# 必须要有这5列Description, GeneRatio, pvalue, qvalue, Count列。
# 读取表格
# dotDatak <- readxl::read_xlsx("./result/proinOR/kegg.xlsx", sheet = 1)

dotDatak <- enrich$KEGGDF$'2'$up
# 筛选数据（按需配合其他筛选）
pic.dotk <- DEGp_Dotplot2(dotDatak, head = 20, title = 'TOP of KEGGpathway', 
                          resultdir = "./result/proteinOR-NC", 
                          filemark = 'KEGG_top', pcompress = 1,
                          pic.save = T)

# 组合图(可选运行，比例不是很好调整，单独出图AI内调整更自由)
gh <- ggplotGrob(pic.dot)
gd <- ggplotGrob(pic.dotk)
cowplot::plot_grid(gh, gd, rel_widths = c(1, 1.2))
ggsave(paste0(dir_out, "/GO_KEGG_top.pdf"), width = 16, height = 10)

3.4.2 GSEA分析

3.4.2.1 R GSEA 批量分析

  • GSEA官网提供了GSEA分析软件和MSigDB数据库中的所有通路下载，如果需要更多的通路集可以自行下载，本包内置了MSigDB的H,C1-8的所有大类集及C2,C5的部分重要子类集，还有整理好的最新版的KEGG官方的PATHWAY通路集合。

  • 本包构建了一个图形界面函数runAPP_GSEA()，可在安装完LZ包后直接通过运行LZ::runAPP_GSEA()启动图形界面，也可在浏览器中打开。详细见R GSEA 图形界面

library(LZ)
library(clusterProfiler)
library(enrichplot)
library(shiny)
library(ggplot2)

# 1. Gene list 排序表
genelist <- all_father$gsealist

# 2. Pathway Gene Set 表
# 内置数据集 gmt.largelist.23.12.Hs.symbols
#   含1. msigDB数据库的全部通路大类 
#     2. msigDB的C2,C5的部分子集通路[这些子集是C2,C5的一部分] 
#     3.最新版本的KEGG全部通路
data("gmt.largelist.23.12.Hs.symbols")   
# 选择KEGG通路集合，把美元符号后面字符删掉，然后按tab键可以选择其他数据集
gmt <- gmt.largelist.23.12.Hs.symbols$kegg.all.23.12.Hs.symbols.gmt

# 全部该GeneSet数据的通路GSEA分析  -----------
gsea <- DEG_runGSEA(genelist=genelist, gmt_set=gmt, pic.save=F)
# 将所有的分析结果导出到本地[gsea.result]文件夹，统计总结表名为gsea_stat
#   导出文件夹名和文件名均可按需修改
DEGp_GSEA_plotALL(gsea, result_dir = "gsea.result", 
                  xl_filename = "gsea_stat")

# 单个GeneSet数据的通路GSEA分析（且从自己准备的gmt文件开始），可从MSigDB网
#  站搜索下载
# 设定文件路径
gmt_filename <- "D:/Team/RNAseq/data/geneset/WP_FERROPTOSIS.v2022.1.Hs.gmt"
gmt_single <- clusterProfiler::read.gmt(gmt_filename)
gsea.single <- DEG_runGSEA(genelist = genelist, gmt_set = gmt_single, 
                           pic.save=T, outdir = "./gsea.result2/", 
                           filename = "ferr")

3.4.2.2 R GSEA 图形界面

图形界面可以用自己的表格数据(上传表格不是必须，方式二启动就不需要)上传来做GSEA分析，表格必须有且仅有两列，分别为Gene列和log2FC列，具体表格形式如下：

Gene log2FC
geneA 9.8
geneX 4.3
geneY 1.2
… …
geneZ 0.3
geneB -0.8
geneZ -5.2

注意：虽然要求第二列名为log2FC,但第二列只要是表示变化倍数就可以，不一定要是log2后数据，没有排序也没关系。

# 启动方法1（不带参数启动），适合自己已经有了差异分析结果的表格的情况，
#  那么运行此句后，不需要在R里写任何代码，如果通过LZ::runAPP_GSEA()甚至
#  都不需要加载包。
#  使用这种方法，则必须要上传上述指定的表格形式的表格后才能点击画图。
runAPP_GSEA() # 或在安装成功LZ包后，直接通过LZ::runAPP_GSEA()来启动

# 启动方法2（带参数启动），适合在R中已经有genelist的情况。
# 已在R里有了genelist的话，运行此句启动图形界面，这样就不需要上传表格
runAPP_GSEA(genelist = genelist)

3.4.2.3 一些进阶操作及技巧(不会没有关系，不做过多解释，自行体会，无R语言基础者慎入)

# 进阶操作 ---------------------------
# 1. 查找指定通路的图
#  如果自己知道通路的名字，可以通过查找来定位到通路的位置，然后单独画图
pathway = "^ABC"
n = grep(pathway, gsea[,"ID"])
grep(pathway, gsea[,"ID"], value = T) # 查看找到的通路名称，必须时唯一值，否则请使用更加详细的查找条件再查找一次。
DEGp_GSEA(gsea, num = n)

# 2. 从当前的gmt文件中获取指定的Pathway Gene set(gmt对象名一定腰围gmt才行)
gmtdf.find <- find_pathway("^Ferr")
gsea.single <- DEG_runGSEA(genelist = genelist, gmt_set=gmtdf.find, 
                           pic.save=T, outdir = "./gsea.ytb3/", 
                           filename = "fer_taget")

# 3. 如果自己的分析中想将自己的数据转化为GSEA要求的genelist数据，
#    例如自己的分析项目中有一个名为re的对象，该对象中有
#    基因名列XX 和 基因变化倍数列 YY
genelist <- re$YY
names(genelist) <- re$XX
genelist <- na.omit(sort(genelist, decreasing = T))

# 4. GSEA图形界面中若使用结束按钮关闭程序，会将最后一次画图的数据保留在
#    R会话中，具体如下：
# 最后画的一幅图
pic_gsea
# 保存图片
pdf('aaa.pdf', width = 6, height = 5)
pic_gsea
dev.off()
# 最后一次选择的通路集详细
sy_gmt_taget
# 最后一次选择的通路名称
sy_pathway_name
# 最后一次选择的大类通路集详细的前六行
sy_gmt %>% head()

# 5. 转换成宽型文件
gmt.w2 <- gmt_longTowide2(gmt)

3.4.2.4 GSEA官方软件

自行搜素方法，网络上有大量图文教程 …

3.4.3 转录组其它可视化

On the way …