懶癌的我全靠Jimmy催著學(xué)習(xí)鉴吹，4月份答應(yīng)交這份作業(yè)的东揣，離中間被催又過了一個(gè)多月，今天完成啦铲敛。
數(shù)據(jù)來源： 擬南芥RNA-seq轉(zhuǎn)錄組分析-salmon（Mac版）
上游分析為salmon直接定量分析（純Mac版）待错，下游分析如下：
1. 軟件安裝
如果有直接library("包名")籽孙，但這一步?jīng)]有安裝的，直接用install.packages("包名安裝")火俄，如果安裝不成功犯建，用BiocManager::install("包名"),詳情參考R包升級(jí)報(bào)錯(cuò)成常態(tài),搜索告訴你。瓜客。适瓦。
#清空環(huán)境
rm(list = ls())
options(stringsAsFactors = F)
#設(shè)置鏡像
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的工作目錄設(shè)置到存儲(chǔ)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 注釋信息獲取
#注釋信息的數(shù)據(jù)庫,目前最新的工具包叫做AnnotationHub
library(AnnotationHub)
#此報(bào)錯(cuò)不定期出現(xiàn)，有時(shí)加載很多次都會(huì)出現(xiàn)如圖ah1的報(bào)錯(cuò)玻熙，可能是網(wǎng)絡(luò)原因否彩，ah2位加載成功，
ah <- AnnotationHub()
#查找包中擬南芥的數(shù)據(jù)
ath <- query(ah,'Arabidopsis thaliana')
#下載最新的注釋ID,此步驟會(huì)出現(xiàn)ah一步類似的報(bào)錯(cuò)揭芍，多試幾次
ath_tx <- ath[['AH52247']]
#提取數(shù)據(jù)庫中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 加載數(shù)據(jù)
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為列，將列都轉(zhuǎn)化為整數(shù)
exprSet=apply(tmp,2, as.integer)
rownames(exprSet)=rownames(tmp)
#查看表達(dá)兩信息
head(exprSet)
dim(exprSet)
exprSet=as.data.frame(exprSet)
save(exprSet,file=paste0('quants-exprSet.Rdata'))
#準(zhǔn)備樣本信息
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和樣品名對(duì)應(yīng)
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. 差異分析
分三步： 構(gòu)建矩陣-均一化-表達(dá)矩陣
3.1 構(gòu)建矩陣
library('DESeq2')
dds<-DESeqDataSetFromTximport(txi,sampleTable,design = ~group_list)
dim(dds)
#數(shù)據(jù)過濾
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 表達(dá)矩陣
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. 時(shí)序表達(dá)聚類分析
用3.2 均一化(normalisation)后的表達(dá)矩陣
library("Mfuzz")
library('dplyr')
count <- data.matrix(exprSet_new)
eset <- new("ExpressionSet",exprs = count)
# 根據(jù)標(biāo)準(zhǔn)差去除樣本間差異太小的基因
eset <- filter.std(eset,min.std=0)
eset <- standardise(eset)
image.png
# 如何決定聚類個(gè)數(shù)筷转？
c <- 6
# 評(píng)估出最佳的m值
m <- mestimate(eset)
# 聚類
cl <- mfuzz(eset, c = c, m = m)
mfuzz.plot(
eset,
cl,
mfrow=c(2,3),
new.window= FALSE)
聚類個(gè)數(shù)位6
image.png
聚類個(gè)數(shù)位16
image.png
基因表達(dá)聚類姑原，黃線和綠線表示隨時(shí)間變化表達(dá)量相差小的基因，紅線和紫線表明隨時(shí)間變化表大量相差大的基因呜舒。
