1.## 加載R包
## 下載數(shù)據(jù)掸绞，如果文件夾中有會(huì)直接讀入
gset = getGEO('GSE32575', destdir=".",getGPL = F)
# 獲取ExpressionSet對(duì)象婆排，包括的表達(dá)矩陣和分組信息
gset=gset[[1]]
2.#通過(guò)pData函數(shù)獲取分組信息
## 獲取分組信息,點(diǎn)開(kāi)查閱信息
pdata=pData(gset)
## 只要后36個(gè),本次選擇的這36個(gè)是配對(duì)的赡茸。
## 所以，別人的芯片我們也不是全要丑勤，我們只要適合自己的數(shù)據(jù)
group_list=c(rep('before',18),rep('after',18))
group_list=factor(group_list)
## 強(qiáng)制限定順序
group_list <- relevel(group_list, ref="before")
3.#通過(guò)exprs函數(shù)獲取表達(dá)矩陣并校正
exprSet=exprs(gset)
可以先簡(jiǎn)單看一下整體樣本的表達(dá)情況，
其中exprSet[,-c(1:12)]的意思是去掉這個(gè)數(shù)據(jù)的前12個(gè)慎恒，因?yàn)槲覀冃枰暮竺?6個(gè).
boxplot(exprSet[,-c(1:12)],outline=FALSE, notch=T,col=group_list, las=2)
需要人工校正一下矛紫，用的方法類似于Quntile Normalization
這里我們用limma包內(nèi)置的一個(gè)函數(shù)
library(limma)
exprSet=normalizeBetweenArrays(exprSet)
boxplot(exprSet,outline=FALSE, notch=T,col=group_list, las=2)
4.#判斷是否需要進(jìn)行數(shù)據(jù)轉(zhuǎn)換
我們通過(guò)分組信息已經(jīng)知道赎瞎，前12個(gè)樣本本次不需要，所以先去掉
exprSet = as.data.frame(exprSet)[,-seq(1,12)]
肉眼看到數(shù)字很大颊咬，這時(shí)候需要log轉(zhuǎn)換(選log2)务甥。
此外還有一個(gè)腳本可以自動(dòng)判斷是否需要log轉(zhuǎn)換，這個(gè)腳本來(lái)自于GEO2R贪染，我改寫了一點(diǎn)點(diǎn)缓呛。
ex <- exprSet
qx <- as.numeric(quantile(ex, c(0., 0.25, 0.5, 0.75, 0.99, 1.0), na.rm=T))
LogC <- (qx[5] > 100) ||
? (qx[6]-qx[1] > 50 && qx[2] > 0) ||
? (qx[2] > 0 && qx[2] < 1 && qx[4] > 1 && qx[4] < 2)
if (LogC) { ex[which(ex <= 0)] <- NaN
exprSet <- log2(ex)
print("log2 transform finished")}else{print("log2 transform not needed")}
這個(gè)語(yǔ)句會(huì)自動(dòng)判斷催享，如果已經(jīng)log話就跳過(guò)杭隙，并提示
"log2 transform not needed"
如果沒(méi)有l(wèi)og話，他自動(dòng)log2因妙，并且提示：
"log2 transform finished"
但是我們發(fā)現(xiàn)一個(gè)問(wèn)題痰憎，這個(gè)表達(dá)數(shù)據(jù)行名是探針，
不是我們熟悉的基因名稱如果TP53攀涵，BRCA1這樣的铣耘，所以我們需要注釋他。
5.#獲取注釋信息
通常情況下我們使用R包來(lái)注釋以故，R包和平臺(tái)的注釋信息對(duì)應(yīng)關(guān)系我整理了一個(gè)表格 platformMap蜗细，
是一個(gè)文本文件，在文末有獲取方法。
platformMap <- data.table::fread("platformMap.txt")
我們根據(jù)上述帖子里面提到的platformMap炉媒，找到對(duì)應(yīng)的R包是：
"illuminaHumanv2.db"
index = gset@annotation
index = gset@annotation
platformDB = paste0(platformMap$bioc_package[grep(index,platformMap$gpl)],".db")
安裝R包并加載
# 第一句是為了增加鏡像
if(length(getOption("BioC_mirror"))==0) options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
## 沒(méi)有這個(gè)包就給你裝踪区，有就不裝
if(!require("illuminaHumanv2.db")) BiocManager::install("illuminaHumanv2.db",update = F,ask = F)
## 加載R包
library(illuminaHumanv2.db)
獲取探針對(duì)應(yīng)的symbol信息
## 獲取表達(dá)矩陣的行名，就是探針名稱
probeset <- rownames(exprSet)
## 使用lookup函數(shù)吊骤，找到探針在illuminaHumanv2.db中的對(duì)應(yīng)基因名稱
## 如果分析別的芯片數(shù)據(jù)缎岗，把illuminaHumanv2.db更換即可
SYMBOL <-? annotate::lookUp(probeset,"illuminaHumanv2.db", "SYMBOL")
## 轉(zhuǎn)換為向量
symbol = as.vector(unlist(SYMBOL))
制作probe2symbol轉(zhuǎn)換文件
probe2symbol = data.frame(probeset,symbol,stringsAsFactors = F)
6.探針轉(zhuǎn)換與基因去重
一個(gè)基因會(huì)對(duì)應(yīng)對(duì)個(gè)探針，有些基因名稱會(huì)是重復(fù)的白粉，這些都需要處理传泊。
對(duì)于，多個(gè)探針鸭巴，我們選取在樣本中平均表達(dá)量最高的探針作為對(duì)應(yīng)基因的表達(dá)量眷细。
一下代碼完成所有事情，而且可以復(fù)用鹃祖。
library(dplyr)
library(tibble)
exprSet <- exprSet %>%
? rownames_to_column(var="probeset") %>%
? #合并探針的信息
? inner_join(probe2symbol_df,by="probeset") %>%
? #去掉多余信息
? select(-probeset) %>%
? #重新排列
? select(symbol,everything()) %>%
? #求出平均數(shù)(這邊的點(diǎn)號(hào)代表上一步產(chǎn)出的數(shù)據(jù))
? mutate(rowMean =rowMeans(.[grep("GSM", names(.))])) %>%
? #去除symbol中的NA
? filter(symbol != "NA") %>%
? #把表達(dá)量的平均值按從大到小排序
? arrange(desc(rowMean)) %>%
? # symbol留下第一個(gè)
? distinct(symbol,.keep_all = T) %>%
? #反向選擇去除rowMean這一列
? select(-rowMean) %>%
? # 列名變成行名
? column_to_rownames(var = "symbol")
7.差異分析
差異分析的矩陣構(gòu)建有兩種方式
我們選取格式比較簡(jiǎn)單的薪鹦。如果沒(méi)有配對(duì)信息，差異分析這樣做：
design=model.matrix(~ group_list)
fit=lmFit(exprSet,design)
fit=eBayes(fit)
allDiff=topTable(fit,adjust='fdr',coef="group_listafter",number=Inf,p.value=0.05)
其中allDiff得到6799行惯豆。
因?yàn)槲覀冞@里實(shí)際上是有配對(duì)信息的池磁，差異分析應(yīng)該這樣做：
pairinfo = factor(rep(1:18,2))
design=model.matrix(~ pairinfo+group_list)
fit=lmFit(exprSet,design)
fit=eBayes(fit)
allDiff_pair=topTable(fit,adjust='BH',coef="group_listafter",number=Inf,p.value=0.05)
得到的差異分析結(jié)果allDiff_pair有7650個(gè)，比直接做差異分析要多接近1000個(gè)楷兽。
8.作圖驗(yàn)證(非必要)
差異分析的結(jié)果地熄，配對(duì)和非配對(duì)理解起來(lái)比較抽象，我們用圖片直觀地感受一下芯杀。
首先我們需要得到ggplot2喜歡的數(shù)據(jù)格式端考，行是觀測(cè)，列是變量揭厚，這也叫clean data却特，tide data, 清潔數(shù)據(jù)。
首先基因名稱需要在列筛圆，所以需要用t函數(shù)裂明，行列轉(zhuǎn)置。
data_plot = as.data.frame(t(exprSet))
data_plot = data.frame(pairinfo=rep(1:18,2),
? ? ? ? ? ? ? ? ? ? ? group=group_list,
? ? ? ? ? ? ? ? ? ? ? data_plot,stringsAsFactors = F)
作圖展示:
? 挑選了一個(gè)基因CAMKK2
library(ggplot2)
ggplot(data_plot, aes(group,CAMKK2,fill=group)) +
? geom_jitter(aes(colour=group), size=2, alpha=0.5)+
? xlab("") +
? ylab(paste("Expression of ","CAMKK2"))+
? theme_classic()+
? theme(legend.position = "none")
看起來(lái)雜亂一片太援，根本得不到有用信息闽晦，我們加入他的配對(duì)信息，再來(lái)作圖：
library(ggplot2)
ggplot(data_plot, aes(group,CAMKK2,fill=group)) +
? geom_boxplot() +
? geom_point(size=2, alpha=0.5) +
? geom_line(aes(group=pairinfo), colour="black", linetype="11") +
? xlab("") +
? ylab(paste("Expression of ","CAMKK2"))+
? theme_classic()+
? theme(legend.position = "none")
再來(lái)看看真正有差異的基因吧提岔，比如：
library(ggplot2)
ggplot(data_plot, aes(group,CH25H,fill=group)) +
? geom_boxplot() +
? geom_point(size=2, alpha=0.5) +
? geom_line(aes(group=pairinfo), colour="black", linetype="11") +
? xlab("") +
? ylab(paste("Expression of ","CH25H"))+
? theme_classic()+
? theme(legend.position = "none")
我們還可以批量地作圖仙蛉，這段不必要，可以自己一張張畫碱蒙，我只是展示一下如何批量畫出差異最大的8個(gè)基因：
library(dplyr)
library(tibble)
allDiff_arrange <- allDiff_pair %>%
? rownames_to_column(var="genesymbol") %>%
? arrange(desc(abs(logFC)))
genes <- allDiff_arrange$genesymbol[1:8]
plotlist <- lapply(genes, function(x){
? data =data.frame(data_plot[,c("pairinfo","group")],gene=data_plot[,x])
? ggplot(data, aes(group,gene,fill=group)) +
? ? geom_boxplot() +
? ? geom_point(size=2, alpha=0.5) +
? ? geom_line(aes(group=pairinfo), colour="black", linetype="11") +
? ? xlab("") +
? ? ylab(paste("Expression of ",x))+
? ? theme_classic()+
? ? theme(legend.position = "none")
})
library(cowplot)
plot_grid(plotlist=plotlist, ncol=4,labels = LETTERS[1:8])
#9.后續(xù)分析
差異分析后還有很多操作荠瘪，
畫一張熱圖
library(pheatmap)
## 設(shè)定差異基因閾值，減少差異基因用于提取表達(dá)矩陣
allDiff_pair=topTable(fit,adjust='BH',coef="group_listafter",number=Inf,p.value=0.05,lfc =0.5)
##提前部分?jǐn)?shù)據(jù)用作熱圖繪制
heatdata <- exprSet[rownames(allDiff_pair),]
##制作一個(gè)分組信息用于注釋
annotation_col <- data.frame(group_list)
rownames(annotation_col) <- colnames(heatdata)
#如果注釋出界，可以通過(guò)調(diào)整格子比例和字體修正
pheatmap(heatdata, #熱圖的數(shù)據(jù)
? ? ? ? cluster_rows = TRUE,#行聚類
? ? ? ? cluster_cols = TRUE,#列聚類哀墓，可以看出樣本之間的區(qū)分度
? ? ? ? annotation_col =annotation_col, #標(biāo)注樣本分類
? ? ? ? annotation_legend=TRUE, # 顯示注釋
? ? ? ? show_rownames = F,# 顯示行名
? ? ? ? show_colnames = F,# 顯示列名
? ? ? ? scale = "row", #以行來(lái)標(biāo)準(zhǔn)化鞭莽，這個(gè)功能很不錯(cuò)
? ? ? ? color =colorRampPalette(c("blue", "white","red"))(100))
我們還可以畫一個(gè)火山圖
library(ggplot2)
library(ggrepel)
library(dplyr)
libr
data <- topTable(fit,adjust='BH',coef="group_listafter",number=Inf)
data$significant <- as.factor(data$adj.P.Val<0.05 & abs(data$logFC) > 0.5)
data$gene <- rownames(data)
ggplot(data=data, aes(x=logFC, y =-log10(adj.P.Val),color=significant)) +
? geom_point(alpha=0.8, size=1.2,col="black")+
? geom_point(data=subset(data, logFC > 0.5),alpha=0.8, size=1.2,col="red")+
? geom_point(data=subset(data, logFC < -0.5),alpha=0.6, size=1.2,col="blue")+
? labs(x="log2 (fold change)",y="-log10 (adj.P.Val)")+
? theme(plot.title = element_text(hjust = 0.4))+
? geom_hline(yintercept = -log10(0.05),lty=4,lwd=0.6,alpha=0.8)+
? geom_vline(xintercept = c(0.5,-0.5),lty=4,lwd=0.6,alpha=0.8)+
? theme_bw()+
? theme(panel.border = element_blank(),
? ? ? ? panel.grid.major = element_blank(),
? ? ? ? panel.grid.minor = element_blank(),?
? ? ? ? axis.line = element_line(colour = "black")) +
? geom_point(data=subset(data, abs(logFC) > 1),alpha=0.8, size=3,col="green")+
? geom_text_repel(data=subset(data, abs(logFC) > 1),
? ? ? ? ? ? ? ? ? aes(label=gene),col="black",alpha = 0.8)
d = data.frame(obs = c(1, 2, 3), treat = c("A", "B", "A" ),
? ? ? ? ? ? ? weight = c(2.3, NA, 9))
d
# 保存為簡(jiǎn)單的文本文件
write.table(d,file="SDDDD.txt")
write.table(d,file="SDDDD.txt",row.names=F,quote = F)
# row.names = F指定行名不寫入文件，quote = F表示變量名不放入雙引號(hào)中
# 保存為逗號(hào)分割的文本文件
write.csv(d,file = "SDD.csv")
# 保存為R格式文件(直接用save）
save(d,file="SDDDD.RData")
# 在經(jīng)過(guò)一段時(shí)間的分析后麸祷，常需要將工作空間的映像保存澎怒，命令為：
save.image()
# 它等價(jià)于 save(list = ls(all=TRUE), file = ".RData")
數(shù)據(jù)的讀取
可以使用read.table()、scan()阶牍、readfwf()函數(shù)讀入存儲(chǔ)在文本文件中的數(shù)據(jù)
A=read.table("SDDDD.txt",header=TRUE)? # header = TRUE 表明讀入的第一行是變量名
A
# read.table有四個(gè)變形
# read.csv(), read.csv2(), read.delim(), read.delim2()
# 使用函數(shù)scan()喷面，它比read.table()更加靈活，唯一的區(qū)別是scan()可以指定變量的類型
mydata=scan("SDD.csv",what=list(obs=0,treat="",weight=0),sep=",",skip = 1)
mydata
另一個(gè)重要區(qū)別是scan()函數(shù)可以創(chuàng)建不同的對(duì)象走孽，在缺省情況下（what被省略）惧辈，將創(chuàng)建一個(gè)數(shù)值型向量。
# 使用函數(shù)read.fwf()磕瓷，它可以用來(lái)讀取文件中一些固定寬度格式的數(shù)據(jù)盒齿。
mydata = read.fwf("SDDDD.txt", widths=c(1,4,3), col.names=c("X", "Y", "Z"))
widths=c(1,4,3)
install.packages("‘BiocManager",repos="s https://github.com/Bioconductor/BiocManager/issues")
library(RODBC)
z = odbcConnectExcel("LYC.xls")
attach(mtcars)
attach(mtcars)
mtcars2 = data.frame(mtcars[, c(1,4)])
head(mtcars2)
load("SDDDD.RData")
demo("graphics")
demo(persp)
dim(Puromycin)
head(Puromycin)
# 對(duì)于狀態(tài)treated，畫出rate關(guān)于cone的散點(diǎn)圖
PuroA = subset(Puromycin, state == "treated")
plot(rate ~ conc, data = PuroA)
有三種方法指明函數(shù)plot()使用的數(shù)據(jù)集
#plot()函數(shù)中使用data選項(xiàng)
? PuroA = subset(Puromycin, state == "treated")
? plot(rate ~ conc, data = PuroA)
#在with()中使用plot()
? with(PuroA, plot(conc, rate))
#使用$直接指向數(shù)據(jù)與變量
? plot(PuroA$rate, PuroA$conc)
# R提供25種不同的符號(hào)和8種不同顏色困食，瀏覽它們的命令是
u = 1:25
plot(u ~ 1, pch = u, col = u, cex = 3)
bg? # 背景色
fg? # 前景色
col # 顏色
col.axis #坐標(biāo)軸顏色
col.lab? #標(biāo)簽顏色
col.main #題目顏色
col.sub? #副題目顏色
和字體相關(guān)的參數(shù)有下面幾個(gè)：
字體形態(tài)
family #全局字體边翁，特指字體的類型，如宋體還是楷體
font? #字體硕盹，特指字體的形態(tài)符匾，如斜體還是粗體
font.axis #坐標(biāo)軸字體
font.lab? #標(biāo)簽字體
font.main #題目字體
font.sub #副題目字體
par(font.axis=1) # 1 普通文本
par(font.lab=2) # 2 粗體
par(font.main=3) # 3 斜體
par(font.sub=4) # 4 粗斜體
字體大小
PS=10 指所有字體
cex指部分
cex.axis #坐標(biāo)軸
cex.lab? #標(biāo)簽
cex.main #題目
cex.sun? #副題目
線條
lty #line style
? 線條的類型：1、2瘩例、3啊胶、4、5垛贤、6? ? 線條的大小可以用lwd調(diào)節(jié)
pch
? 線條符號(hào)種類：0-24 包括正方形焰坪、三角、圓形等24種. 符號(hào)的大小可以用cex調(diào)節(jié)(圖形有多大)聘惦。
# R提供25種不同的符號(hào)和8種不同顏色某饰，瀏覽它們的命令是
u = 1:25
plot(u ~ 1, pch = u, col = u, cex = 1)
# 選擇合適的符號(hào)及其大小與顏色
plot(rate ~ conc, data = PuroA, pch = 2, col = 4, cex = 2.5) 選第二種圖形，第四種顏色部凑，圖形大小為2.5
# 坐標(biāo)軸與標(biāo)題設(shè)定
plot(rate ~ conc, data = PuroA, pch = 2, col = 4,
? ? cex = 2.5, xlim = c(0, 1.2), ylim = c(40, 210),? #xlim 橫軸1~1.2 #ylim 縱軸區(qū)間40~210
? ? ylab = "Concentration",? ? ? ? ? ? ? ? ? ? ? ? ? #xlab 橫軸標(biāo)簽? cex.lab 標(biāo)簽字體大小
? ? xlab = "Rate", cex.lab = 3)
title(main = "Puromycin", cex.main = 3)? ? ? ? ? ? ?
# 連接數(shù)據(jù)點(diǎn)
library(doBy)
PuroA.mean = summaryBy(rate ~ conc, data = PuroA, FUN = mean)
PuroA.mean
conc rate.mean
1 0.02? ? ? 61.5
2 0.06? ? 102.0
3 0.11? ? 131.0
4 0.22? ? 155.5
5 0.56? ? 196.0
6 1.10? ? 203.5
plot(rate ~ conc, data = PuroA, pch = 16, col = 4, cex = 1.5)? #散點(diǎn)圖
points(rate.mean? ~ conc, data = PuroA.mean, col = "cyan", lwd = 10, pch = "x")? #cyan藍(lán)綠色
lines(rate.mean ~ conc, data = PuroA.mean, col = "blue")? ? ? ? #添加趨勢(shì)線
# 下面命令給出兩條光滑曲線
plot(rate ~ conc, data = PuroA)
smooth1 = with(PuroA, lowess(rate ~ conc, f = 0.9))? #f表示曲線的光滑程度 越大轉(zhuǎn)折點(diǎn)越少
smooth2 = with(PuroA, lowess(rate ~ conc, f = 0.3))
lines(smooth1, col = "black")
lines(smooth2, col = "blue")
# 添加多項(xiàng)式擬合線露乏。下面命令給出一次碧浊、二次和三次多項(xiàng)式擬合
plot(rate ~ conc, data = PuroA)
m1 = lm(rate ~ conc, data = PuroA)
m2 = lm(rate ~ conc + I(conc^2), data = PuroA)
m3 = lm(rate ~ conc + I(conc^2) + I(conc^3), data = PuroA)
lines(fitted(m1) ~ conc, data = PuroA, col = "red");lines(fitted(m2) ~ conc, data = PuroA, col = "blue");lines(fitted(m3) ~ conc, data = PuroA, col = "cyan
# 表達(dá)式由對(duì)象和函數(shù)組成涂邀。我們可以用;對(duì)同一行不同的表達(dá)式進(jìn)行分隔
"this expression will be printed";7+13;exp(0+1i*pi)
x <- c(1,2,3,4,5)
x
class(x)
typeof(x)
x[2]="hat"
x
typeof(x)
class(x)
if (x > 1) "orange" else "apple"
typeof(quote(if (x > 1) "orange" else "apple"))
sessionInfo()
sqrt(-1)
`%myop%` <- function(a, b){2*a + 2*b}
1 %myop% 1
x <- 1; y <- 2; z <- 3
function(3)
x=function(3)
x={3,function(),1}
str(x)
a <- rep("a", times=5)
b <- rep("b", times=5)
ifelse(c(TRUE, TRUE, TRUE, FALSE, TRUE), a, b)
example(glm)
help(glm)
?? regression
vignette("affy")
vignette(all=FALSE)，axis(2,tick = FALSE))
rug("deoxy time")
rug("appotosisi",side=3)
title(main="流式結(jié)果",cex.main=2)
lines(c(0:5),col="black",lwd=3,lty=3)
plot(c(0:5), col = "red")
text(4,5, labels = 'dian', font =2)
grid(nx=NA, ny=4, lwd=2, col='skyblue')
par(mfrow=c(3,3),mar=c(2,2,2,2))
dose <- c(20, 30, 40, 45, 60)
drugA <- c(16, 20, 27, 40, 60)
drugB <- c(15, 18, 25, 31, 40)
plot(dose,drugA,type = "b")
opar<-par(no.readonly=TRUE)
par(mfrow=c(2,2))
plot(dose,drugA,type = "b")
par(lty = 2, pch = 17)
plot(dose, drugA, type = "b")
par(opar)
plot(dose,drugA,type="b",main = "第三幅")
a=1:10
dim(a)=c(2,5)
a[2,3]
b=as.data.frame(a)
pheatmap::pheatmap(b)
? ? install.packages("pheatmap")? ?
BiocManager::install("GEOquery")
installed.packages("pheatmap")
a=1:10
? dim(a)=c(2,5)
? ? pheatmap::pheatmap(a)
name<-c("Zhangshan","Lisi","Wangwu","Zhaoliu")
> age<-c(20,30,40,50)
name<-c("Zhangshan","Lisi","Wangwu","Zhaoliu")
age<-c(20,30,40,50)
onedata.frame<-data.frame(name,age)
onedata.frame
attach(onedata.frame)
最后介紹一下hist函數(shù)的返回值
quantile(mpg)
Summary(mpg)
plot(hp, mpg,pch=cyl)
barplot(table(cyl))
hist(cyl)
legend(250, 30,pch=c(4,6,8))
legend=c("4 cylinders", "6cylinders", "8 cylinders"))
text.legend=c("上周pv","本周pv","pv同比增長(zhǎng)","pv環(huán)比增長(zhǎng)")
col2<-c("black","blue")
legend("topleft",pch=c(15,15,16,16),legend=text.legend,col=c(col,col2),bty="n",horiz=TRUE)
tu=par(no.readonly=TRUE)
runif(10,2,3)##生成10個(gè)2到3之間的隨機(jī)數(shù) 前面的r表示隨機(jī)箱锐，unif表示服從均勻分布比勉，
? ? 類似的還有rnorm()，它的功能類似，只是它生成的是服從正態(tài)分布的隨機(jī)數(shù)浩聋，其中norm就是正態(tài)分布的意思观蜗。
# set.seed(1000)
? runif(10,2,3)? 前面用set.seed(1000)限定，可以保證每次取的2到3之間的10個(gè)數(shù)都相同衣洁。
#w<-c(75.0,64.0,47.4,66.9,62.2,62.2,58.7,63.5,66.6,64.0,57.0,69.0,56.9,50.0,72.0)
# quantile(w,probs=seq(0,1,0.1),na.rm=FALSE)
? 0%? 10%? 20%? 30%? 40%? 50%? 60%? 70%? 80%? 90%? 100%
47.40 52.76 56.98 59.40 62.20 63.50 64.00 66.08 67.32 70.80 75.00
