转录组差异表达分析和火山图可视化

利用R包DEseq2进行差异表达分析和可视化

- count数矩阵
- 差异分析
- - 1. 安装并载入R包
  - 2. count数矩阵导入并对矩阵进行数据处理
  - 3. 查看样本相关性并采用热图展示
  - 4. hclust对样本进行聚类分析
  - 5. 构建原始dds矩阵并保存为Rdata对象
  - 6. 原始dds矩阵标准化并保存
  - 7. 提取差异分析的结果
  - 8. 绘制火山图
  - 9. 简单gene ID转换
- 参考文件

首先附上文献中的坚定差异基因的流程图。

count数矩阵

在Linux下，通过HISAT2 对fastq数据文件进行比对，FeatureCounts软件进行基因水平定量，得到count数矩阵。之后便可以载入R语言中进行差异分析。

差异分析

第一次分析RNA-seq数据，走到这一步相对容易了许多。转录组数据分析主要参考了生信技能树Jimmy老师的相关课程及推文。
RNA-seq的read count普遍认为符合泊松分布，但是之前分析过的芯片数据符合正态分布，所以筛选DEGs的方法有一定差别。

1. 安装并载入R包

# 设置R语言镜像# options(BioC_mirror="http://mirrors.tuna.tsinghua.edu.cn/bioconductor/")# options("repos" = c(CRAN="http://mirrors.tuna.tsinghua.edu.cn/CRAN/"))# 安装R包# if(!require(c("ggthemes","ggpubr","ggthemes","ggrepel"))) install.packages(c("ggthemes","ggpubr","ggthemes","ggrepel"))# BiocManager::install("DESeq2")

#载入R包suppressPackageStartupMessages(library(DESeq2))suppressPackageStartupMessages(library(ggpubr))suppressPackageStartupMessages(library(ggplot2))suppressPackageStartupMessages(library(ggrepel))suppressPackageStartupMessages(library(ggthemes))

2. count数矩阵导入并对矩阵进行数据处理

exprset <- read.table("RNA-seq_counts_matrix.csv",sep = ",",header = T,check.names = F)rownames(exprset) <- exprset[,1]exprset <- exprset[,-1]exprset <- exprset[apply(exprset,1,function(x) sum(x > 1) > 1),] #先判断值是否为0，得到逻辑值，再取和，判断0的个数是否小于3dim(exprset)# 7428    4head(exprset)

head(exprset)

	control1	control2	treat1	treat2
ENSMUSG00000000028	27	0	0	6
ENSMUSG00000000088	124	268	87	313
ENSMUSG00000000094	5	12	2	0
ENSMUSG00000000131	17	5	6	5
ENSMUSG00000000134	23	79	0	1
ENSMUSG00000000142	6	10	0	0

3. 查看样本相关性并采用热图展示

expcor <- cor(exprset, method = "spearman")head(expcor)pheatmap::pheatmap(expcor, clustering_method = "average",                   treeheight_row = 0,treeheight_col = 0,                   display_numbers = T)

expcor data

	control1	control2	treat1	treat2
control1	1.0000000	0.7089970	0.2366665	0.0209855
control2	0.7089970	1.0000000	0.2990182	0.0866515
treat1	0.2366665	0.2990182	1.0000000	0.4533486
treat2	0.0209855	0.0866515	0.4533486	1.0000000

热图展示

4. hclust对样本进行聚类分析

# t_exprset <- t(exprset)# t_exprset <- t_exprset[,names(sort(apply(t_exprset,2,mad),decreasing = T))[1:500]]# out.dist <- dist(t_exprset,method = 'euclidean')# out.hclust <- hclust(out.dist,method = 'complete')# rect.hclust(out.hclust,k=3)# plot(out.hclust,xlab = "",main = "")

5. 构建原始dds矩阵并保存为Rdata对象

group_list <- factor(c(rep("untrt",2),rep("trt",2))) #因子型变量group_listtable(group_list)## group_list##   trt untrt ##     2     2colData <- data.frame(row.names = colnames(exprset),                        group_list = group_list)colDatadds <- DESeqDataSetFromMatrix(countData = exprset,                               colData = colData,                               design = ~group_list) #~在R里面用于构建公式对象，~左边为因变量，右边为自变量。head(dds)## class: DESeqDataSet ## dim: 6 4 ## metadata(1): version## assays(1): counts## rownames(6): ENSMUSG00000000028 ENSMUSG00000000088 ...##   ENSMUSG00000000134 ENSMUSG00000000142## rowData names(0):## colnames(4): control1 control2 treat1 treat2## colData names(1): group_listtem_f <- 'RNA-seq_DESeq2-dds.Rdata'

colData

	group_list
control1	untrt
control2	untrt
treat1	trt
treat2	trt

6. 原始dds矩阵标准化并保存

