生信技能树R作业答案-高级

为达到更好的浏览效果,可转至RPubs进行查看

1 安装一些R包

数据包: ALL, CLL, pasilla, airway 软件包:limma,DESeq2,clusterProfiler 工具包:reshape2 绘图包:ggplot2

1
2
3
4
5
6
7
8
9
10
rm(list = ls())
if(F){
source("http://bioconductor.org/biocLite.R")
options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
options(BioC_mirror="http://mirrors.ustc.edu.cn/bioc/")#修改镜像,安装会加速
BiocManager::install("clusterProfiler")
BiocManager::install("ComplexHeatmap")
BiocManager::install("maftools")
BiocManager::install("reshape2")
}

2了解ExpressionSet对象.

比如CLL包里面就有data(sCLLex) 找到它包含的元素,提取其表达矩阵(使用exprs函数),查看其大小

A.参考:http://www.bio-info-trainee.com/bioconductor_China/software/limma.html B.参考:https://github.com/bioconductor-china/basic/blob/master/ExpressionSet.md

1
2
3
4
library(CLL)
data("sCLLex")
View(sCLLex)#查看其包含的元素
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
1
2
exp <- exprs(sCLLex)#提取表达矩阵
dim(exp)#查看其大小
## [1] 12625    22
1
2
group <- pData(sCLLex)#提取分组信息
table(group)#设计矩阵
##         Disease
## SampleID progres. stable
##    CLL11        1      0
##    CLL12        0      1
##    CLL13        1      0
##    CLL14        1      0
##    CLL15        1      0
##    CLL16        1      0
##    CLL17        0      1
##    CLL18        0      1
##    CLL19        1      0
##    CLL2         0      1
##    CLL20        0      1
##    CLL21        1      0
##    CLL22        0      1
##    CLL23        1      0
##    CLL24        0      1
##    CLL3         1      0
##    CLL4         1      0
##    CLL5         1      0
##    CLL6         1      0
##    CLL7         1      0
##    CLL8         1      0
##    CLL9         0      1

3了解 str,head,help函数

作用于 第二步提取到的表达矩阵

1
str(exp)
##  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" ...
1
head(exp)
##           CLL11.CEL CLL12.CEL CLL13.CEL CLL14.CEL CLL15.CEL CLL16.CEL
## 1000_at    5.743132  6.219412  5.523328  5.340477  5.229904  4.920686
## 1001_at    2.285143  2.291229  2.287986  2.295313  2.662170  2.278040
## 1002_f_at  3.309294  3.318466  3.354423  3.327130  3.365113  3.568353
## 1003_s_at  1.085264  1.117288  1.084010  1.103217  1.074243  1.073097
## 1004_at    7.544884  7.671801  7.474025  7.152482  6.902932  7.368660
## 1005_at    5.083793  7.610593  7.631311  6.518594  5.059087  4.855161
##           CLL17.CEL CLL18.CEL CLL19.CEL CLL20.CEL CLL21.CEL CLL22.CEL
## 1000_at    5.325348  4.826131  5.212387  5.285830  5.581859  6.251678
## 1001_at    2.350796  2.325163  2.432635  2.256547  2.348389  2.263849
## 1002_f_at  3.502440  3.394410  3.617099  3.279726  3.391734  3.306811
## 1003_s_at  1.091264  1.076470  1.132558  1.096870  1.138386  1.061184
## 1004_at    6.456285  6.824862  7.304803  8.757756  6.695295  7.372185
## 1005_at    5.176975  4.874563  4.097580  9.250585  7.657463  7.683677
##           CLL23.CEL CLL24.CEL CLL2.CEL CLL3.CEL CLL4.CEL CLL5.CEL CLL6.CEL
## 1000_at    5.480752  5.216033 5.966942 5.397508 5.281720 5.414718 5.460626
## 1001_at    2.264434  2.344079 2.350073 2.406846 2.341961 2.372928 2.356978
## 1002_f_at  3.341444  3.798335 3.427736 3.453564 3.412944 3.411922 3.396466
## 1003_s_at  1.046188  1.082169 1.501367 1.191339 1.139510 1.153548 1.135671
## 1004_at    7.616749  6.839187 7.545673 8.802729 7.307752 7.491507 8.063202
## 1005_at    8.700667  5.546996 9.611601 5.732182 6.483191 6.072558 9.919125
##           CLL7.CEL CLL8.CEL CLL9.CEL
## 1000_at   5.897821 5.253883 5.214155
## 1001_at   2.222276 2.254772 2.358544
## 1002_f_at 3.247276 3.255148 3.365746
## 1003_s_at 1.074631 1.166247 1.151184
## 1004_at   7.014543 8.019108 7.432331
## 1005_at   5.149411 4.989604 5.339848
1
help(exp)
## starting httpd help server ... done

