20个R语言习题

发表于 | 分类于
| 字数统计: 5,729 | 阅读时长 ≈ 26

原题目

生信人的20个R语言习题。

题目代码

1. 安装R语言包

安装以下包：

1
2
3
4
数据包： ALL, CLL, pasilla, airway 
软件包：limma，DESeq2，clusterProfiler 
工具包：reshape2
绘图包：ggplot2

代码如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
source("http://bioconductor.org/biocLite.R")
biocLite("ALL")
biocLite("CLL")
biocLite("pasilla")
biocLite("airway")
biocLite("limma")
biocLite("DESeq2")
biocLite("clusterProfiler ")
biocLite("airway")
library(ALL)
library(CLL)
library(pasilla)
library(airway)
library(DESeq2)
library(clusterProfiler )
library(reshape2)
library(ggplot2)

2. ExpressionSet对象

注：补充ExpressionSet知识

CLL包中里面有一个sCLLex数据，它是ExpressionSet格式的，如下所示：

1
2
3
4
5
data("sCLLex")
class(sCLLex)
# [1] "ExpressionSet"
# attr(,"package")
# [1] "Biobase"

现在查看一下这个对象，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
> sCLLex
ExpressionSet (storageMode: lockedEnvironment)
assayData: 12625 features, 22 samples 
  element names: exprs 
protocolData: none
phenoData
  sampleNames: CLL11.CEL CLL12.CEL ... CLL9.CEL (22 total)
  varLabels: SampleID Disease
  varMetadata: labelDescription
featureData: none
experimentData: use 'experimentData(object)'
Annotation: hgu95av2

看一下这个数据的文档，可以得到这些信息：

这个数据集中有22个样本，12625个基因；使用exprs可以查看这个数据集的表达水平；使用phenoData可以查看这个数据集的分组信息，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
data_expression <- exprs(sCLLex) 
#提取出来表达矩阵，赋给data_expression
head(data_expression)
# CLL11.CEL CLL12.CEL CLL13.CEL CLL14.CEL CLL15.CEL CLL16.CEL CLL17.CEL
# 1000_at    5.743132  6.219412  5.523328  5.340477  5.229904  4.920686  5.325348
# 1001_at    2.285143  2.291229  2.287986  2.295313  2.662170  2.278040  2.350796
...
phenoData(sCLLex)
# 查看分组信息
# An object of class 'AnnotatedDataFrame'
# sampleNames: CLL11.CEL CLL12.CEL ... CLL9.CEL (22 total)
# varLabels: SampleID Disease
# varMetadata: labelDescription

使用View(data_expression)可以直观看出来这个表达矩阵，如下所示：

mark

再来查看一个sCCLex这个数据集的其它信息，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
sampleNames(sCLLex)
#查看样本编号
# [1] "CLL11.CEL" "CLL12.CEL" "CLL13.CEL" "CLL14.CEL" "CLL15.CEL" "CLL16.CEL"
# [7] "CLL17.CEL" "CLL18.CEL" "CLL19.CEL" "CLL20.CEL" "CLL21.CEL" "CLL22.CEL"
# [13] "CLL23.CEL" "CLL24.CEL" "CLL2.CEL"  "CLL3.CEL"  "CLL4.CEL"  "CLL5.CEL" 
# [19] "CLL6.CEL"  "CLL7.CEL"  "CLL8.CEL"  "CLL9.CEL" 
varMetadata(sCLLex) 
#查看所有表型变量
# labelDescription
# SampleID          Sample ID
# Disease  Stable/Progressive
data_phenotype=pData(sCLLex) 
#提取表型信息
# data_phenotype
# SampleID  Disease
# CLL11.CEL    CLL11 progres.
# CLL12.CEL    CLL12   stable
# CLL13.CEL    CLL13 progres.
# CLL14.CEL    CLL14 progres.
# CLL15.CEL    CLL15 progres.
# CLL16.CEL    CLL16 progres.
# CLL17.CEL    CLL17   stable
# CLL18.CEL    CLL18   stable
# CLL19.CEL    CLL19 progres.
# CLL20.CEL    CLL20   stable
# CLL21.CEL    CLL21 progres.
# CLL22.CEL    CLL22   stable
# CLL23.CEL    CLL23 progres.
# CLL24.CEL    CLL24   stable
# CLL2.CEL      CLL2   stable
# CLL3.CEL      CLL3 progres.
# CLL4.CEL      CLL4 progres.
# CLL5.CEL      CLL5 progres.
# CLL6.CEL      CLL6 progres.
# CLL7.CEL      CLL7 progres.
# CLL8.CEL      CLL8 progres.
# CLL9.CEL      CLL9   stable
featureNames(sCLLex)[1:10] 
#查看基因芯片编码
# [1] "1000_at"   "1001_at"   "1002_f_at" "1003_s_at" "1004_at"   "1005_at"  
# [7] "1006_at"   "1007_s_at" "1008_f_at" "1009_at" 
featureNames(sCLLex) %>% unique() %>% length() 
# 查看是否有重复的编码
# [1] 12625

