2019-07-26

1.## 加载R包

## 下载数据，如果文件夹中有会直接读入

gset = getGEO('GSE32575', destdir=".",getGPL = F)

# 获取ExpressionSet对象，包括的表达矩阵和分组信息

gset=gset[[1]]

2.#通过pData函数获取分组信息

## 获取分组信息,点开查阅信息

pdata=pData(gset)

## 只要后36个,本次选择的这36个是配对的。

## 所以，别人的芯片我们也不是全要，我们只要适合自己的数据

group_list=c(rep('before',18),rep('after',18))

group_list=factor(group_list)

## 强制限定顺序

group_list <- relevel(group_list, ref="before")

3.#通过exprs函数获取表达矩阵并校正

exprSet=exprs(gset)

可以先简单看一下整体样本的表达情况，

其中exprSet[,-c(1:12)]的意思是去掉这个数据的前12个，因为我们需要的后面36个.

boxplot(exprSet[,-c(1:12)],outline=FALSE, notch=T,col=group_list, las=2)

需要人工校正一下，用的方法类似于Quntile Normalization

这里我们用limma包内置的一个函数

library(limma)

exprSet=normalizeBetweenArrays(exprSet)

boxplot(exprSet,outline=FALSE, notch=T,col=group_list, las=2)

4.#判断是否需要进行数据转换

我们通过分组信息已经知道，前12个样本本次不需要，所以先去掉

exprSet = as.data.frame(exprSet)[,-seq(1,12)]

肉眼看到数字很大，这时候需要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")}

这个语句会自动判断，如果已经log话就跳过，并提示

"log2 transform not needed"

如果没有log话，他自动log2，并且提示：

"log2 transform finished"

但是我们发现一个问题，这个表达数据行名是探针，

不是我们熟悉的基因名称如果TP53，BRCA1这样的，所以我们需要注释他。

5.#获取注释信息

通常情况下我们使用R包来注释，R包和平台的注释信息对应关系我整理了一个表格 platformMap，

是一个文本文件，在文末有获取方法。

platformMap <- data.table::fread("platformMap.txt")

我们根据上述帖子里面提到的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函数，找到探针在illuminaHumanv2.db中的对应基因名称

## 如果分析别的芯片数据，把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()) %>%

  #求出平均数(这边的点号代表上一步产出的数据)

  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喜欢的数据格式，行是观测，列是变量，这也叫clean data，tide data, 清洁数据。

首先基因名称需要在列，所以需要用t函数，行列转置。

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.后续分析

差异分析后还有很多操作，

画一张热图

library(pheatmap)

## 设定差异基因阈值，减少差异基因用于提取表达矩阵

allDiff_pair=topTable(fit,adjust='BH',coef="group_listafter",number=Inf,p.value=0.05,lfc =0.5)

##提前部分数据用作热图绘制

heatdata <- exprSet[rownames(allDiff_pair),]

##制作一个分组信息用于注释

annotation_col <- data.frame(group_list)

rownames(annotation_col) <- colnames(heatdata)

#如果注释出界，可以通过调整格子比例和字体修正

pheatmap(heatdata, #热图的数据

        cluster_rows = TRUE,#行聚类

        cluster_cols = TRUE,#列聚类，可以看出样本之间的区分度

        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")

# 在经过一段时间的分析后，常需要将工作空间的映像保存，命令为：

save.image()

# 它等价于 save(list = ls(all=TRUE), file = ".RData")

数据的读取

可以使用read.table()、scan()、readfwf()函数读入存储在文本文件中的数据

A=read.table("SDDDD.txt",header=TRUE)  # header = TRUE 表明读入的第一行是变量名

A

# read.table有四个变形

# read.csv(), read.csv2(), read.delim(), read.delim2()

# 使用函数scan()，它比read.table()更加灵活，唯一的区别是scan()可以指定变量的类型

mydata=scan("SDD.csv",what=list(obs=0,treat="",weight=0),sep=",",skip = 1)

mydata

另一个重要区别是scan()函数可以创建不同的对象，在缺省情况下（what被省略），将创建一个数值型向量。

# 使用函数read.fwf()，它可以用来读取文件中一些固定宽度格式的数据。

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)

# 对于状态treated，画出rate关于cone的散点图

PuroA = subset(Puromycin, state == "treated")

plot(rate ~ conc, data = PuroA)

有三种方法指明函数plot()使用的数据集

#plot()函数中使用data选项

  PuroA = subset(Puromycin, state == "treated")

  plot(rate ~ conc, data = PuroA)

#在with()中使用plot()

  with(PuroA, plot(conc, rate))

#使用$直接指向数据与变量

  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  #副题目颜色

和字体相关的参数有下面几个：

字体形态

family #全局字体，特指字体的类型，如宋体还是楷体

font  #字体，特指字体的形态，如斜体还是粗体

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调节

pch

  线条符号种类：0-24 包括正方形、三角、圆形等24种. 符号的大小可以用cex调节(图形有多大)。

# 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 纵轴区间40~210

    ylab = "Concentration",                          #xlab 横轴标签  cex.lab 标签字体大小

    xlab = "Rate", cex.lab = 3)

title(main = "Puromycin", cex.main = 3)             

# 连接数据点

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

# 表达式由对象和函数组成。我们可以用;对同一行不同的表达式进行分隔

"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函数的返回值

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环比增长")

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之间的随机数 前面的r表示随机，unif表示服从均匀分布，

    类似的还有rnorm()，它的功能类似，只是它生成的是服从正态分布的随机数，其中norm就是正态分布的意思。

# set.seed(1000)

  runif(10,2,3)  前面用set.seed(1000)限定，可以保证每次取的2到3之间的10个数都相同。

#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