4安装并了解 hgu95av2.db 包.

看看 ls(“package:hgu95av2.db”) 后 显示的那些变量

1
2
library(hgu95av2.db)
ls("package:hgu95av2.db")
##  [1] "hgu95av2"              "hgu95av2.db"          
##  [3] "hgu95av2_dbconn"       "hgu95av2_dbfile"      
##  [5] "hgu95av2_dbInfo"       "hgu95av2_dbschema"    
##  [7] "hgu95av2ACCNUM"        "hgu95av2ALIAS2PROBE"  
##  [9] "hgu95av2CHR"           "hgu95av2CHRLENGTHS"   
## [11] "hgu95av2CHRLOC"        "hgu95av2CHRLOCEND"    
## [13] "hgu95av2ENSEMBL"       "hgu95av2ENSEMBL2PROBE"
## [15] "hgu95av2ENTREZID"      "hgu95av2ENZYME"       
## [17] "hgu95av2ENZYME2PROBE"  "hgu95av2GENENAME"     
## [19] "hgu95av2GO"            "hgu95av2GO2ALLPROBES" 
## [21] "hgu95av2GO2PROBE"      "hgu95av2MAP"          
## [23] "hgu95av2MAPCOUNTS"     "hgu95av2OMIM"         
## [25] "hgu95av2ORGANISM"      "hgu95av2ORGPKG"       
## [27] "hgu95av2PATH"          "hgu95av2PATH2PROBE"   
## [29] "hgu95av2PFAM"          "hgu95av2PMID"         
## [31] "hgu95av2PMID2PROBE"    "hgu95av2PROSITE"      
## [33] "hgu95av2REFSEQ"        "hgu95av2SYMBOL"       
## [35] "hgu95av2UNIGENE"       "hgu95av2UNIPROT"
1
ls(hgu95av2.db)
## [1] "conn"        "packageName"

5利用注释包找到某基因对应的探针

理解 head(toTable(hgu95av2SYMBOL)) 的用法,找到 TP53 基因对应的探针ID

1
head(toTable(hgu95av2SYMBOL))
##    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
1
ids <- toTable(hgu95av2SYMBOL)

5.1 subset函数

1
a1 <- subset(ids, ids$symbol=='TP53')

5.2 grep函数

1
a4 <- ids[grep("^TP53$", ids$symbol),]

5.3 filter函数(依赖dplyr包)

1
a5 <- filter(ids, ids$symbol=='TP53')

6理解基于与探针的对应关系

总共多少个基因,基因最多对应多少个探针, 是哪些基因,是不是因为这些基因很长,所以在其上面设计多个探针呢?总共多少gene

1
length(unique(ids$symbol))
## [1] 8585

对应探针多的是那些gene

1
2
ids_fre <- data.frame(table(ids$symbol))
tail(ids_fre[order(ids_fre$Freq),])
##        Var1 Freq
## 8365 YME1L1    7
## 2809  GAPDH    8
## 3638 INPP4A    8
## 4720    MYB    8
## 6010 PTGER3    8
## 7339  STAT1    8

或者用sort命令

1
tail(sort(table(ids$symbol)))
## 
## YME1L1  GAPDH INPP4A    MYB PTGER3  STAT1 
##      7      8      8      8      8      8

不管是Agilent芯片,还是Affymetrix芯片,上面设计的探针都非常短。 最长的如Agilent芯片上的探针,往往都是60bp,但是往往一个基因的长 度都好几Kb。因此一般多个探针对应一个基因,取最大表达值探针来作为 基因的表达量。 # 7.找到芯片有而hgu95av2.db中没有对应基因名的探针 第二步提取到的表达矩阵是12625个探针在22个样本的表达量矩阵, 找到那些不在 hgu95av2.db 包收录的对应着SYMBOL的探针。 提示:有1165个探针是没有对应基因名字的。 ## 第一种方法:match函数

7找到芯片有而注释包没有symbole的探针

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

1
dim(exp)
## [1] 12625    22
1
dim(hgu95av2SYMBOL)
## [1] 11460     2
1
2
b1 <- exp[-match(ids$probe_id, rownames(exp)),]
dim(b1)
## [1] 1165   22

第二种方法:%in%函数

1
2
b2 <- exp[!(rownames(exp) %in% ids$probe_id),]
b1==b2

8删除注释包没有symbol的探针