上面的%>%这个符号相当于Linux中的管道命令。

接着，提取表型信息，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
pdata <- pData(sCLLex)
head(pdata)
# SampleID  Disease
# CLL11.CEL    CLL11 progres.
# CLL12.CEL    CLL12   stable
# CLL13.CEL    CLL13 progres.
# CLL14.CEL    CLL14 progres.
# CLL15.CEL    CLL15 progres.
# CLL16.CEL    CLL16 progres.
group_list=as.character(pdata$Disease) 
#从数据框中只要表型信息
table(group_list) 
#统计表型信息
# group_list
# progres.   stable 
# 14        8

绘制芯片数据的质量值，如下所示：

1
2
3
4
5
6
y <- melt(as.data.frame(data_expression))
p <- ggplot(data=y,aes(x=variable,y=value))
p <- p + geom_boxplot(aes(fill=variable))
p <- p + theme(axis.text.x = element_text(angle=45,hjust=1,vjust=1))
p <- p + xlab("分组信息") + ylab("表达值") + guides(fill=FALSE)
p

3.了解str,head,help函数

略，第2题中已经有了。

4.了解并安装hgu95av2.db

注：补充db数据格式的相关知识

hgu95av2.db是一个注释包，它为hgu95av2平台的芯片提供注释，这个包中有很多注释文件，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
ls("package:hgu95av2.db")
# [1] "hgu95av2"              "hgu95av2.db"           "hgu95av2_dbconn"      
# [4] "hgu95av2_dbfile"       "hgu95av2_dbInfo"       "hgu95av2_dbschema"    
# [7] "hgu95av2ACCNUM"        "hgu95av2ALIAS2PROBE"   "hgu95av2CHR"          
# [10] "hgu95av2CHRLENGTHS"    "hgu95av2CHRLOC"        "hgu95av2CHRLOCEND"    
# [13] "hgu95av2ENSEMBL"       "hgu95av2ENSEMBL2PROBE" "hgu95av2ENTREZID"     
# [16] "hgu95av2ENZYME"        "hgu95av2ENZYME2PROBE"  "hgu95av2GENENAME"     
# [19] "hgu95av2GO"            "hgu95av2GO2ALLPROBES"  "hgu95av2GO2PROBE"     
# [22] "hgu95av2MAP"           "hgu95av2MAPCOUNTS"     "hgu95av2OMIM"         
# [25] "hgu95av2ORGANISM"      "hgu95av2ORGPKG"        "hgu95av2PATH"         
# [28] "hgu95av2PATH2PROBE"    "hgu95av2PFAM"          "hgu95av2PMID"         
# [31] "hgu95av2PMID2PROBE"    "hgu95av2PROSITE"       "hgu95av2REFSEQ"       
# [34] "hgu95av2SYMBOL"        "hgu95av2UNIGENE"       "hgu95av2UNIPROT"

安装hgu95av2.db，如下所示：

1
2
biocLite("hgu95av2.db")
library(hgu95av2.db)

5.理解toTable的用法

原是描述：理解 head(toTable(hgu95av2SYMBOL)) 的用法，找到 TP53 基因对应的探针ID。

先看一下hgu95av2SYMBOL是什么，如下所示：

1
2
?hgu95av2SYMBOL
# hgu95av2SYMBOL is an R object that provides mappings between manufacturer identifiers and gene abbreviations.

大致意思就是说，hgu95av2SYMBOL是一个R对象，它描述厂家的标记符(identifiers)与基因缩写之间的映射关系。

再看一下toTable的用法，如下所示：

1
2
?toTable
These methods are part of the Bimap interface (see ?Bimap for a quick overview of the Bimap objects and their interface).

