数据标准化—异常样本于重复性检测—差异表达分析

一、数据标准化

  1. rm(list = ls())
  2. options(stringsAsFactors = F)
  3. library(stringr)
  5. ## ====================1.读取数据
  6. # 读取raw count表达矩阵
  7. rawcount <- read.table("data/raw_counts.txt",row.names = 1, 
  8.                        sep = "\t", header = T)
  9. colnames(rawcount)
  11. # 查看表达谱
  12. rawcount[1:4,1:4]
  14. # 去除前的基因表达矩阵情况
  15. dim(rawcount)
  17. # 获取分组信息
  18. group <- read.table("data/filereport_read_run_PRJNA229998_tsv.txt",
  19.                     header = T,sep = "\t", quote = "\"")
  20. colnames(group)
  22. # 提取表达矩阵对应的样本表型信息
  23. group <- group[match(colnames(rawcount), group$run_accession),
  24.                c("run_accession","sample_title")]
  25. group
  27. # 差异分析方案为：Dex vs untreated
  28. group$sample_title <- str_split_fixed(group$sample_title,pattern  = "_", n=2)[,2]
  29. group
  30. write.table(group,file = "data/group.txt",row.names = F,sep = "\t",quote = F)
  33. ## =================== 2.表达矩阵预处理
  34. # 过滤低表达基因
  35. keep <- rowSums(rawcount>0) >= floor(0.75*ncol(rawcount))
  36. table(keep)
  38. filter_count <- rawcount[keep,]
  39. filter_count[1:4,1:4]
  40. dim(filter_count)
  42. # 加载edgeR包计算counts per millio(cpm) 表达矩阵
  43. library(edgeR)
  44. express_cpm <- cpm(filter_count)
  45. express_cpm[1:6,1:6]
  47. # 保存表达矩阵和分组结果
  48. save(filter_count, express_cpm, group, file = "data/Step01-airwayData.Rdata")

