image.png

流程

image.png

数据分析可视化

image.pngimage.png

本地分析脚本

  1. rm(list = ls())
  2. options(stringsAsFactors = F)
  4. # 加载包
  5. library(edgeR)
  6. library(ggplot2)
  8. # 读取基因表达矩阵信息并查看分组信息和表达矩阵数据
  9. lname <- load(file = "data/Step01-airwayData.Rdata");lname
  11. # 表达谱
  12. filter_count[1:4,1:4]
  14. # 分组信息
  15. group_list <- group[match(colnames(filter_count),group$run_accession),2]
  16. group_list
  18. comp <- unlist(strsplit("Dex_vs_untreated",split = "_vs_"))
  19. group_list <- factor(group_list,levels = comp)
  20. group_list
  21. table(group_list)
  24. # 构建线性模型。0代表x线性模型的截距为0
  25. design <- model.matrix(~0+group_list)
  26. rownames(design) <- colnames(filter_count)
  27. colnames(design) <- levels(factor(group_list))
  28. design
  30. # 构建edgeR的DGEList对象
  31. DEG <- DGEList(counts=filter_count, 
  32.                group=factor(group_list))
  34. # 归一化基因表达分布
  35. DEG <- calcNormFactors(DEG)
  37. # 计算线性模型的参数
  38. DEG <- estimateGLMCommonDisp(DEG,design)
  39. DEG <- estimateGLMTrendedDisp(DEG, design)
  40. DEG <- estimateGLMTagwiseDisp(DEG, design)
  42. # 拟合线性模型
  43. fit <- glmFit(DEG, design)
  45. # 进行差异分析
  46. lrt <- glmLRT(fit, contrast=c(1,-1)) 
  48. # 提取过滤差异分析结果
  49. DEG_edgeR <- as.data.frame(topTags(lrt, n=nrow(DEG)))
  50. head(DEG_edgeR)
  52. # 筛选上下调，设定阈值
  53. fc_cutoff <- 1.5
  54. pvalue <- 0.05
  56. DEG_edgeR$regulated <- "normal"
  58. loc_up <- intersect(which( DEG_edgeR$logFC > log2(fc_cutoff) ),
  59.                     which( DEG_edgeR$PValue < pvalue) )
  61. loc_down <- intersect(which(DEG_edgeR$logFC < (-log2(fc_cutoff))),
  62.                       which(DEG_edgeR$PValue<pvalue))
  64. DEG_edgeR$regulated[loc_up] <- "up"
  65. DEG_edgeR$regulated[loc_down] <- "down"
  67. table(DEG_edgeR$regulated)
  68. ## 添加一列gene symbol
  69. # 方法1：使用包
  70. library(org.Hs.eg.db)
  71. keytypes(org.Hs.eg.db)
  73. library(clusterProfiler)
  74. id2symbol <- bitr(rownames(DEG_edgeR), 
  75.                   fromType = "ENSEMBL", 
  76.                   toType = "SYMBOL", 
  77.                   OrgDb = org.Hs.eg.db)
  78. head(id2symbol)
  80. DEG_edgeR <- cbind(GeneID=rownames(DEG_edgeR),DEG_edgeR)
  81. DEG_edgeR_symbol <- merge(id2symbol,DEG_edgeR,
  82.                           by.x="ENSEMBL",by.y="GeneID",all.y=T)
  83. head(DEG_edgeR_symbol)
  85. # 方法2：gtf文件中得到的id与name关系
  86. # Assembly: GRCh37(hg19) Release: ？
  87. # 使用上课测试得到的count做
  88. # 选择显著差异表达的结果
  89. library(tidyverse)
  90. DEG_edgeR_symbol_Sig <- filter(DEG_edgeR_symbol,regulated!="normal")
  92. # 保存
  93. write.csv(DEG_edgeR_symbol,"result/4.DEG_edgeR_all.csv", row.names = F)
  94. write.csv(DEG_edgeR_symbol_Sig,"result/4.DEG_edgeR_Sig.csv", row.names = F)
  95. save(DEG_edgeR_symbol,file = "data/Step03-edgeR_nrDEG.Rdata")
  99. ##====== 检查是否上下调设置错了
  100. # 挑选一个差异表达基因
  101. head(DEG_edgeR_symbol_Sig)
  103. exp <- c(t(express_cpm[match("ENSG00000001626",rownames(express_cpm)),]))
  104. test <- data.frame(value=exp, group=group_list)
  106. ggplot(data=test,aes(x=group,y=value,fill=group)) + geom_boxplot()