大致意思是说，toTable是一种能够以数据框的形式来操作一个Bimap对象的方法，也就是把Bimap对象转换为一个数据框，这些方法是Bimap interface方法的一部分。Bimap指的是一种映射关系，例如探针的编号与基因名称之间的映射，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
summary(hgu95av2SYMBOL)
# 查看hgu95av2SYMBOL大致信息
# SYMBOL map for chip hgu95av2 (object of class "ProbeAnnDbBimap")
# |
#   | Lkeyname: probe_id (Ltablename: probes)
# |    Lkeys: "1000_at", "1001_at", ... (total=12625/mapped=11471)
# |
#   | Rkeyname: symbol (Rtablename: gene_info)
# |    Rkeys: "A1BG", "A2M", ... (total=59682/mapped=8595)
# |
#   | direction: L --> R
gene_id <- toTable(hgu95av2SYMBOL)
head(gene_id) #这是一个数据框
# probe_id  symbol
# 1   1000_at   MAPK3
# 2   1001_at    TIE1
# 3 1002_f_at CYP2C19
# 4 1003_s_at   CXCR5
# 5   1004_at   CXCR5
# 6   1005_at   DUSP1
str(gene_id)
# 'data.frame':	11471 obs. of  2 variables:
#   $ probe_id: chr  "1000_at" "1001_at" "1002_f_at" "1003_s_at" ...
# $ symbol  : chr  "MAPK3" "TIE1" "CYP2C19" "CXCR5" ...
filter(gene_id,symbol=="TP53")
# 找到 TP53 基因对应的探针ID
# probe_id symbol
# 1   1939_at   TP53
# 2 1974_s_at   TP53
# 3  31618_at   TP53
library(dplyr)
filter(gene_id,symbol=="TP53")
# 找到 TP53 基因对应的探针ID
# probe_id symbol
# 1   1939_at   TP53
# 2 1974_s_at   TP53
# 3  31618_at   TP53

6.探针与基因的关系

原题目描述：

理解探针与基因的对应关系，总共多少个基因，基因最多对应多少个探针，是哪些基因，是不是因为这些基因很长，所以在其上面设计多个探针呢？

答：不管是Agilent芯片，还是Affymetrix芯片，上面设计的探针都非常短。最长的如Agilent芯片上的探针，往往都是60bp，但是往往一个基因的长度都好几Kb。因此一般多个探针对应一个基因，取最大表达值探针来作为基因的表达量。

7.提取目的基因

原题目描述：

第二步提取到的表达矩阵是12625个探针在22个样本的表达量矩阵，找到那些不在 hgu95av2.db 包收录的对应着SYMBOL的探针。

芯片表达的数据存放在data_expression这个变量中，它有22个样本的12625个探针数据，如下所示：

1
2
3
4
5
6
7
class(data_expression)
# [1] "matrix"
str(data_expression)
# num [1:12625, 1:22] 5.74 2.29 3.31 1.09 7.54 ...
# - attr(*, "dimnames")=List of 2
# ..$ : chr [1:12625] "1000_at" "1001_at" "1002_f_at" "1003_s_at" ...
# ..$ : chr [1:22] "CLL11.CEL" "CLL12.CEL" "CLL13.CEL" "CLL14.CEL" ...

现在要解决的一个问题是：在data_expression中有一些基因没有收录在hug95av2.db包中，现在要找到这些基因，把它们都去掉。

思路：找到hug95av2.db包中能够与基因映射起来的探针，再找到这些探针对应的基因SYMBOL，然后在data_expression中取出这些基因即可。此时要使用到hug95av2.db包中的hug95av2SYMBOL这个文件，如下所示：

1
?hug95av2SYMBOL

查阅其文档，了解到如下信息：hug95av2SYMBOL是一个R对象，它提供的是芯片生产厂家与基因缩写之间的映射信息。这个映射的信息主要依据Entrez Gene数据库。现在我们通过mappedkeys()这个函数，得到映射到基因上的探针信息，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
probe_map <- hgu95av2SYMBOL
length(probe_map) 
#全部的探针数目
# [1] 12625
probe_info <- mappedkeys(probe_map)
length(probe_info) 
#探针与基因产生映射的数目
# [1] 11471
gene_info <- as.list(probe_map[probe_info])
# 转化为数据表
length(gene_info)
# [1] 11471
gene_symbol <- toTable(probe_map[probe_info])
# 从hgu95av2SYMBOL文件中，取出有映射关系的探针，并生成数据框给gene_symbol
head(gene_symbol)
# probe_id  symbol
# 1   1000_at   MAPK3
# 2   1001_at    TIE1
# 3 1002_f_at CYP2C19
# 4 1003_s_at   CXCR5
# 5   1004_at   CXCR5
# 6   1005_at   DUSP1