二、异常样本与重复性检测

  1. rm(list = ls())
  2. options(stringsAsFactors = F)
  4. # 加载包，设置绘图参数
  5. library(ggplot2)
  6. library(ggsci)
  7. mythe <- theme_bw() + theme(panel.grid.major=element_blank(),
  8.                             panel.grid.minor=element_blank())
  11. # 加载原始表达的数据
  12. lname <- load(file = "data/Step01-airwayData.Rdata")
  13. lname
  14. exprSet <- log10(as.matrix(express_cpm)+1)
  15. exprSet[1:6,1:6]
  19. ## 1.样本表达总体分布-箱式图
  20. # 构造绘图数据
  21. data <- data.frame(expression=c(exprSet),
  22.                    sample=rep(colnames(exprSet),each=nrow(exprSet)))
  23. head(data)
  25. p <- ggplot(data = data, aes(x=sample,y=expression,fill=sample))
  26. p1 <- p + geom_boxplot() + 
  27.   mythe+ theme(axis.text.x = element_text(angle = 90)) + 
  28.   xlab(NULL) + ylab("log10(CPM+1)") + scale_fill_lancet()
  30. p1
  32. # 保存图片
  33. png(file = "result/1.sample_boxplot.png",width = 800, height = 900,res=150)
  34. print(p1)
  35. dev.off()
  39. ## 2.样本表达总体分布-小提琴图
  40. p2 <- p + geom_violin() +  mythe +
  41.   theme(axis.text = element_text(size = 12),
  42.         axis.text.x = element_text(angle = 90)) + 
  43.   xlab(NULL) + ylab("log10(CPM+1)")+scale_fill_lancet()
  44. p2
  46. # 保存图片
  47. png(file = "result/1.sample_violin.png",width = 800, height = 900,res=150)
  48. print(p2)
  49. dev.off()
  53. ## 3.样本表达总体分布-概率密度分布图
  54. m <- ggplot(data=data, aes(x=expression)) 
  55. p3 <- m +  geom_density(aes(fill=sample, colour=sample),alpha = 0.1) + 
  56.   xlab("log10(CPM+1)") + mythe +scale_fill_lancet()
  57. p3
  59. # 保存图片
  60. png(file = "result/1.sample_density.png",width = 800, height = 700, res=150)
  61. print(p3)
  62. dev.off()
  65. # 魔幻操作，一键清空
  66. rm(list = ls()) 
  67. options(stringsAsFactors = F)
  69. library(FactoMineR)
  70. library(factoextra)
  71. library(corrplot)
  72. library(pheatmap)
  73. library(tidyverse)
  75. # 加载数据并检查
  76. lname <- load(file = 'data/Step01-airwayData.Rdata')
  77. lname
  81. ## 1.样本之间的相关性-层次聚类树----
  82. dat <- log10(express_cpm+1)
  83. dat[1:4,1:4]
  84. dim(dat)
  85. sampleTree <- hclust(dist(t(dat)), method = "average")
  86. plot(sampleTree)
  88. # 提取样本聚类信息
  89. temp <- as.data.frame(cutree(sampleTree,k = 2)) %>% 
  90.   rownames_to_column(var="sample")
  91. temp1 <- merge(temp,group,by.x = "sample",by.y="run_accession")
  92. table(temp1$`cutree(sampleTree, k = 2)`,temp1$sample_title)
  94. # 保存结果
  95. pdf(file = "result/2.sample_Treeplot.pdf",width = 7,height = 6)
  96. plot(sampleTree)
  97. dev.off()
  101. ## 2.样本之间的相关性-PCA----
  102. # 第一步，数据预处理
  103. dat <- log10(express_cpm+1)
  104. dat[1:4,1:4]
  106. dat <- as.data.frame(t(dat))
  107. dat_pca <- PCA(dat, graph = FALSE)
  109. group_list <- group[match(group$run_accession,rownames(dat)), 2]
  110. group_list
  112. # geom.ind: point显示点，text显示文字
  113. # palette: 用不同颜色表示分组
  114. # addEllipses: 是否圈起来
  115. mythe <- theme_bw() + 
  116.   theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank()) +
  117.   theme(plot.title = element_text(hjust = 0.5))
  119. p <- fviz_pca_ind(dat_pca,
  120.                   geom.ind = "text",  #point
  121.                   col.ind = group_list, 
  122.                   palette = c("#00AFBB", "#E7B800"), 
  123.                   addEllipses = T,  
  124.                   legend.title = "Groups") + mythe
  125. p
  127. # 保存结果
  128. pdf(file = "result/2.sample_PCA.pdf",width = 6.5,height = 6)
  129. plot(p)
  130. dev.off()
  135. ## 3.样本之间的相关性-cor----
  136. # 选择差异变化大的基因算样本相关性
  137. exprSet <- express_cpm
  138. exprSet = exprSet[names(sort(apply(exprSet, 1, mad),decreasing = T)[1:800]),]
  139. dim(exprSet)
  141. # 计算相关性
  142. M <- cor(exprSet,method = "spearman")
  143. M
  145. # 构造注释条
  146. anno <- data.frame(group=group$sample_title,row.names = group$run_accession )
  148. # 保存结果
  149. pheatmap(M,display_numbers = T,
  150.          annotation_col = anno,
  151.          fontsize = 10,cellheight = 30,
  152.          cellwidth = 30,cluster_rows = T,
  153.          cluster_cols = T,
  154.          filename = "result/2.sample_Cor.pdf",width = 7.5,height = 7)