删除注释包中没有对应基因名的探针

过滤表达矩阵,删除那1165个没有对应基因名字的探针。

1
2
exp_filter_no_symbol <- as.data.frame(exp[match(ids$probe_id, rownames(exp)),])
exp_filter_no_symbol2 <- as.data.frame(exp[rownames(exp) %in% ids$probe_id,])

9多个探针保留一个

多个探针对应一个基因的情况下,只保留在所有样本里面平均表达量最大的那个探针。
提示,理解 tapply,by,aggregate,split 函数 , 首先对每个基因找到最大表达量的探针。 然后根据得到探针去过滤原始表达矩阵 这种方法能够一定程度上增加差异基因的数目,但也容易造成假阳性的结果 https://blog.csdn.net/tommyhechina/article/details/80468361 合并探针,先看有无NA值,若有,可以删除或填充 先检查有无,若返回值>0,说明有NA值。

1
length(which(is.na(exp_filter_no_symbol)))
## [1] 0

结果是0,说明没有NA值。如果有则用impute包来进行填充 exp_symbol<- impute.knn(exp_symbol)$data 下面开始进行合并,用aggregate命令

1
2
3
4
5
ids <- toTable(hgu95av2SYMBOL)
ids$median <- apply(exp_filter_no_symbol, 1, median)
ids <- ids[order(ids$symbol, ids$median, decreasing = T),]
ids_filtered <- ids[!duplicated(ids$symbol),]
dim(ids_filtered)
## [1] 8585    3

10过滤后的表达矩阵行名改为SYMBOL

把过滤后的表达矩阵更改行名为基因的symbol,因为这个时候探针和基因是一对一关系了。

1
2
3
4
exp_filter_no_symbol$probe_id <- rownames(exp_filter_no_symbol)
exp_ids <- merge(ids_filtered, exp_filter_no_symbol, by ='probe_id' )
rownames(exp_ids) <- exp_ids$symbol
exp_sym <- exp_ids[,-c(1:3)]

11画第一个和所有样本的基于表达量的图

对第10步得到的表达矩阵进行探索,先画第一个样本的所有基因的表达量的boxplot,hist,density , 然后画所有样本的 这些图 参考:http://bio-info-trainee.com/tmp/basic_visualization_for_expression_matrix.html 理解ggplot2的绘图语法,数据和图形元素的映射关系

1
boxplot(exp_sym$CLL11.CEL)

1
boxplot(exp_sym)

用ggplot2来画 ggplot2画图要改变数据形式

1
2
3
4
5
library(reshape2)
exp_g <- melt(exp_ids[,-c(1,3)], id.vars = 'symbol')
exp_g$group <- rep(group$Disease, each = nrow(exp_sym))
colnames(exp_g)[2] <- 'sample'
head(exp_g)
##    symbol    sample    value    group
## 1   MAPK3 CLL11.CEL 5.743132 progres.
## 2    TIE1 CLL11.CEL 2.285143 progres.
## 3 CYP2C19 CLL11.CEL 3.309294 progres.
## 4   CXCR5 CLL11.CEL 7.544884 progres.
## 5   DUSP1 CLL11.CEL 5.083793 progres.
## 6   MMP10 CLL11.CEL 3.252213 progres.

第一个样本的

1
library(ggplot2)
## Registered S3 methods overwritten by 'ggplot2':
##   method         from 
##   [.quosures     rlang
##   c.quosures     rlang
##   print.quosures rlang
1
2
sample_1 <- subset(exp_g, exp_g$sample=='CLL2.CEL')
head(sample_1)
##         symbol   sample    value  group
## 120191   MAPK3 CLL2.CEL 5.966942 stable
## 120192    TIE1 CLL2.CEL 2.350073 stable
## 120193 CYP2C19 CLL2.CEL 3.427736 stable
## 120194   CXCR5 CLL2.CEL 7.545673 stable
## 120195   DUSP1 CLL2.CEL 9.611601 stable
## 120196   MMP10 CLL2.CEL 3.380044 stable
1
2
p_s1 <- ggplot(sample_1, aes(x = sample, y = value, fill = 'lightblue'))
p_s1 + geom_boxplot()

1
p_s1 + geom_violin()+geom_boxplot(fill = 'white')

1
2
ggplot(sample_1, aes(x = value, fill = 'lightblue'))+
geom_histogram(bins = 500)

1
2
ggplot(sample_1, aes(x = value, fill = 'lightblue'))+
geom_density()

下面是所有样本的

1
2
p <- ggplot(exp_g, aes(x = sample, y = value, fill = group))
p + geom_boxplot()