其中，gene_symbol就是有映射关系的基因，它里面的数据是探针的编码以及探针对应的基因，在data_expression中它们挑出来即可。这里使用到了mappedkeys()函数，这个函数用于处理映射文件(bimap)，它的参数就是一个Bimap对象，例如上面的hgu95av2SYMBOL这个对象，关于Bimap对象，可以看这篇笔记《Bimap对象笔记》。

8. 过滤表达矩阵

原问题描述：过滤表达矩阵，删除那1165个没有对应基因名字的探针。这时再看一下原始数据，如下所示：

1
2
3
4
5
6
7
8
head(data_expression)
# CLL11.CEL CLL12.CEL CLL13.CEL CLL14.CEL CLL15.CEL CLL16.CEL CLL17.CEL
# 1000_at    5.743132  6.219412  5.523328  5.340477  5.229904  4.920686  5.325348
# 1001_at    2.285143  2.291229  2.287986  2.295313  2.662170  2.278040  2.350796
# 1002_f_at  3.309294  3.318466  3.354423  3.327130  3.365113  3.568353  3.502440
# 1003_s_at  1.085264  1.117288  1.084010  1.103217  1.074243  1.073097  1.091264
# 1004_at    7.544884  7.671801  7.474025  7.152482  6.902932  7.368660  6.456285
# 1005_at    5.083793  7.610593  7.631311  6.518594  5.059087  4.855161  5.176975

思路：在这个原始数据中，有12625个探针，下面通过isNA函数可以得到一个逻辑向量，其中TRUE表示没有映射关系的探针，FALSE表示有映射关系的探针。然后使用seq生成一个整数向量，并提取那些有映射关系的探针所在的列，提取原始数据即可，如下所示：data_filter <- data_expression[rownames(data_expression)%in%probe_info]

1
2
3
4
5
6
7
8
9
10
11
12
13
14
data_filter <- data_expression[rownames(data_expression)%in%probe_info,]
# rownames(data_expression)是原始数据表达值的名称，这里提探针的编号
# probe_info是有映射的探针的编号
# A%in%B，表示A中有哪些成员在B中，有的返回TRUE，无的返回FALSE
# 这行命令表示，提取那些在probe_info中的探针的表达值，赋值给变量data_filter
table(rownames(data_expression) %in% probe_info)
# FALSE  TRUE 
# 1154 11471 
# 从上面的数据可以看出来，原始表达值中有11471个数据在probe_info中
table(rownames(data_expression) %in% rownames(data_filter))
# FALSE  TRUE 
# 1154 11471

现在就得了过滤后的表达矩阵，即data_filter。

9.整合表达矩阵

原问题描述：整合表达矩阵，多个探针对应一个基因的情况下，只保留在所有样本里面平均表达量最大的那个探针。

思路：基因芯片的的多个探针会针对同一个基因，因此会有重复，先看一下表达矩阵data_filter的样子：

mark

其中列是不同的样本，行是不同的探针，此时我们解决问题的思路应该是这个样子的：第一，把表达矩阵中的探针名称转换为基因名称；第二，计算每行的均值；第三，在相同基因的表达值中，取最大的那行。

第一步，要使用到的数据是gene_symbol，它的内容是，去掉了无映射关系的探针编号后，剩余的探针与相应基因的映射关系，如下所示：

1
2
3
4
5
6
7
8
head(gene_symbol)
# probe_id  symbol
# 1   1000_at   MAPK3
# 2   1001_at    TIE1
# 3 1002_f_at CYP2C19
# 4 1003_s_at   CXCR5
# 5   1004_at   CXCR5
# 6   1005_at   DUSP1

现在将表达矩阵data_filter中的探针名称排序，如下所示：

1
2
3
4
5
6
7
8
data_filter_order <- data_filter[order(rownames(data_filter)),]
# 对过滤后的表达值的数据根据探针名称进行排序
gene_symbol_order <- gene_symbol[order(gene_symbol$probe_id),]
# 对gene_symbol按照探针名称进行排序
table(gene_symbol_order$probe_id == rownames(data_filter_order))
# 查看地表达值的矩阵的名称是否与gene_symbol中的探针名称是否一一对应
# TRUE 
# 1147

排序的目的在于，将探针与基因名一一对应起来，排序后的表达矩阵如下所示：

mark

第二步，计算每行的均值，使用apply函数即可，如下所示：

1
2
3
4
row_mean <- apply(data_filter_order,1,mean)
head(row_mean)
#  100_g_at   1000_at   1001_at 1002_f_at 1003_s_at   1004_at 
#  2.289762  5.456561  2.334098  3.402072  1.126168  7.439034