三、差异表达分析

  1. rm(list = ls())
  2. options(stringsAsFactors = F)
  4. # 加载包
  5. library(edgeR)
  6. library(ggplot2)
  8. # 读取基因表达矩阵信息并查看分组信息和表达矩阵数据
  9. lname <- load(file = "data/Step01-airwayData.Rdata")
  10. lname
  12. # 表达谱
  13. filter_count[1:4,1:4]
  15. # 分组信息
  16. group_list <- group[match(colnames(filter_count),group$run_accession),2]
  17. group_list
  19. # treat vs control
  20. comp <- unlist(strsplit("Dex_vs_untreated",split = "_vs_"))
  21. group_list <- factor(group_list,levels = comp)
  22. group_list
  23. table(group_list)
  26. # 构建线性模型。0代表x线性模型的截距为0
  27. design <- model.matrix(~0+group_list)
  28. rownames(design) <- colnames(filter_count)
  29. colnames(design) <- levels(factor(group_list))
  30. design
  32. # 构建edgeR的DGEList对象
  33. DEG <- DGEList(counts=filter_count, 
  34.                group=factor(group_list))
  36. # 归一化基因表达分布
  37. DEG <- calcNormFactors(DEG)
  39. # 计算线性模型的参数
  40. DEG <- estimateGLMCommonDisp(DEG,design)
  41. DEG <- estimateGLMTrendedDisp(DEG, design)
  42. DEG <- estimateGLMTagwiseDisp(DEG, design)
  44. # 拟合线性模型
  45. fit <- glmFit(DEG, design)
  47. # 进行差异分析
  48. lrt <- glmLRT(fit, contrast=c(1,-1)) 
  50. # 提取过滤差异分析结果
  51. DEG_edgeR <- as.data.frame(topTags(lrt, n=nrow(DEG)))
  52. head(DEG_edgeR)
  54. # 筛选上下调，设定阈值
  55. fc_cutoff <- 1.5
  56. pvalue <- 0.05
  58. DEG_edgeR$regulated <- "normal"
  60. loc_up <- intersect(which( DEG_edgeR$logFC > log2(fc_cutoff) ),
  61.                     which( DEG_edgeR$PValue < pvalue) )
  63. loc_down <- intersect(which(DEG_edgeR$logFC < (-log2(fc_cutoff))),
  64.                       which(DEG_edgeR$PValue<pvalue))
  66. DEG_edgeR$regulated[loc_up] <- "up"
  67. DEG_edgeR$regulated[loc_down] <- "down"
  69. table(DEG_edgeR$regulated)
  72. ## 添加一列gene symbol
  73. # 方法1：使用包
  74. library(org.Hs.eg.db)
  75. keytypes(org.Hs.eg.db)
  77. library(clusterProfiler)
  78. id2symbol <- bitr(rownames(DEG_edgeR), 
  79.                   fromType = "ENSEMBL", 
  80.                   toType = "SYMBOL", 
  81.                   OrgDb = org.Hs.eg.db)
  82. head(id2symbol)
  84. DEG_edgeR <- cbind(GeneID=rownames(DEG_edgeR),DEG_edgeR)
  85. DEG_edgeR_symbol <- merge(id2symbol,DEG_edgeR,
  86.                           by.x="ENSEMBL",by.y="GeneID",all.y=T)
  87. head(DEG_edgeR_symbol)
  90. # 方法2：gtf文件中得到的id与name关系
  91. # Assembly: GRCh37(hg19) Release: ？
  92. # 使用上课测试得到的count做
  96. # 选择显著差异表达的结果
  97. library(tidyverse)
  98. DEG_edgeR_symbol_Sig <- filter(DEG_edgeR_symbol,regulated!="normal")
  100. # 保存
  101. write.csv(DEG_edgeR_symbol,"result/4.DEG_edgeR_all.csv", row.names = F)
  102. write.csv(DEG_edgeR_symbol_Sig,"result/4.DEG_edgeR_Sig.csv", row.names = F)
  103. save(DEG_edgeR_symbol,file = "data/Step03-edgeR_nrDEG.Rdata")
  107. ##====== 检查是否上下调设置错了
  108. # 挑选一个差异表达基因
  109. head(DEG_edgeR_symbol_Sig)
  111. exp <- c(t(express_cpm[match("ENSG00000001626",rownames(express_cpm)),]))
  112. test <- data.frame(value=exp, group=group_list)
  114. ggplot(data=test,aes(x=group,y=value,fill=group)) + geom_boxplot()
  117. ##可视化
  118. rm(list = ls())
  119. options(stringsAsFactors = F)
  120. library(ggplot2)
  121. library(tidyverse)
  123. # 读差异分析结果
  124. lname <- load(file = "data//Step03-edgeR_nrDEG.Rdata")
  126. # 根据需要修改DEG的值
  127. data <- DEG_edgeR_symbol
  128. colnames(data)
  131. # 绘制火山图
  132. colnames(data)
  133. p <- ggplot(data=data, aes(x=logFC, y=-log10(PValue),color=regulated)) + 
  134.   geom_point(alpha=0.5, size=1.2) + 
  135.   theme_set(theme_set(theme_bw(base_size=20))) + theme_bw() +
  136.   theme(panel.grid.major=element_blank(),panel.grid.minor=element_blank()) +
  137.   xlab("log2FC") + ylab("-log10(Pvalue)") +
  138.   scale_colour_manual(values = c(down='blue',normal='grey',up='red')) +
  139.   geom_vline(xintercept=c(-(log2(1.5)),log2(1.5)),lty=2,col="black",lwd=0.6) +
  140.   geom_hline(yintercept = -log10(0.05),lty=2,col="black",lwd=0.6)
  141. p
  144. # 添加top基因
  145. # 通过FC选取TOP10
  146. label <- data[order(abs(data$logFC),decreasing = T)[1:10],]
  147. # 通过pvalue选取TOP10
  148. #label <- data[order(abs(data$PValue),decreasing = F)[1:10],]
  149. label <- na.omit(label)
  151. p1 <- p + geom_point(size = 3, shape = 1, data = label) +
  152.   ggrepel::geom_text_repel( aes(label = SYMBOL), data = label, color="black" )
  154. p1
  157. # 保存结果
  158. png(file = "result/5.Volcano_Plot.png",width = 900, height = 800, res=150)
  159. plot(p1)
  160. dev.off()
  162. rm(list = ls())
  163. options(stringsAsFactors = F)
  164. # 加载包
  165. library(pheatmap)
  166. library(tidyverse)
  168. # 加载原始表达矩阵
  169. lname <- load(file = "data/Step01-airwayData.Rdata")
  170. lname
  172. express_cpm1 <- rownames_to_column(as.data.frame(express_cpm) ,var = "ID")
  174. # 读取差异分析结果
  175. lname <- load(file = "data/Step03-edgeR_nrDEG.Rdata")
  176. lname
  178. # 提取所有差异表达的基因名
  179. edgeR_sigGene <- DEG_edgeR_symbol[DEG_edgeR_symbol$regulated!="normal",]
  180. head(edgeR_sigGene)
  182. data <- merge(edgeR_sigGene,express_cpm1,by.x = "ENSEMBL",by.y = "ID")
  183. data <- na.omit(data)
  184. data <- data[!duplicated(data$SYMBOL),]
  186. # 绘制热图
  187. dat <- select(data,starts_with("SRR"))
  188. rownames(dat) <- data$SYMBOL
  189. dat[1:4,1:4]
  190. anno <- data.frame(group=group$sample_title,row.names = group$run_accession)
  192. pheatmap(dat,scale = "row",show_colnames =T,
  193.          show_rownames = F, cluster_cols = T,
  194.          annotation_col=anno,
  195.          main = "edgeR's DEG")
  198. # 显示指定symbol,这里随便展示10个基因symbol
  199. labels <- rep(x = "",times=nrow(dat))
  200. labels[1:10] <- rownames(dat)[1:10]
  201. pheatmap(dat,scale = "row",show_colnames =T,
  202.          show_rownames = T, 
  203.          cluster_cols = T,
  204.          annotation_col=anno,
  205.          labels_row = labels,
  206.          fontsize_row = 8,
  207.          main = "edgeR's DEG")
  211. # 按照指定顺序绘制热图
  212. dex_exp <- express_cpm[,match(rownames(anno)[which(anno$group=="Dex")],
  213.                               colnames(express_cpm))]
  215. untreated_exp <- express_cpm[,match(rownames(anno)[which(anno$group=="untreated")],
  216.                               colnames(express_cpm))]
  218. data_new <- cbind(dex_exp, untreated_exp)
  219. dat1 <- data_new[match(edgeR_sigGene$ENSEMBL,rownames(data_new)),]
  221. pheatmap(dat1, scale = "row",show_colnames =T,show_rownames = F, 
  222.               cluster_cols = F, 
  223.               annotation_col=anno,
  224.               main = "edgeR's DEG")