if (!file.exists(tem_f)) {    dds <- DESeq(dds) # 标准化    save(dds,file = tem_f)  }load(file = tem_f)# 结果用`result()`函数提取res <- results(dds,              contrast = c("group_list","untrt","trt")) # 差异分析结果resOrdered <- res[order(res$padj),] # 对结果按照调整后的p值进行排序head(resOrdered)summary(res)## ## out of 7428 with nonzero total read count## adjusted p-value < 0.1## LFC > 0 (up)       : 465, 6.3%## LFC < 0 (down)     : 507, 6.8%## outliers [1]       : 0, 0%## low counts [2]     : 2160, 29%## (mean count < 4)## [1] see 'cooksCutoff' argument of ?results## [2] see 'independentFiltering' argument of ?results

head(resOrdered)

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj
ENSMUSG00000061787	1308.2358	-9.456575	1.564545	-6.044298	0e+00	3.60e-06
ENSMUSG00000064370	1304.1697	-13.689071	2.284209	-5.992916	0e+00	3.60e-06
ENSMUSG00000096745	667.1955	-12.722138	2.066186	-6.157306	0e+00	3.60e-06
ENSMUSG00000096363	320.2598	-11.663243	2.067930	-5.640056	0e+00	2.24e-05
ENSMUSG00000031504	229.8465	-11.184637	2.077845	-5.382805	1e-07	4.24e-05
ENSMUSG00000038900	583.4616	-8.543657	1.597311	-5.348775	1e-07	4.24e-05

7. 提取差异分析的结果

DEG <- as.data.frame(resOrdered)DESeq2_DEG <- na.omit(DEG)diff <- subset(DESeq2_DEG,pvalue < 0.05) #先筛选P值up <- subset(diff,log2FoldChange > 2) #上调down <- subset(diff,log2FoldChange < -2) #下调#可利用`write.csv()`函数保存文件

8. 绘制火山图

DEG_data <- DESeq2_DEGDEG_data$logP <- -log10(DEG_data$padj) # 对差异基因矫正后p-value进行log10()转换dim(DEG_data)## [1] 5268    7#将基因分为三类：not-siginficant，up，dowm#将adj.P.value小于0.05，logFC大于2的基因设置为显著上调基因#将adj.P.value小于0.05，logFC小于-2的基因设置为显著上调基因DEG_data$Group <- "not-siginficant"DEG_data$Group[which((DEG_data$padj < 0.05) & DEG_data$log2FoldChange > 2)] = "up-regulated"DEG_data$Group[which((DEG_data$padj < 0.05) & DEG_data$log2FoldChange < -2)] = "down-regulated"table(DEG_data$Group)## ##  down-regulated not-siginficant    up-regulated ##             336            4659             273DEG_data <- DEG_data[order(DEG_data$padj),]#对差异表达基因调整后的p值进行排序#火山图中添加点(数据构建)up_label <- head(DEG_data[DEG_data$Group == "up-regulated",],1)down_label <- head(DEG_data[DEG_data$Group == "down-regulated",],1)deg_label_gene <- data.frame(gene = c(rownames(up_label),rownames(down_label)),                                label = c(rownames(up_label),rownames(down_label)))DEG_data$gene <- rownames(DEG_data)DEG_data <- merge(DEG_data,deg_label_gene,by = 'gene',all = T)#不添加labelggscatter(DEG_data,x = "log2FoldChange",y = "logP",          color = "Group",          palette = c("green","gray","red"),          repel = T,          ylab = "-log10(Padj)",          size = 1) +   theme_base()+  scale_y_continuous(limits = c(0,8))+  scale_x_continuous(limits = c(-18,18))+  geom_hline(yintercept = 1.3,linetype = "dashed")+  geom_vline(xintercept = c(-2,2),linetype = "dashed")

#添加特定基因labelggscatter(DEG_data,x = "log2FoldChange",y = "logP",          color = "Group",          palette = c("green","gray","red"),          label = DEG_data$label,          repel = T,          ylab = "-log10(Padj)",          size = 1) +   theme_base()+  theme(element_line(size = 0),element_rect(size = 1.5))+ #坐标轴线条大小设置  scale_y_continuous(limits = c(0,8))+  scale_x_continuous(limits = c(-18,18))+  geom_hline(yintercept = 1.3,linetype = "dashed")+  geom_vline(xintercept = c(-2,2),linetype = "dashed")

9. 简单gene ID转换

对于初学者来说如果要对gene ID进行转换，可利用Ensembl数据库的BioMart工具。因为相对于R包biomaRt，界面化的操作更加易懂，快捷。BioMart网页工具的原始界面如下所示：

  其中左侧菜单栏分别是Dataset--选择相关物种参考基因组；  Filters--选择数据gene ID的类型，并输入gene ID，也存在其他类型的ID输入；  Attributes--选择需要输出的ID类型；  点击Result可以输出结果，并且支持文件下载。

第一次写推文，请大家多提宝贵意见！
##如有侵权请联系作者删除！

参考文件

[1] https://mp.weixin.qq.com/s/uDnFJC0szOHtO2NqREz2wA

[2] https://www.jianshu.com/p/3a0e1e3e41d0

[3] https://www.bioconductor.org/help/workflows/RNAseq123/

[4] https://www.bioconductor.org/help/workflows/rnaseqGene/

[5] http://www.biotrainee.com/forum.phpmod=viewthread&tid=1750#lastpost

[6] https://mp.weixin.qq.com/s/ZYB06Yudck2hD0qWJKJcwQ

本站仅提供存储服务，所有内容均由用户发布，如发现有害或侵权内容，请点击举报。