第三步，保留相同基因中表达值最高的那一行，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
index_probe <- function(x){
  names(x)[which.max(x)]
}
# 构建函数
max_index <- tapply(row_mean,gene_symbol$symbol,index_probe)
# index_probe这个函数的功能在于，用gene_symbol$symol（基因名）对row_mean进行分组，然后挑出相同基因中表达值最大的那个值所对应的探针编号
data_filter_order_unique <- data_filter_order[(gene_symbol_order$probe_id %in% max_index),]
# gene_symbol_order与data_filter_order都按照探针名称进行排序，再挑出max_index中的探针哪些在gene_symbol_order$probe_id中
unique_map_gene_probe <- gene_symbol_order[(gene_symbol_order$probe_id %in% rownames(data_filter_order_unique)),]
# 将data_filter_order_unique那些与gene_symbol_order$probe_id对应起来的探针与基因名提取出来
#table(unique_map_gene_probe$probe_id == rownames(data_filter_order_unique))
TRUE 
8595 
# 可以发现，unique_map_gene_probe$probe_id中的探针顺序与data_filter_order_unique中的探针顺序是致的
rownames(data_filter_order_unique) <- unique_map_gene_probe$symbol
# 将最终过滤后的矩阵的名称换为基因名

mark

10.更改矩阵名称

上题已经实现。

11.绘制基本图形

原题描述：对第10步得到的表达矩阵进行探索，先画第一个样本的所有基因的表达量的boxplot,hist,density ， 然后画所有样本的 这些图

思路：使用ggplot2绘图，要把数据转换为ggplot2可识别的形式，也就是第1列是基因名，第2列是分组信息，第3列是表达值，如下所示：

1
2
3
4
5
library(ggplot2)
library(reshape2)
data_draw <- melt(data_filter_order_unique)
colnames(data_draw) <- c("Gene","Group","Value")
data_draw$status <- rep(data_phenotype$Disease,each=nrow(data_filter_order_unique))

箱线图

1
2
3
p <- ggplot(data=data_draw,aes(x=Group,y=Value,fill=status))+ geom_boxplot()
p <- p + theme(axis.text.x = element_text(angle=45,hjust=1,vjust=1))
p

mark

小提琴图

1
2
3
p <- ggplot(data=data_draw,aes(x=Group,y=Value,fill=status))+ geom_violin()
p <- p + theme(axis.text.x = element_text(angle=45,hjust=1,vjust=1))
p

mark

密度图

1
2
3
p <- ggplot(data=data_draw,aes(x=Value,fill=status))+ geom_histogram(bins=500)
p <- p + facet_wrap(~Group,nrow=5)
p

mark

密度图

1
2
3
p <- ggplot(data=data_draw,aes(x=Value,fill=status))+ geom_density()
p <- p + facet_wrap(~Group,nrow=5)
p

mark

PCA图

1
2
3
4
5
6
pca_data <- prcomp(t(data_expression),scale=TRUE)
pcx <- data.frame(pca_data$x)
pcr <- cbind(samples=rownames(pcx),group_list, pcx) 
p  <- ggplot(pcr, aes(PC1, PC2))+geom_point(size=5, aes(color=group_list)) +
  geom_text(aes(label=samples),hjust=-0.1, vjust=-0.3)
p

mark

聚类图

聚类的相关知识参考这篇笔记《聚类分析笔记》。

1
2
3
4
5
6
7
8
9
data_clust <- t(data_expression)
rownames(data_clust) <- colnames(data_expression)
data_clust_dist <- dist(data_clust,method="euclidean")
hc <-hclust(data_clust_dist,"ward")
library(factoextra)
fviz_dend(hc, k = 4,
          cex = 0.5,
          k_colors = c("#2E9FDF", "#00AFBB", "#E7B800", "#FC4E07"),
          color_labels_by_k = TRUE, rect = TRUE)

mark

12.统计学指标

原题描述：理解统计学指标mean,median,max,min,sd,var,mad并计算出每个基因在所有样本的这些统计学指标，最后按照mad值排序，取top 50 mad值的基因，得到列表。

mean,median,max,min,sd,var,mad中除了mad外，其余的好理解，下面只介绍一下mad。

mad指的是绝对中位数差(median absolute deviation)，在R中的函数是mad()。它的用法如下所示：

1
mad(x, center = median(x), constant = 1.4826, na.rm = FALSE, low = FALSE, high = FALSE)

其中，x数值向量；center可选，是中心点，默认是中位数；constant是缩放因子(scale factor)，na.rm如果是TRUE，在计算前将x中的NA删除，low如果是TRUE，计算lo-median，也就是说，对于个数为偶数的样本，在最后计算中位数时不使用两个中间值的均值，而是采用其中较小的值。high如果为TRUE，计算‘hi-median’，也就是对于偶数样本，采用两个中间值的较大者作为中位数。