四、功能注释

学习KEGG

五、功能富集

1、GO_KEGG

rm(list = ls())
options(stringsAsFactors = F)

library(clusterProfiler)
library(org.Hs.eg.db)
library(GSEABase)
library(ggplot2)
library(tidyverse)


## Error in download.KEGG.Path(species)
# https://github.com/YuLab-SMU/clusterProfiler/pull/471
getOption("clusterProfiler.download.method")
#R.utils::setOption("clusterProfiler.download.method",'auto')
options(clusterProfiler.download.method = "wininet")
#options(clusterProfiler.download.method = "wget")
getOption("clusterProfiler.download.method")



# 读取差异分析结果
load(file = "data/Step03-edgeR_nrDEG.Rdata")
ls()

# 提取所有差异表达的基因名
DEG <- DEG_edgeR_symbol[DEG_edgeR_symbol$regulated!="normal",2]
head(DEG)

## ===GO数据库, 输出所有结果，后续可根据pvalue挑选结果
ego_CC <- enrichGO(gene=DEG, OrgDb= 'org.Hs.eg.db', keyType='SYMBOL', ont="CC", pvalueCutoff= 1,qvalueCutoff= 1)
ego_MF <- enrichGO(gene=DEG, OrgDb= 'org.Hs.eg.db', keyType='SYMBOL', ont="MF", pvalueCutoff= 1,qvalueCutoff= 1)
ego_BP <- enrichGO(gene=DEG, OrgDb= 'org.Hs.eg.db', keyType='SYMBOL', ont="BP", pvalueCutoff= 1,qvalueCutoff= 1)

