懒癌的我全靠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
基因表达聚类，黄线和绿线表示随时间变化表达量相差小的基因，红线和紫线表明随时间变化表大量相差大的基因。