在R中，mad的计算公式为constant * cMedian(abs(x - center))，其中，center默认为median(x)，cMedian为“低”中位数或“高”中位数，constant默认值为1.4826，它近似等于$1/\Phi^{-1}(3/4)=1/qnorm(3/4)$，这样就可以确保对于服从$N(\mu,\sigma^2)$分布的$X_{i}$和大的n值，存在以下一致性：

$E[mad(X_{1},…,X_{n})]=\sigma$

这里取3/4是因为MAD包含标准正态累积分布函数的50%，具体的推导过程可以看这个帖子《绝对中位差Median Absolute Deviation》。

如果na.rm为TRUE，在执行计算之前，将会删除x中的NA值；否则，只要x中存在任何NA值，mad函数将会返回NA

我们来看一个案例：

1
2
3
4
5
6
7
8
9
10
11
12
x <- c(1,2,3,5,7,8)
x_median <- median(x)
x_median
# [1] 4
y <- median(abs(x - x_median))
y
# [1] 2.5
mad(x)
# [1] 3.7065
# 这个值也就是y乘以1.4826
mad(x)/y
# [1] 1.4826

13.提取 TOP 50 mad值的基因列表

原是描述：根据第12步骤得到top 50 mad值的基因列表来取表达矩阵的子集，并且热图可视化子表达矩阵。试试看其它5种热图的包的不同效果，先看一下原表，如下所示：

mark

从上面的这个表可以知道，行是基因名，列是样本名，现在进行计算，如下所示：

注：补充apply系列函数的相关知识

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
data_mad <- apply(data_filter_order_unique,1,mad)
head(data_mad)
# RABGGTA      MAPK3    CYP2C19      CXCR5      DUSP1    EIF2AK2 
# 0.03886242 0.22779776 0.08294023 0.38497146 1.43854685 0.24478690 
data_mad_30 <- head(sort(data_mad,decreasing = TRUE),30)
data_mad_30
# FAM30A  IGF2BP3      DMD     TCF7   SLAMF1      FOS      HBB   LGALS1    IGLC1 
# 3.982191 3.234011 3.071541 2.993352 2.944105 2.938078 2.732776 2.600604 2.590895 
# ZAP70     FCN1   LHFPL2      HBB   S100A8  GUSBP11   COBLL1    VIPR1    PCDH9 
# 2.579046 2.371590 2.317045 2.267136 2.220970 2.204212 2.179666 2.171912 2.144223 
# SLC25A1     JUND      IGH     H1FX  ZNF804A     OAS1     CCL3     GNLY     CYBB 
# 2.115238 2.105234 2.090866 2.052788 1.986163 1.883509 1.862311 1.814364 1.811973 
# VAMP5   RNASE6     RGS2 
# 1.791017 1.775430 1.770151 
data_mad_30_names <- names(data_mad_30)
data_mad_30_names
# [1] "FAM30A"  "IGF2BP3" "DMD"     "TCF7"    "SLAMF1"  "FOS"     "HBB"     "LGALS1" 
# [9] "IGLC1"   "ZAP70"   "FCN1"    "LHFPL2"  "HBB"     "S100A8"  "GUSBP11" "COBLL1" 
# [17] "VIPR1"   "PCDH9"   "SLC25A1" "JUND"    "IGH"     "H1FX"    "ZNF804A" "OAS1"   
# [25] "CCL3"    "GNLY"    "CYBB"    "VAMP5"   "RNASE6"  "RGS2" 
data_mad_30_heatmap <- data_filter_order_unique[data_mad_30_names,]
# 提取mad 30的子集

以下为提取后的数据集，这个数据集是一个矩阵，行是基因名，列是样本名，如下所示：

mark

绘制这个数据集的热图，目前已经绘制热图的函数有：

heatmap

heatmap(stat)这个是stat包中自带的热图工具，用法如下所示：

1
2
3
4
5
6
7
8
9
10
heatmap(x, Rowv = NULL, Colv = if(symm)"Rowv" else NULL,
        distfun = dist, hclustfun = hclust,
        reorderfun = function(d, w) reorder(d, w),
        add.expr, symm = FALSE, revC = identical(Colv, "Rowv"),
        scale = c("row", "column", "none"), na.rm = TRUE,
        margins = c(5, 5), ColSideColors, RowSideColors,
        cexRow = 0.2 + 1/log10(nr), cexCol = 0.2 + 1/log10(nc),
        labRow = NULL, labCol = NULL, main = NULL,
        xlab = NULL, ylab = NULL,
        keep.dendro = FALSE, verbose = getOption("verbose"), ...)