p_BP <- barplot(ego_BP,showCategory = 10) + ggtitle("Biological process")
p_CC <- barplot(ego_CC,showCategory = 10) + ggtitle("Cellular component")
p_MF <- barplot(ego_MF,showCategory = 10) + ggtitle("Molecular function")
plotc <- p_BP/p_CC/p_MF
plotc
ggsave('result/6.enrichGO.png', plotc, width = 10,height = 16)

ego_BP <- data.frame(ego_BP)
ego_CC <- data.frame(ego_CC)
ego_MF <- data.frame(ego_MF)
write.csv(ego_BP,'result/6.enrichGO_BP.csv')
write.csv(ego_CC,'result/6.enrichGO_CC.csv')
write.csv(ego_MF,'result/6.enrichGO_MF.csv')



## === KEGG
genelist <- bitr(gene=DEG, fromType="SYMBOL", toType="ENTREZID", OrgDb='org.Hs.eg.db')
genelist <- pull(genelist,ENTREZID)               
ekegg <- enrichKEGG(gene = genelist, organism = 'hsa', pvalueCutoff = 1, qvalueCutoff = 1)
p1 <- barplot(ekegg, showCategory=10)
p2 <- dotplot(ekegg, showCategory=10)
plotc = p1/p2
plotc
ggsave('result/6.enrichKEGG.png', plot = plotc, width = 8, height = 10)

ekegg <- data.frame(ekegg)
write.csv(ekegg,'result/6.enrichKEGG.csv')



## === 其他数据库通路
geneset <- read.gmt("data/MsigDB/v7.4/h.all.v7.4.symbols.gmt")
table(geneset$term)
geneset$term <- gsub(pattern = "HALLMARK_","", geneset$term)
geneset$term <- str_to_title(geneset$term)

my_path <- enricher(gene=DEG, pvalueCutoff = 1, qvalueCutoff = 1, TERM2GENE=geneset)
p1 <- barplot(my_path, showCategory=15,color = "pvalue")
p1
ggsave("result/6.enrich_HALLMARK.png", plot = p1, width = 8, height = 7)

my_path <- data.frame(my_path)
write.csv(my_path,"result/6.enrich_HALLMARK.csv")

2、GSEA

# 清空当前环境变量
rm(list = ls())
options(stringsAsFactors = F)

# 加载包
library(GSEABase)
library(clusterProfiler)
library(enrichplot)
library(ggplot2)
library(stats)

# 加载数据
load("data/Step03-edgeR_nrDEG.Rdata")
DEG <- DEG_edgeR_symbol

## 构造GSEA分析数据
# 去掉没有配对上symbol的行
DEG <- DEG[!is.na(DEG$SYMBOL),]

# 去掉重复行
DEG <- DEG[!duplicated(DEG$SYMBOL),]

geneList <- DEG$logFC
names(geneList) <- DEG$SYMBOL
head(geneList)
geneList <- sort(geneList,decreasing = T)
head(geneList)
tail(geneList)


# 选择gmt文件
geneset <- read.gmt("data/MsigDB/v7.4/c5.go.bp.v7.4.symbols.gmt")

# 运行,输出全部结果
egmt <- GSEA(geneList, TERM2GENE=geneset, pvalueCutoff = 1)

#出点图
dotplot(egmt)

#按p值出点图
dotplot(egmt,color="pvalue")   

# 单个通路图
# 按照通路名
gseaplot2(egmt, "GOBP_MITOCHONDRIAL_GENOME_MAINTENANCE",  
          title = "GOBP_MITOCHONDRIAL_GENOME_MAINTENANCE")
# 按照行数
gseaplot2(egmt, 10, color="red", pvalue_table = T)

#按第一到第十个出图，不显示p值
gseaplot2(egmt, 1:10, color="red") 


# 保存结果
go_gsea <- as.data.frame(egmt@result)
write.csv(go_gsea,"result/6.gsea_go_fc.csv",row.names = F)

代码均来自于生信技能树张娟老师