Mfuzz进行时序表达聚类分析（Mac）

懒癌的我全靠Jimmy催着学习，4月份答应交这份作业的，离中间被催又过了一个多月，今天完成啦。

数据来源： 拟南芥RNA-seq转录组分析-salmon（Mac版）
上游分析为salmon直接定量分析（纯Mac版），下游分析如下：

1. 软件安装

如果有直接library("包名")，但这一步没有安装的，直接用install.packages("包名安装")，如果安装不成功，用BiocManager::install("包名"),详情参考R包升级报错成常态,搜索告诉你。。。

#清空环境
rm(list = ls()) 
options(stringsAsFactors = F)
#设置镜像
options()$repos 
options()$BioC_mirror
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
options()$repos 
options()$BioC_mirror

install.packages("tidyverse") ; library(tidyverse)
install.packages("optparse") ; library(optparse)
install.packages("UpSetR") ; library(UpSetR)
install.packages("rjson") ; library(rjson)
install.packages("Mfuzz"); library(Mfuzz)
# https://bioconductor.org/packages/release/bioc/html/GEOquery.html
if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
BiocManager::install(c("tximport","KEGG.db","DESeq2","edgeR" ,"org.At.tair.db","pheatmap","AnnotationHub"),ask = F,update = F)
BiocManager::install(c("clusterProfiler","limma","DOSE","GenomicFeatures","RColorBrewer"),ask = F,update = F)

2. 打开Rstudio，加载salmon得到的counts

2.1 加载到counts所在的路径

#setwd将Rstudio的工作目录设置到存储counts的路径
setwd("~/data/")
getwd()
dir=getwd()
files=list.files(pattern = "*sf",dir,recursive = T)
files=file.path(dir,files)
all(file.exists(files))

2.2 注释信息获取

#注释信息的数据库,目前最新的工具包叫做AnnotationHub
library(AnnotationHub)
#此报错不定期出现，有时加载很多次都会出现如图ah1的报错，可能是网络原因，ah2位加载成功， 
ah <- AnnotationHub()
#查找包中拟南芥的数据
ath <- query(ah,'Arabidopsis thaliana')
#下载最新的注释ID,此步骤会出现ah一步类似的报错，多试几次
ath_tx <- ath[['AH52247']]
#提取数据库中GENEID信息
columns(ath_tx)
k <- keys(ath_tx,keytype = "GENEID")
df <- select(ath_tx, keys=k, keytype = "GENEID",columns = "TXNAME")
tx2gene <- df[,2:1]

2.3 加载数据

library('tximport')
library('readr')
txi=tximport(files ,type = "salmon",tx2gene = tx2gene)
head(txi)
head(txi$length)
files
# 加载stringr包,为了得到将txi的列名定义为样本名
library(stringr)
# 以'\'为分隔符，将含有quant.sf的路径分割，并提第七列---取样本名（ERR1698194_quant）
t1=sapply(strsplit(files,'\\/'),function(x)x[6])
# txi 的列counts的列名为样本名,
colnames(txi$counts)=sapply(strsplit(t1,'_'),function(x)x[1])
tmp=txi$counts
head(tmp)
# 1为行，2为列，将列都转化为整数
exprSet=apply(tmp,2, as.integer)
rownames(exprSet)=rownames(tmp)
#查看表达两信息
head(exprSet)
dim(exprSet)
exprSet=as.data.frame(exprSet)
save(exprSet,file=paste0('quants-exprSet.Rdata'))

#准备样本信息
sampleTable=read.csv("sampleTable.txt",sep="\t",header = FALSE)
colnames(sampleTable)=c("sample","group_list")
rownames(sampleTable)=sampleTable[,1]
sampleTable=sampleTable[,-1,drop=FALSE]
sampleTable$group_list=paste0("day",sampleTable$group_list)

# 将counts和样品名对应
names(txi)
head(txi$counts)
head(txi$length)
colnames(txi$length)=colnames(txi$counts)
colnames(txi$abundance)=colnames(txi$counts)
txi$counts[1:4,1:4]

3. 差异分析

分三步： 构建矩阵-均一化-表达矩阵

3.1 构建矩阵

library('DESeq2')
dds<-DESeqDataSetFromTximport(txi,sampleTable,design = ~group_list)
dim(dds)
#数据过滤
dds <- dds[rowSums(counts(dds))>1,]
dim(dds)

suppressMessages(dds2 <- DESeq(dds))
res <- results(dds2)
summary(res)

3.2 均一化

rld=rlogTransformation(dds2)
exprSet_new=assay(rld)
head(exprSet_new)
dim(exprSet_new)

3.3 表达矩阵

resultsNames(dds2)
res0vs1= results(dds2,contrast = c("group_list","day1","day0"))
resOrdered1=res0vs1[order(res0vs1$padj),]
resOrdered1=as.data.frame(resOrdered1)
head(resOrdered1)

library("org.At.tair.db")
library("KEGG.db")
library("clusterProfiler")
library('ggplot2')
resOrdered1$gene_id=rownames(resOrdered1)
id2symbol=toTable(org.At.tairSYMBOL)
resOrdered1=merge(resOrdered1,id2symbol,by='gene_id')

DEG=resOrdered1
colnames(DEG)=colnames(DEG)=c('gene_id' ,'baseMean','logFC','lfcSE','stat','pvalue' , 'P.Value' , 'symbol')

DEG$symbol=as.character(DEG$symbol)
DEG_filter=DEG[nchar(DEG$symbol)>1,]
DEG_filter=DEG_filter[!is.na(DEG_filter$symbol),]

4. 时序表达聚类分析

用3.2 均一化(normalisation)后的表达矩阵

library("Mfuzz")
library('dplyr')

count <- data.matrix(exprSet_new)
eset <- new("ExpressionSet",exprs = count)
# 根据标准差去除样本间差异太小的基因
eset <- filter.std(eset,min.std=0)
eset <- standardise(eset)

image.png

# 如何决定聚类个数？
c <- 6
#  评估出最佳的m值
m <- mestimate(eset)
# 聚类
cl <- mfuzz(eset, c = c, m = m)

mfuzz.plot(
  eset,
  cl,
  mfrow=c(2,3),
  new.window= FALSE)

聚类个数位6

image.png

聚类个数位16

image.png

基因表达聚类，黄线和绿线表示随时间变化表达量相差小的基因，红线和紫线表明随时间变化表大量相差大的基因。