如下所示：

1
2
3
4
heatmap(data_mad_30_heatmap,scale=c("row"))
# 通常我们要查看某个基因在不同样本中的表达情况，就采用row进行均一化
# 上面的这个数据行是基因名，列是样本名，因此采用行均一化，就可以看相同的某个基因在不同样本中的表达情况
# 如果要想看一个样本中不同基因的表达情况，就要采用column进行均一化

mark

heatmap.2

heatmap.2是gplots包中的函数（文档太长，就不列了），绘制热图如下所示：

1
2
library(gplots)
heatmap.2(data_mad_30_heatmap,scale="row")

mark

d3heatmap

d3heatmap包可用于生成交互式热图绘制,可通过以下代码生成： 

1
2
3
4
5
6
if (!require("devtools"))
install.packages("devtools")
devtools::install_github("rstudio/d3heatmap")
library(d3heatmap)
d3heatmap(data_mad_30_heatmap, colors = "RdBu", k_row = 4, k_col = 2)
# colors是颜色，k_row是行的聚类数目，用不同的颜色标记了出来，k_col是列的聚类数目

mark

Heatmap

这个函数是ComplexHeatmap包中的函数。用法如下所示：

1
2
3
4
install.packages("ComplexHeatmap")
biocLite("ComplexHeatmap")
library(ComplexHeatmap)
Heatmap(data_mad_30_heatmap,name="Mad 30")

mark

pheatmap

此函数是pheatmap包中的函数，用法如下所示：

1
2
library(pheatmap)
pheatmap(data_mad_30_heatmap)

mark

上述几种绘制热图的函数是比较常用的，还有很多参数，需要自己调一下。

14.不同统计指标的交集

原题描述：取不同统计学指标mean,median,max,mean,sd,var,mad的各top50基因列表，使用UpSetR包来看他们之间的overlap情况。先计算一下各个统计学指标，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
data_mean = tail(sort(apply(data_filter_order_unique,1,mean)),30)
data_median = tail(sort(apply(data_filter_order_unique,1,median)), 30)
data_max <- tail(sort(apply(data_filter_order_unique,1,max)),30)
data_min <- tail(sort(apply(data_filter_order_unique,1,min)),30)
data_sd <- tail(sort(apply(data_filter_order_unique,1,sd)),30)
data_var <- tail(sort(apply(data_filter_order_unique,1,var)),30)
data_mad <- tail(sort(apply(data_filter_order_unique,1,mad)),30) 
data_all <- unique(c(names(data_mean),names(data_median),names(data_max),names(data_min),
              names(data_sd),names(data_var),names(data_mad)))
# 取并集
data_upset <- data.frame(data_all,
                         data_mean=ifelse(data_all %in% names(data_mean) ,1,0),
                         data_median=ifelse(e_all %in% names(data_median) ,1,0),
                         data_max=ifelse(e_all %in% names(data_max) ,1,0),
                         data_min=ifelse(e_all %in% names(data_min) ,1,0),
                         data_sd=ifelse(e_all %in% names(data_sd) ,1,0),
                         data_var=ifelse(e_all %in% names(data_var) ,1,0),
                         data_mad=ifelse(e_all %in% names(data_mad) ,1,0))
# 生成一个数据框，内容是不同的统计量是在并集中的重合，如果在为1，不在为0
install.packages("UpSetR")
library("UpSetR")
upset(data_upset,nsets = 7,sets.bar.color = "#56B4E9")

mark

15.提取表型数据

已经在前面完成。

16. 聚类

已经在第11题中完成。

17.PCA

已经在第11题中完成。

18.批量t检验

注：这里要补充一下统计学的相关知识。