服务器分析传参脚本

image.png
原因：
1，可以在服务器上解决数据分析，实现环境统一
2，可以批量跑很多方案
脚本

## 帮助文档
library(getopt)
spec <- matrix(c(    
  'help',   'h',  0,  "logical",
  'count',  'c',  1,  "character", 
  'group',  'g',  1,  "character", 
  'comp',   'C',  1,  "character", 
  'fc',     'f',  1,  "numeric", 
  'pvalue', 'p',  1,  "numeric", 
  'od',     'o',  1,  "character"
), byrow = TRUE, ncol = 4)
opt <- getopt(spec)

print_usage <- function(spec=NULL){
  cat("Usage example: Rscript edgeR.R --count rawdata.txt --group group.txt --od ./
  Options: 
    --help  NULL
    --count   character  gene expression , raw count [forced]
    --group   character  group information [forced]
    --comp    character  A vs B, example: Dex_vs_untreated
    --fc      numeric    fold change cutoff, default: 1.5
    --pvalue  numeric    pvalue cutoff, default:0.05
    --od      character  outdir
  \n")
  q(status=1)
}
if (!is.null(opt$help) |is.null(opt$count)  |is.null(opt$od)) { print_usage(spec) }
if (is.null(opt$group) ) { print_usage(spec) }
if (is.null(opt$fc) )    { opt$fc <- 1.5 }
if (is.null(opt$pvalue) ) { opt$pvalue <- 0.05 }
if (!file.exists(opt$od)) { dir.create(opt$od) }

library(edgeR)
library(org.Hs.eg.db)
library(ggplot2)
library(tidyverse)

# 读取基因表达矩阵信息并查看分组信息和表达矩阵数据
# 读取基因表达矩阵
exprSet <- read.delim(opt$count,row.names = 1,sep = "\t",header = T,check.names = F)
# 检查表达谱
dim(exprSet)
exprSet[1:6,1:6]

group <- read.table(opt$group,header = T,sep = "\t",row.names = 1)
head(group)
group <- group[match(rownames(group),colnames(exprSet)),]

# 处理分组信息
comp <- unlist(strsplit(opt$comp,split = "_vs_"))
group_list <- factor(group,levels = comp)
group_list
table(group_list)


# 构建线性模型。0代表x线性模型的截距为0
design <- model.matrix(~0+group_list)
rownames(design) <- colnames(exprSet)
colnames(design) <- comp
design

# 构建edgeR的DGEList对象
DEG <- DGEList(counts=exprSet,  group=factor(group_list))

# 归一化基因表达分布
DEG <- calcNormFactors(DEG)

# 计算线性模型的参数
DEG <- estimateGLMCommonDisp(DEG,design)
DEG <- estimateGLMTrendedDisp(DEG, design)
DEG <- estimateGLMTagwiseDisp(DEG, design)

# 拟合线性模型
fit <- glmFit(DEG, design)

# 进行差异分析，1,-1意味着前比后
lrt <- glmLRT(fit, contrast=c(1,-1))

# 提取过滤差异分析结果
DEG_edgeR <- as.data.frame(topTags(lrt, n=nrow(DEG)))
head(DEG_edgeR)

# 筛选上下调，设定阈值
fc_cutoff <- opt$fc
pvalue <- opt$pvalue

DEG_edgeR$regulated <- "normal"

loc_up <- intersect(which(DEG_edgeR$logFC>log2(fc_cutoff)),
                    which(DEG_edgeR$PValue<pvalue))
loc_down <- intersect(which(DEG_edgeR$logFC < (-log2(fc_cutoff))),
                      which(DEG_edgeR$PValue<pvalue))

DEG_edgeR$regulated[loc_up] <- "up"
DEG_edgeR$regulated[loc_down] <- "down"

table(DEG_edgeR$regulated)

DEG_edgeR <- cbind(GeneID=rownames(DEG_edgeR),DEG_edgeR)
DEG_edgeR_Sig <- filter(DEG_edgeR,regulated!="normal")

# 保存
setwd(opt$od)
write.table(DEG_edgeR,paste0(opt$comp,"_edgeR_all.xls"), row.names = F, sep="\t",quote = F)
write.table(DEG_edgeR_Sig,paste0(opt$comp,"_edgeR_Sig.xls"), row.names = F, sep="\t",quote = F)