1.## 加載R包
## 下載數(shù)據(jù)，如果文件夾中有會直接讀入
gset = getGEO('GSE32575', destdir=".",getGPL = F)
# 獲取ExpressionSet對象，包括的表達矩陣和分組信息
gset=gset[[1]]
2.#通過pData函數(shù)獲取分組信息
## 獲取分組信息,點開查閱信息
pdata=pData(gset)
## 只要后36個,本次選擇的這36個是配對的。
## 所以，別人的芯片我們也不是全要，我們只要適合自己的數(shù)據(jù)
group_list=c(rep('before',18),rep('after',18))
group_list=factor(group_list)
## 強制限定順序
group_list <- relevel(group_list, ref="before")
3.#通過exprs函數(shù)獲取表達矩陣并校正
exprSet=exprs(gset)
可以先簡單看一下整體樣本的表達情況，
其中exprSet[,-c(1:12)]的意思是去掉這個數(shù)據(jù)的前12個，因為我們需要的后面36個.
boxplot(exprSet[,-c(1:12)],outline=FALSE, notch=T,col=group_list, las=2)
需要人工校正一下，用的方法類似于Quntile Normalization
這里我們用limma包內(nèi)置的一個函數(shù)
library(limma)
exprSet=normalizeBetweenArrays(exprSet)
boxplot(exprSet,outline=FALSE, notch=T,col=group_list, las=2)
4.#判斷是否需要進行數(shù)據(jù)轉換
我們通過分組信息已經(jīng)知道，前12個樣本本次不需要，所以先去掉
exprSet = as.data.frame(exprSet)[,-seq(1,12)]
肉眼看到數(shù)字很大，這時候需要log轉換(選log2)。
此外還有一個腳本可以自動判斷是否需要log轉換，這個腳本來自于GEO2R，我改寫了一點點。
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")}
這個語句會自動判斷，如果已經(jīng)log話就跳過，并提示
"log2 transform not needed"
如果沒有l(wèi)og話，他自動log2，并且提示：
"log2 transform finished"
但是我們發(fā)現(xiàn)一個問題，這個表達數(shù)據(jù)行名是探針，
不是我們熟悉的基因名稱如果TP53，BRCA1這樣的，所以我們需要注釋他。
5.#獲取注釋信息
通常情況下我們使用R包來注釋，R包和平臺的注釋信息對應關系我整理了一個表格 platformMap，
是一個文本文件，在文末有獲取方法。
platformMap <- data.table::fread("platformMap.txt")
我們根據(jù)上述帖子里面提到的platformMap，找到對應的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/")
## 沒有這個包就給你裝，有就不裝
if(!require("illuminaHumanv2.db")) BiocManager::install("illuminaHumanv2.db",update = F,ask = F)
## 加載R包
library(illuminaHumanv2.db)
獲取探針對應的symbol信息
## 獲取表達矩陣的行名，就是探針名稱
probeset <- rownames(exprSet)
## 使用lookup函數(shù)，找到探針在illuminaHumanv2.db中的對應基因名稱
## 如果分析別的芯片數(shù)據(jù)，把illuminaHumanv2.db更換即可
SYMBOL <-? annotate::lookUp(probeset,"illuminaHumanv2.db", "SYMBOL")
## 轉換為向量
symbol = as.vector(unlist(SYMBOL))
制作probe2symbol轉換文件
probe2symbol = data.frame(probeset,symbol,stringsAsFactors = F)
6.探針轉換與基因去重
一個基因會對應對個探針，有些基因名稱會是重復的，這些都需要處理。
對于，多個探針，我們選取在樣本中平均表達量最高的探針作為對應基因的表達量。
一下代碼完成所有事情，而且可以復用。
library(dplyr)
library(tibble)
exprSet <- exprSet %>%
? rownames_to_column(var="probeset") %>%
? #合并探針的信息
? inner_join(probe2symbol_df,by="probeset") %>%
? #去掉多余信息
? select(-probeset) %>%
? #重新排列
? select(symbol,everything()) %>%
? #求出平均數(shù)(這邊的點號代表上一步產(chǎn)出的數(shù)據(jù))
? mutate(rowMean =rowMeans(.[grep("GSM", names(.))])) %>%
? #去除symbol中的NA
? filter(symbol != "NA") %>%
? #把表達量的平均值按從大到小排序
? arrange(desc(rowMean)) %>%
? # symbol留下第一個
? distinct(symbol,.keep_all = T) %>%
? #反向選擇去除rowMean這一列
? select(-rowMean) %>%
? # 列名變成行名
? column_to_rownames(var = "symbol")
7.差異分析
差異分析的矩陣構建有兩種方式
我們選取格式比較簡單的。如果沒有配對信息，差異分析這樣做：
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行。
因為我們這里實際上是有配對信息的，差異分析應該這樣做：
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)
得到的差異分析結果allDiff_pair有7650個，比直接做差異分析要多接近1000個。
8.作圖驗證(非必要)
差異分析的結果，配對和非配對理解起來比較抽象，我們用圖片直觀地感受一下。
首先我們需要得到ggplot2喜歡的數(shù)據(jù)格式，行是觀測，列是變量，這也叫clean data，tide data, 清潔數(shù)據(jù)。
首先基因名稱需要在列，所以需要用t函數(shù)，行列轉置。
data_plot = as.data.frame(t(exprSet))
data_plot = data.frame(pairinfo=rep(1:18,2),
? ? ? ? ? ? ? ? ? ? ? group=group_list,
? ? ? ? ? ? ? ? ? ? ? data_plot,stringsAsFactors = F)
作圖展示:
? 挑選了一個基因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")
看起來雜亂一片，根本得不到有用信息，我們加入他的配對信息，再來作圖：
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")
再來看看真正有差異的基因吧，比如：
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個基因：
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)
## 設定差異基因閾值，減少差異基因用于提取表達矩陣
allDiff_pair=topTable(fit,adjust='BH',coef="group_listafter",number=Inf,p.value=0.05,lfc =0.5)
##提前部分數(shù)據(jù)用作熱圖繪制
heatdata <- exprSet[rownames(allDiff_pair),]
##制作一個分組信息用于注釋
annotation_col <- data.frame(group_list)
rownames(annotation_col) <- colnames(heatdata)
#如果注釋出界，可以通過調(diào)整格子比例和字體修正
pheatmap(heatdata, #熱圖的數(shù)據(jù)
? ? ? ? cluster_rows = TRUE,#行聚類
? ? ? ? cluster_cols = TRUE,#列聚類，可以看出樣本之間的區(qū)分度
? ? ? ? annotation_col =annotation_col, #標注樣本分類
? ? ? ? annotation_legend=TRUE, # 顯示注釋
? ? ? ? show_rownames = F,# 顯示行名
? ? ? ? show_colnames = F,# 顯示列名
? ? ? ? scale = "row", #以行來標準化，這個功能很不錯
? ? ? ? color =colorRampPalette(c("blue", "white","red"))(100))
我們還可以畫一個火山圖
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
# 保存為簡單的文本文件
write.table(d,file="SDDDD.txt")
write.table(d,file="SDDDD.txt",row.names=F,quote = F)
# row.names = F指定行名不寫入文件，quote = F表示變量名不放入雙引號中
# 保存為逗號分割的文本文件
write.csv(d,file = "SDD.csv")
# 保存為R格式文件(直接用save）
save(d,file="SDDDD.RData")
# 在經(jīng)過一段時間的分析后，常需要將工作空間的映像保存，命令為：
save.image()
# 它等價于 save(list = ls(all=TRUE), file = ".RData")
數(shù)據(jù)的讀取
可以使用read.table()、scan()、readfwf()函數(shù)讀入存儲在文本文件中的數(shù)據(jù)
A=read.table("SDDDD.txt",header=TRUE)? # header = TRUE 表明讀入的第一行是變量名
A
# read.table有四個變形
# 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
另一個重要區(qū)別是scan()函數(shù)可以創(chuàng)建不同的對象，在缺省情況下（what被省略），將創(chuàng)建一個數(shù)值型向量。
# 使用函數(shù)read.fwf()，它可以用來讀取文件中一些固定寬度格式的數(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)
# 對于狀態(tài)treated，畫出rate關于cone的散點圖
PuroA = subset(Puromycin, state == "treated")
plot(rate ~ conc, data = PuroA)
有三種方法指明函數(shù)plot()使用的數(shù)據(jù)集
#plot()函數(shù)中使用data選項
? 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種不同的符號和8種不同顏色，瀏覽它們的命令是
u = 1:25
plot(u ~ 1, pch = u, col = u, cex = 3)
bg? # 背景色
fg? # 前景色
col # 顏色
col.axis #坐標軸顏色
col.lab? #標簽顏色
col.main #題目顏色
col.sub? #副題目顏色
和字體相關的參數(shù)有下面幾個：
字體形態(tài)
family #全局字體，特指字體的類型，如宋體還是楷體
font? #字體，特指字體的形態(tài)，如斜體還是粗體
font.axis #坐標軸字體
font.lab? #標簽字體
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 #坐標軸
cex.lab? #標簽
cex.main #題目
cex.sun? #副題目
線條
lty #line style
? 線條的類型：1、2、3、4、5、6? ? 線條的大小可以用lwd調(diào)節(jié)
pch
? 線條符號種類：0-24 包括正方形、三角、圓形等24種. 符號的大小可以用cex調(diào)節(jié)(圖形有多大)。
# R提供25種不同的符號和8種不同顏色，瀏覽它們的命令是
u = 1:25
plot(u ~ 1, pch = u, col = u, cex = 1)
# 選擇合適的符號及其大小與顏色
plot(rate ~ conc, data = PuroA, pch = 2, col = 4, cex = 2.5) 選第二種圖形，第四種顏色，圖形大小為2.5
# 坐標軸與標題設定
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 橫軸標簽? cex.lab 標簽字體大小
? ? xlab = "Rate", cex.lab = 3)
title(main = "Puromycin", cex.main = 3)? ? ? ? ? ? ?
# 連接數(shù)據(jù)點
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)? #散點圖
points(rate.mean? ~ conc, data = PuroA.mean, col = "cyan", lwd = 10, pch = "x")? #cyan藍綠色
lines(rate.mean ~ conc, data = PuroA.mean, col = "blue")? ? ? ? #添加趨勢線
# 下面命令給出兩條光滑曲線
plot(rate ~ conc, data = PuroA)
smooth1 = with(PuroA, lowess(rate ~ conc, f = 0.9))? #f表示曲線的光滑程度 越大轉折點越少
smooth2 = with(PuroA, lowess(rate ~ conc, f = 0.3))
lines(smooth1, col = "black")
lines(smooth2, col = "blue")
# 添加多項式擬合線。下面命令給出一次、二次和三次多項式擬合
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
# 表達式由對象和函數(shù)組成。我們可以用;對同一行不同的表達式進行分隔
"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="流式結果",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同比增長","pv環(huán)比增長")
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個2到3之間的隨機數(shù) 前面的r表示隨機，unif表示服從均勻分布，
? ? 類似的還有rnorm()，它的功能類似，只是它生成的是服從正態(tài)分布的隨機數(shù)，其中norm就是正態(tài)分布的意思。
# set.seed(1000)
? runif(10,2,3)? 前面用set.seed(1000)限定，可以保證每次取的2到3之間的10個數(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