代码如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
group_list
# [1] "progres." "stable"   "progres." "progres." "progres." "progres." "stable"  
# [8] "stable"   "progres." "stable"   "progres." "stable"   "progres." "stable"  
# [15] "stable"   "progres." "progres." "progres." "progres." "progres." "progres."
# [22] "stable"  
gl=as.factor(group_list) 
#将表型数据因子化
group1 = which(group_list == levels(gl)[1]) 
# levels(group_list)[1]返回第一个因子progres
# 从group_list中选出progres的元素
# 用which来获取他们所在的索引
group2 = which(group_list == levels(gl)[2]) 
#返回第二个因子stable
data_t1 = data_expression[, group1] 
#将表型为progres的样本选出来
data_t2 = data_expression[, group2] 
#将表型为stable的样本选出来
data_t = cbind(data_t1, data_t2) 
pvals = apply(data_expression, 1, function(x){
  t.test(as.numeric(x)~group_list)$p.value
})
p.adj = p.adjust(pvals, method = "BH")  
data_mean_1 = rowMeans(data_t1) 
#progres是对照组
data_mean_2 = rowMeans(data_t2) 
#stable是使用药物处理后的——处理组
logFC = data_mean_2-data_mean_1
DEG_t.test = cbind(data_mean_1, data_mean_2, log2FC, pvals, p.adj)
DEG_t.test=DEG_t.test[order(DEG_t.test[,4]),] #从小到大排序
DEG_t.test=as.data.frame(DEG_t.test)
head(DEG_t.test)
# data_mean_1 data_mean_2     log2FC        pvals     p.adj
# 36129_at    7.875615    8.791753  0.9161377 1.629755e-05 0.2057566
# 37676_at    6.622749    7.965007  1.3422581 4.058944e-05 0.2436177
# 33791_at    7.616197    5.786041 -1.8301554 6.965416e-05 0.2436177
# 39967_at    4.456446    2.152471 -2.3039752 8.993339e-05 0.2436177
# 34594_at    5.988866    7.058738  1.0698718 9.648226e-05 0.2436177
# 32198_at    4.157971    3.407405 -0.7505660 2.454557e-04 0.3516678

19.limma

原题描述：使用limma包对表达矩阵及样本分组信息进行差异分析，得到差异分析表格，重点看logFC和P值，画个火山图(就是logFC和-log10(P值)的散点图。)。

注：这里需要补充一下基因芯片原理的相关知识。

limma包是常用的分析芯片数据的R包，分析过程如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
design1=model.matrix(~factor(group_list)) 
# 这是一个设计矩阵
colnames(design1)=levels(factor(group_list))
rownames(design1)=colnames(data_expression)
fit1 = lmFit(data_expression,design1)
fit1=eBayes(fit1)
options(digits = 3) #设置结果的小数位数为3
mtx1 = topTable(fit1,coef=2,adjust='BH',n=Inf) #coef要么必须等于2， 要么是个字符串；关于adjust的设置，说明书中13.1章有描述，BH是最流行的设置
# topTable 默认显示前10个基因的统计数据；使用选项n可以设置,n=Inf就是不设上限，全部输出
DEG_mtx1 = na.omit(mtx1) #去除缺失值
head(DEG_mtx1)

解释一下代码：design1=model.matrix(~factor(group_list))这是一个设计矩阵，关于设计矩阵的一些知识，可以看这篇笔记《线性模型》。

注：这里需要补充一下线性代数的相关知识。

lmFit函数是用于构建一个线性模型；它的参数一个是表达矩阵，一个是分组对象

eBays函数是用于构建一个微阵列线性模型。

20.火山图

设定绘图数据的阈值，如下所示：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
DEG=DEG_mtx1
logFC_cutoff <- with(DEG,mean(abs(logFC)) + 2*sd(abs(logFC)) ) 
DEG$result = as.factor(ifelse(DEG$P.Value < 0.05 & abs(DEG$logFC) > logFC_cutoff,
                              ifelse(DEG$logFC > logFC_cutoff ,'UP','DOWN'),'NOT'))
this_tile <- paste0('Cutoff for logFC is ',round(logFC_cutoff,3), #round保留小数位数
                    '\nThe number of up gene is ',nrow(DEG[DEG$result =='UP',]) ,
                    '\nThe number of down gene is ',nrow(DEG[DEG$result =='DOWN',])
)
ggplot(data=DEG, aes(x=logFC, y=-log10(P.Value), color=result)) +
  geom_point(alpha=0.4, size=1.75) +
  theme_set(theme_set(theme_bw(base_size=20)))+
  xlab("log2 fold change") + ylab("-log10 p-value") +
  ggtitle( this_tile ) + theme(plot.title = element_text(size=15,hjust = 0.5))+
  scale_colour_manual(values = c('blue','black','red'))

mark

结论

与R语言有关的这20道题中涉及的知识不少，代码背后的数理统计与线性代数知识目前比较欠缺，需要补充。

参考资料

  1. ExpressionSet 对象简单讲解
  2. 生信人的20个R语言习题
  3. 用limma对芯片数据做差异分析
  4. basic visualization for expression matrix
  5. R语言20练习题【完整版】
  6. 绝对中位差Median Absolute Deviation