logo
搜索
我的图书馆

发文章

发文工具

其他工具

留言交流
搜索
分享

【原】生信人应该这样学R语言系列视频学习心得笔记分享

 健明 2021-07-14
耗费半年的时间精心制作了成套的生物信息学入门视频教程，并且在生信技能树联盟平台发布了这个长达74个小时全套生物信息学入门视频：生信技能树视频课程学习路径，这么好的视频还免费！
在B站看了看，大家学的热火朝天， 接下来我们就一个个知识点进行专题介绍，主要是一些优秀学生的笔记分享，希望大家在学习的过程中也能吸收到我传达的学习经验，人生感悟，只要你发给我笔记（邮箱 jmzeng1314@163.com），就有惊喜！
专题历史目录：
3个学生的linux视频学习笔记
接下来介绍R语言：

【生信技能树】生信人应该这样学R语言


1.课前准备：周末生物信息学培训班准备工作【http://www./3727.html】

2.有好奇心，主动学习的能力

3.不可一下运行全部代码，避免不知道哪一步报错。一步步运行，慢+理解=进步

4.逻辑思考代码含义，比较，修改参数推理含义

5.敢于折腾

6.找规律

7.代码错误，看变量类型 class(), mode(), typeof()

8.R语言高手认识每个配置元素

9.通过project组织项目

10.报错的两种情况：自己的逻辑思维出问题(你以为的不是你以为的)；系统提示报错

11.函数初学看参数，进阶把常用参数记在心里

12.python,R,Linux本质上木有区别，取决于对哪种语法更熟悉，目的都是处理数据

1. 介绍R语言及Rstudio

  • 了解R，Rstudio及R包

    • 安装的包在packages中检查

      .libPaths()  #找安装路径

    • 帮助文档，帮忙看路径

      ?substring

    • 学会代码所蕴含的含义

      > tmp='abcde'
      > tmp
      [1] "abcde"
      > substring(tmp,1,3)
      [1] "abc"

    • 定位文件、设置文件位置

      getwd()
      setwd()  

    • History- to console, to source

    • 环境变量：赋值后会被记录

    • plot画图，如果没有出现所画的图形，要关掉之前开的窗口，再画图

      plot(1:10)
      png('tmp.png')
      dev.off()
      plot(1:10)

    • 下载R语言软件

    • 下载Rstudio编辑器

    • 安装一些必要的包，了解CRAN和bioconductor, 了解镜像并设置

2. R语言基础变量讲解

1.使用源代码或者重用R包，报错是因为没有基础知识，即变量类型

2.看五本以上的R书

  • 从无到有创建各种变量，掌握很多函数

    • 理解索引

    • 理解行列

    • 理解取子集的两种方式

    • grep()搜索函数

      index1 = grep('RNA-Seq', a$Assay_Type)
      index2 = grepl('RNA-Seq', a$Assay_Type)
      b = a[index1,]  # 下标
      b = a[index2,]  # 索引

    • list取元素用双括号来取，不然只能取到子list

    • 向量型(Vector)

      > a=c(1,2,3)
      > a
      [1] 1 2 3
      > class(a) # 数值类型不相容的时候查变量类型是否不符合函数类型
      [1] "numeric"

    • 双整型double

    • 整型integer

    • 字符型character

    • 逻辑型logical

    • 复数类型complex

    • 原始类型raw

    • 矩阵型(Matrix)和数组型(Array)

    • 数据框架型(Data Frame)及列表型(List)

    • 创建字符时，有字符有数值，创建的向量就全变成字符  #向量有优先级

      > a=c(1,'b',2)
      > class(a)
      [1] "character"

    • 字符串要加引号，单双引号没啥区别。特殊情况下有区别，单双引号都在

    • 二维和多维

    • 每个元素可以用$符号取

    • 概括来说，R可以识别六种基本类型的原子型向量，可以具有名称属性或者维度属性

    • 用函数或者内置变量创造变量

      > a=1:10
      > a
      [1]  1  2  3  4  5  6  7  8  9 10
      > a = seq(1:10)
      > a
      [1]  1  2  3  4  5  6  7  8  9 10
      > a = LETTERS
      > a
      [1] "A" "B" "C" "D" "E" "F" "G" "H" "I" "J" "K" "L" "M" "N"
      [15] "O" "P" "Q" "R" "S" "T" "U" "V" "W" "X" "Y" "Z"
      > a = LETTERS[1:10]
      > a
      [1] "A" "B" "C" "D" "E" "F" "G" "H" "I" "J"

    • 五种变量结构（class属性）

    • 变量的类型判断及转换

      is. X
      as. X

    • 数据变量的各种操作

    • a=data.frame(n=LETTERS[1:10],v=1:10)

    • rownames(a)=paste0('row',1:10)

  • 为什么画个热图这么难

3. 外部数据导入导出

  • Excel表格到R语言转换

    a = read.table('', head = T, sep = '\t')
    b = read.table('', comment.char = '!', head = T, sep = '\t')  # 带感叹号的不读取
    write.csv(b,'.csv')
    d = read.csv('.csv')

    read.table('', head = T, sep = '\t')
    save(b,file = 'b_input.Rdata')
    load(file = 'b_input.Rdata')
    b = b[,-1]  # 取第一行
    b = log2(b)
    pheatmap::pheatmap(b[1:10,])

  • 建议把excel转成csv来读取

4. 中级变量操作

变量类型会贯穿所有分析、报错之中（你以为的不是你以为的）

打开Project

进阶之后，一个函数通用多个功能

几乎所有英文单词都是一个函数，函数无穷无尽，找规律

  • 所有函数有参数，很多是互通的， eg sort, max, min, fivenum(四分位值)

  • 向量化操作

不断练习

几乎所有英文单词都是一个函数, 如果不是，可以把它变成一个函数

变量可以叫任何名字，关键要知道它指代是啥

函数可以满足各种需求，自己可以创造函数

  • 从无到有创建各种变量，掌握很多函数

    • 五种变量结构(class属性)

    • 变量类型判断及转换

    • 数据变量的各种操作

    • a.data.frame(n=LETTWRS[1:10], v = 1:10)

    • rownames(a)=paste0('row', 1:10)

  • 变量，数据类型和结构，函数运算，条件控制和循环

    • 习惯看str()，head()

    • 改变变量类型 as.numeric(b[1,])

    • 批量操作(多种方法可以做同一件事)

    • 应用函数

    • 循环的思想

      for(i in 1: nrow(n){                
         print(mean(as.numeric(b[i,])))
      })                                      # 类似于C语言
         apply(b,1，function(x){
             mean(x)
         })                                 #  矩阵b, 行叫1，列叫2
      apply(b,1,max)
      rowMax = function(x){
         apply(x,1,max)
      }

    • 有意义的统计量

      sort(apply(b,1,sd),decreasing = T)[1:50]
      cg = names(sort(apply(b,1,sd),decreasing = T)[1:50])
      cg
      sample(1:nrow(b)，50)    # 随机数
      pheatmap::pheatmap(b[sample(1:nrow(b)，50,])    # 随机数
      pheatmap::pheatmap(b[cg,])

    • abs, sqrt :戒对治，平方根

    • Log, log10, log2, exp: 对数与指数函数

    • Sin, cos, tan, acos, atan, atan2: 三角函数

    • sinh, cosh, tanh, asinha, acosh, aranh：双曲函数

    • 集合运算，reshape, merge总结

    • 思考一下excel表格里面有变量类型吗

      a = read.table('XXtable.txt', head = T
                   sep = '\t')
      b = read.table('GSE17215_series_matrix.txt.gz',
                   comment,char = '!', head = T,
                   sep = '\t')
      write.csv(b, 'GSE17215_series_matrix.csv')
      write.table(b,'tmp.csv', sep = ',')
      ##把行名去掉
      d = read.csv('GSE17215_series_matrix.csv')
      # readline 读入之后拆分

    • 读取一个txt文档看看

    • 简单运算（加减乘除）

    • 逻辑运算和布尔运算

    • 数学函数（一通百通）

    • 简单统计量（表达矩阵）

    • 进阶技能

5. 热图

  • 为什么画个热图这么难

  • 热图—不同大小的数据，代表不同颜色的深浅

    paste('a1', 1:20, sep = '_')

    • 跟着帮助文档的example学，边学边看

    • 先看最基础的

    • 看代码，慢慢理解

    • 模拟数据，内置数据

    • 整理好的真实数据

    • 原始数据提取热图需要的数据

    • 学R语言的动力

    • 画热图的包有十几个

6. 选取差异明显的基因的表达量矩阵绘制热图

rm(list = ls())   #魔幻操作，一键清空
library(pheatmap)
a1 = rnorm(100)
dim(a1) = c(5,20)
pheatmap(a1)
a2 = rnow(100)+2
dim(a2) = c(5,20)
library(pheatmap)
pheatmap(a1, cluster_rows = F, cluster_cols = F)
pheatmap(cbind(a1,a2))
pheatmap(cbind(a1,a2), show_rownames = F, show_colnames = F)

  • 拉平极差值

7. ID转换

  • 学会R语言操作的高级技巧

    • R语言批量处理

  • 找对应表格写代码转换

  • 一些函数

    strsplit('','[.]')  #根据点号分割
    duplicated()  #去重

  • 一些包

    org.Hs.eg.db    #在包里有基因注释关系


8. 任意基因任意癌症表达量分组的生存分析

  • 在线制作生存曲线的工具，oncolnc【http://www./】

  • 学习包   survival

    • 微调参数，大部分都不是自己写的

    • 学会R语言基础变量结构，基础函数操作，应用高级代码

    • 统计学、机器学习——数据做好，R包进，再进行理解


9. 任意基因任意癌症表达量和临床性状关联

  • 某个基因在某个癌症的表达量，关联临床信息

  • Cbioportal【http://www./】

  • 读取文本文件

  • 尽量避免用函数名来命名


10. 表达矩阵样本的相关性

  • 看两个变量的相关性

    > cor(1:10,1:10)
    [1] 1
    > a = rnorm(10)
    > b = rnorm(10)
    > cor(a,b)
    [1] -0.1555608
    > a = rnorm(10)
    > b = 10*a+rnorm(10)
    > cor(a,b)
    [1] 0.9971822

  • R包分类

    • 数据包（加载数据）

      view()
      dim() #看维度

    • 功能函数包

    • 注释包（各种芯片、基因间转换）

  • 知道代码背后的结构是啥，知道啥啥意思


11. 芯片表达矩阵下游分析

  • 表达矩阵/分组信息

  • DEG by Lima

    • 加入了回归相关的校验

    • 构造比较矩阵

    • 画火山图

    • 挑选变化比较大大基因

    • 富集分析

  • 做芯片高级分析时看包的说明书


12. RNA-seq表达矩阵差异分析

芯片的表达矩阵  exprSet = exprs(sCLLex)
RNA-seq表达矩阵 exprSet = assay(ariway)  # assay获取表达矩阵
# 对象时一锅粥都有，用@符号表示，用$符号来取，用大的包来分析
大部分包从基础变量开始，变量为矩阵、数据框、列表、数组、
exprSet[1:4,1:4]
boxplot(log(exprSet+1))  # normalization的一种方式

  • 两种写包的流派

    • 包装越好，越合适初学者用

    • 尽量不用任何高级的东西，纯粹自己写

  • RNA-Seq  NES分析

  • 分析方法

    • 法一：DESeq2

    • 多个探针对应同一个基因取最大表达量探针极简代码【http://www./3693.html】


13. R语言小作业

14. R学习资源介绍

学习资源（第二集）

  • 学习资源大全

  • R语言基础博客【https://www.jianshu.com/nb/22007361】

  • 统计学基础【https://mp.weixin.qq.com/s/OtB2h6f00U2SRZLzveJKfQ?token=1714554631&lang=zh_CN】

  • 可视化基础【http://www./english/】

  • 数据处理基础

  • bioconductor基础【https://bioconductor./BiocWorkshops/】

  • 通读一本书：R语言实战，R语言编程艺术

  • 初学者之友

    • 多打view(), str(), ?

    • 使用Rstudio编辑器提供的project形式组织隔离各个项目

    • cheatsheet大全

    • R语言初学笔记吾日三省吾身

    • R语言练习题大全

    • R黑魔法合集

    • 代码调试技巧

R语言为解释器

Rstudio为编辑器

  • 保留代码记录，提高写代码效率

    list.files()  #查看文件


15. R-2与excel的区别

打开R和R studio的正确方式：建文件夹—R project-打开R studio和R语言

Excel通过鼠标操纵数据，R语言通过编程操纵数据

R语言在很多情况下没有循环这个概念了

  • 定位getwd()

    a = read.table('   .txt', header = T, sep = '\t')

  • 排序

  • 求和

  • a[行号，列号]

    • b = a[1:4, 1:4]

    • seq(1, 10)

    • (d = b[.1])   #加括号是赋值的同时，打印出来

    • b[1, ]

    • b[, 1]

    • 取行, 取列

    • 固定表头

  • 变量，数据类型和结构，函数运算，条件控制和循环

    • 正则表达式

    • 合并与分列

    • 匹配与替换

    • 缺失值插补

    • 去重与排序

    • 控制流、循环与判断

    • 不停地取不同的元素做同样的事情

    • sum, mean, var, sd, min, max, range, median, IQR(四分位间距)等为统计量

    • Sort, order, rank与排序有关

    • 其他还有ave, fivenum, mad, quantity, stem等

    • 拿到数字自己算

    • 自己从头到尾看R语言实战

    • 向量型（vector)

    • 矩阵型（matrix）  #线性代数中用，一般不用

    • 数组型（array)

    • 列表型（list)

    • 数据框架型（data frame)：多行多列

      > c(1,2,3,4)==4    #c()代表构造一个向量
      [1] FALSE FALSE FALSE  TRUE
      > temp=c(1,2,3,4)
      > temp
      [1] 1 2 3 4
      > temp==4
      [1] FALSE FALSE FALSE  TRUE

      Excel可做的，R也可以做

      table(a[,18])   等同于   table(a$tissue)

    • 逻辑型，整型，数值型，复数型，字符型

    • table(a[,1]=='WXS')    #计数

    • which(a[,1]=='WXS')   #行号

    • 读取一个txt文档看看

    • 五种变量类型（class(), view())

    • 简单运算（加减乘除），集合运算，reshape以及merge总结

    • 数学函数

    • 简单统计量

    • 进阶技能

  • apply系列函数，因子的高阶用法，字符串的操作

    • 内置函数

    • Hadley Wickham开发的一个灵魂的字符串处理包stringr

      str_split(a$Library_Name,',',simplify = T)[,2]
      do.call(cbind, strsplit(a$Library_Name,','))[,2]
      # 殊路同归，多练习，同一问题有多种解决方案
      # troubleshooting能力

    • 字符串分割语数 strsplit()

    • 字符串连接语数 paste()

    • 计算字符串长度 nchar()

    • 字符串截取语数 substr()及substring()

    • 字符串替换函数 chartr()

    • 大小写转换语数 toupper(),tolower(), casefold()

    • 看一个包的help语法

    • vignette('stringr')

    • 创建因子

    • 更改level, 差异分析case和control的问题，绘图顺序，chr1-22, X, Y, MT等

    • cut语数，table函数

    • 函数名 含义

    • apply  apply

    • Apply list apply

    • sapply  simplify list apply

    • tapply  table apply

    • mapply. mutiple  list apply

    • rapply  recursively  apply

    • vapply. vector. apply

    • apply. Environment  apply

    • 可以省略掉所有循环

    • apply函数族的向量化运算是基于R语言函数来实现的，八大函数，参考张丹老师博客【http://blog./r-apply/】

    • 大部分新手错误都来自因子

    • 字符串的操作

  • 基础绘图及ggplot绘图语法    basic plot in R for biotrainee

    绘图不会的用goole搜代码，看测试数据用到啥数据，把自己的数据想办法弄成范例的样子，不用自己设计新的图形

    • 一张统计图就是从数据到几何对象（点、线、条形等）的图形属性（颜色、形状、大小等）的一个映射

    • 数据（data)，最基础的是可视化的数据和一系列图形映射（aesthetic mappings), 该映射描述了数据中的变量如何映射到可见到图形属性

    • 几何对象（geometric objects, geoms) 代表在图中实际看到的点，现，多边形等

    • 统计转换（statistical trassformations, stats)是对数据进行某种汇总，例如将数据分组创建直方图，或将一个二维的关系用线性模型进行解释

    • 标度（scales), 是将数据的取值映射到图形空间，例如颜色、大小或形状来表示不同的取值，展现标度的常见做法是绘制图例和坐标轴

    • 坐标系（coordination system, coord)描述数据是如何映射到图形所在的平面，同时提供看图所需的坐标轴和网格线

    • 分图（faceting）如何将数据分解为子集，以及如何对子集做图并展示

    • 主题（theme)控制细节显示，例如字体大小和图形的背景色

    • 基础语法看书即可【http:///jmzeng/markdown/basic-plot-R.html】

    • ggplot练习【http:///jmzeng/markdown/ggplot-in-R.html】

    • boxplot

    • 加title

  • 路还很长

    • 层次聚类会先计算样本之间的距离，每次将距离最近的点合并到同一个类，然后再计算类与类之间的距离，将距离最近的类合并为一个大类，不停合并，直到合并成了一个类

    • 动态聚类则是先抽几个点，把周围的点聚集起来，然后算每个类的重心或平均值什么的，以算出来的结果为分类点，不断重复，直到分类的结果收敛为止。

    • r语言使用dist(x, method = 'euclidean', drag = FALSE, upper = FALSE, p = 2)来计算距离。其中X是样本矩阵或者数据框，method表示计算哪种距离

    • r语言中使用hclust(d, method = 'complete', members = NULL)来进行层次聚类，其中d为距离矩阵，method为类的合并方法

    • r语言中主要使用kmeans, centers, tierce_max = 10, nastart = 1, algorithm=c("Hartigsn-WONG", "Lloyd", "Forgy", "MacQueen")来进行聚类

    • 百分号（%）符号，鬼画符系列

    • 你以为使用的是一个函数

  • 墙内群众特别配置

  • 作业【http://www./3415.html】

    20题做完，则R语言过关了


16. R-3简单统计及数学函数

  • Excel文件，另存为csv格式

    a = read.csv('    .csv')

    所有英文单词，在R中都是一个函数

    view(a)
    mean(a$MBases)
    t.test(rna$MBases, wxs$MBases)
    fivenum(rna$MBases)


17. R-4基础语法

差异表达分析，其他包本质上与limma包没有区别

所有的语言是殊路同归的，把linux学到极致，R语言可以不用；R语言学到极致，可以不用python

grep('WXS', a[,1]).  #返回行号

grepl('WXS', a[,1])  #返回true  or false

wxs = a[a[,1]=='WXS',]
wxs = a[grepl('WXS',a[,1]),]
wxs = a[which(a[,1])=='WXS',]
wxs = a[grep('WXS',a[,1]),]
# 以上四个殊路同归

  • 背下来lapply的结构

lapply(strsplit(tmp,','),function(x){cat(x)})
lapply(l[1:4],function(x){cat(x)})
for(x in l(1:4)){
   cat(x[2])
}

  • 管道的概念

  • 匹配与替换

18. 高级数据处理技巧

  • 使用R包


19. R-6-绘图该如何学

  • 善用table自动补全功能

  • 掌握R语言基本思维，可以猜出怎么用

    • ggpuber- 针对publication的R

    • STHDA【http://www./english/】

    • 没事多逛一下绘图的R网站，了解一下代码，跑一遍代码

    • R语言绘图设计哲学


20. R-7-作业

  1. 矩阵的数据结构中只能包含相同的类型或模式（数值型、字符型或逻辑型）的组成

  2. 在R中，对数据框进行操作时，若需要创建在with()结构以外存在的对象（保存在全局环境中），应该使用<<-

  3. NA-未定义，NaN-not a number，inf-正无穷，-inf-负无穷，NULL-缺失值

  4. 实操

diabetes = rep(c('t1','t2'),2)

  1. 查看长度、维度、类型、模式、行名、列名

    length, dim, str, model, colname, rowname

  2. head(), tail()

  3. 删除rm()

安装python用豆瓣

dog用阿里云

R 包清华， bioconductor用中科大？

  1. 要完完全全理解变量，注意细节问题

####21. windows电脑安装R及R studio

  • R【https://www./】

  • Rstidio【https://www./】


22. 史上最贴心R包安装示范

不是所有的R包都需要安装成功的

  1. 内陆的小伙伴修改镜像

  2. 所有的软件都安装在C盘

  3. 系统用户名最好不要用中文，写代码最怕中文字符串

  4. 不管是什么电脑，都务必安装好R及R studio

  5. 初次安装R包会有很多依赖包一起被安装

  • 安装下面的包，用于临床三线表

    • tableone

    • rJava

    • ReporteRs

  • 从bioconductor安装包

    • 不同的R版本需要对应不同的biconductor版本


23. 系统性的入门R语言

  • why学R语言

    • 各种统计分析（生存分析，差异分析，lasso回归）

    • 想画一个热图

    • 想做GEO数据分析

  • 各种搜索渠道

    • R语言实战

    • 百度云资料

    • 腾讯课堂

    • 网易课堂

    • 入门书籍

  • 了解并安装R及Rstudio

    • 安装bioconductor包

  • 了解R语言与excel异同

  • 创建数据框、列表、数组

  • 变量的基本操作

    • 凡是英文单词都是英文函数

    • 变量有向量因子、数据框、列表、数组

    • 数据的高级操作

  • 高级分支

    科比退役的第。。。天  想他

    会创建R包，就代表R学好了

    • 统计学

    • 可视化

    • bioconductor与生物信息学

    • R语言操作数据库

    • shiny写网页

    • R语言操作图片

    • R语言写爬虫


24. R语言可视化基础

  • basic visualization for expression matrix【http:///tmp/basic_visualization_for_expression_matrix.html】

    • 用R markdown HTML制作网页

    • 查看R markdown cheatsheet

  • 画图学基本代码。把画的内容记录好

    • 慢慢熟悉，跟着做然后理解，数据对应成图形

    • 用测试数据练习

  • 如何通过Google来使用ggplot2可视化(上)【http://www./2332.html】

  • ggplotcheatsheet


25. R语言可视化进阶

  • google画图，会搜索

  • 推荐文献代码


26. 用shiny写一个计算器

编程很容易


27. R语言操作图片

  • 用什么学什么，现搜


28. R语言写爬虫

  • 用R语言的RCurl包结合XML包批量下载生信课件【http://www./799.html】

  • 解析规律

  • 爬虫高阶（自己探索）


29. 用R语言爬去生信软件列表到思维导图

  • 单细胞转录组数据分析软件大全

  • 了解网页基础知识

    • 背后的dom结构

    • 期刊探索【https://cloud.tencent.com/developer/article/1050838】

    • peerJ期刊探索【https://cloud.tencent.com/developer/article/1050850】

    • H2, H3,

    • HTML 5 教程【http://www.w3school.com.cn/html5/index.asp】

    • 了解你要爬虫的目标网站

  • 幕布软件学习


30. 用R语言进行ID转换

  • ID转换大全【https://cloud.tencent.com/developer/article/1054744】

  • 生物信息各类ID转换攻略【https://www.jianshu.com/p/5bef111ae3ba】

  • R语言中的排序，集合运算，reshape，以及merge总结【https://mp.weixin.qq.com/s/0u1Ms5cy3KonoC14w9GWgg】

    在R语言中做数据转换非常简单，找到想要的注释关系在哪里

    自己学基础语法，了解是啥，自己查如何用

  • 对entrez ID 跟symbol表示的基因进行KEGG注释【http://www./thread-409-1-1.html】

  • 了解生信技能树    ID+转换    公众号内容

  • 生信编程实战5个月传送门【http://www./thread-1075-1-1.html】


31. 多个差异分析结果直接量取交集并集

自己看编程的基本语法

For, while

  • 多个差异分析结果直接量量取交集并集

  • R，Per, Python处理相同任务进行比较

上面是优秀学员惠美的R视频教程学习笔记！

~~~~~~假装这里是分割线~~~~~~~

下面是小洁的 R语言小练习题的答案：

准备工作--安装所需的包

cran_packages <- c('tidyverse',
                   'ggpubr',
                   'ggstatsplot')
Biocductor_packages <- c('org.Hs.eg.db',
                         'hgu133a.db',
                         'CLL',
                         'hgu95av2.db',
                         'survminer',
                         'survival',
                         'hugene10sttranscriptcluster',
                         'limma')
if(length(getOption("CRAN"))==0) options(CRAN="https://mirrors.tuna./CRAN/")
for (pkg in cran_packages){
  if (! require(pkg,character.only=T) ) {
    install.packages(pkg,ask = F,update = F)
    require(pkg,character.only=T) 
  }
}
# first prepare BioManager on CRAN
if(length(getOption("CRAN"))==0) options(CRAN="https://mirrors.tuna./CRAN/")
if(!require("BiocManager")) install.packages("BiocManager",update = F,ask = F)
if(length(getOption("BioC_mirror"))==0) options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
# use BiocManager to install
for (pkg in Biocductor_packages){
  if (! require(pkg,character.only=T) ) {
    BiocManager::install(pkg,ask = F,update = F)
    require(pkg,character.only=T) 
  }
}

1.找到指定ensembl ID与symbol的对应关系

ENSG00000000003.13
ENSG00000000005.5
ENSG00000000419.11
ENSG00000000457.12
ENSG00000000460.15
ENSG00000000938.11

思路：在注释包中有gene_id与symbol、gene_id与ensembl_id的对应关系。

#将以上基因id保存在a.txt，存放于工作目录下。
rm(list=ls())
options(stringsAsFactors = F)

a=read.table('e1.txt')

g2s <- toTable(org.Hs.egSYMBOL)
g2e <- toTable(org.Hs.egENSEMBL)
a$V1 = apply(a[1], 1,function(x){
  str_split(x,'[.]')[[1]][1]
}) %>% unlist
colnames(a) <- 'ensembl_id'
tmp <- merge(a,g2e, by="ensembl_id")
result <- merge(tmp,g2s, by="gene_id")[-1]

2.根据探针名找对应symbol ID

1053_at
117_at
121_at
1255_g_at
1316_at
1320_at
1405_i_at
1431_at
1438_at
1487_at
1494_f_at
1598_g_at
160020_at
1729_at
177_at

思路：找到注释包中探针与symbol的对应关系然后取子集

rm(list=ls())
options(stringsAsFactors = F)
library(hgu133a.db)
p2s=toTable(hgu133aSYMBOL)
a=read.table('e2.txt')
colnames(a) <- colnames(p2s)[1]
# 方法一：利用merge
tmp1 <- merge(a,p2s, by="probe_id")
# 方法二：利用match得到第一组向量在第二组中的坐标
tmp2 <- p2s[match(a$probe_id,p2s$probe_id),]
## 附：判断得到的两组结果是否一致
# 法一：
identical(tmp1,tmp2) #返回逻辑值
# 法二：
dplyr::setdiff(tmp1,tmp2) #返回两组的差别【没差就返回空】

3.根据symbol找基因表达量并作图

找到R包CLL内置的数据集的表达矩阵里面的TP53基因的表达量，并且绘制在 progres.-stable分组的boxplot图,并通过 ggpubr 进行美化。
探针的三大内容：表达矩阵assay/exprs、探针信息featureData、样本信息phenoData

# 从内置数据集的表达矩阵中找TP53基因的表达量
rm(list=ls())
options(stringsAsFactors = F)
suppressMessages(library(CLL))
data(sCLLex)
# sCLLex
exprSet <- exprs(sCLLex) #探针的表达量
pdata <- pData(sCLLex) #sampleID与disease的对应关系
p2s <- toTable(hgu95av2SYMBOL) #探针与symbol的对应关系
p3 <- filter(p2s,symbol=='TP53')
# boxplot [find TP53 has 3 probe IDs]
probe_tp53 <- c("1939_at","1974_s_at","31618_at")
i = 3 #可换1，2
boxplot(exprSet[probe_tp53[i],] ~ pdata$Disease)

#用ggpubr作图
#http://www./english/articles/24-ggpubr-publication-ready-plots/
exp_tab <- rownames_to_column(as.data.frame(exprSet))
exp_tab2 <- gather(exp_tab,
                   key = 'sample',
                   value = 'exp',-1)
pdata <- rownames_to_column(pdata)
exp_tab3 <- merge(exp_tab2,pdata,by.x='sample',by.y='rowname')
i=1 ###可换1，2
dev.off()
p <- ggboxplot(filter(exp_tab3,rowname==probe_tp53[i]), 
               x = 'Disease',
               y = 'exp',
               color = "Disease", palette =c("#00AFBB", "#E7B800", "#FC4E07"),
               add = "jitter", shape = "Disease")
p

4.BRCA1基因表达量

找到BRCA1基因在TCGA数据库的乳腺癌数据集(Breast Invasive Carcinoma (TCGA, PanCancer Atlas))的表达情况

提示：使用http://www./index.do 定位数据集：http://www./datasets
该基因有四个亚型，用ggbetweenstats作图比较一下。

rm(list=ls())
options(stringsAsFactors = F)
#ID，四个亚型，表达量
f <- read.csv("e4-plot.txt", sep = "\t")
## boxplot
colnames(f) <- c("id", "subtype", "expression", "mut")
da <- f
library(ggstatsplot)
ggbetweenstats(data = da, 
               x = subtype,
               y = expression)
library(ggplot2)
ggsave("e4-BRCA1-TCGA.png")

5 .TP53 生存分析

找到TP53基因在TCGA数据库的乳腺癌数据集的表达量分组看其是否影响生存

提示使用：http://www./
思路：生存分析，TP53表达量分为高低两组做图比较

# Use http://www./ to get raw csv da

rm(list=ls())
options(stringsAsFactors = F)

f <- read.csv('e5-BRCA_7157_50_50.csv')
library(ggstatsplot)
ggbetweenstats(data = da,
               x = Group,
               y = Expression)
da <- f

library(ggplot2)
library(survival)
library(survminer)
table(da$Status)

da$Status <- ifelse(da$Status == "Dead", 1, 0)
survf <- survfit(Surv(Days,Status)~Group, data=da)
ggsurvplot(survf, conf.int = F, pval = T)

# complex survplot
ggsurvplot(survf,palette = c("orange", "navyblue"),
           risk.table = T, pval = T,
           conf.int = T, xlab = "Time of days",
           ggtheme = theme_light(),
           ncensor.plot = T)
ggsave("survival_TP53_in_BRCA_taga.png")

6.从表达矩阵中提取指定基因画热图

下载数据集GSE17215的表达矩阵并且提取下面的基因画热图

ACTR3B ANLN BAG1 BCL2 BIRC5 BLVRA CCNB1 CCNE1 CDC20 CDC6 CDCA1 CDH3 CENPF CEP55 CXXC5 EGFR ERBB2 ESR1 EXO1 FGFR4 FOXA1 FOXC1 GPR160 GRB7 KIF2C KNTC2 KRT14 KRT17 KRT5 MAPT MDM2 MELK MIA MKI67 MLPH MMP11 MYBL2 MYC NAT1 ORC6L PGR PHGDH PTTG1 RRM2 SFRP1 SLC39A6 TMEM45B TYMS UBE2C UBE2T

提示：根据基因名拿到探针ID，缩小表达矩阵绘制热图，没有检查到的基因直接忽略即可。

## Exercise 5: Retrive genes from GEO to plot heatmap

rm(list=ls())
options(stringsAsFactors = F)
#下载和表达矩阵
library(GEOquery)
GSE <- "GSE17215"
if(!file.exists(GSE)){
    geo <- getGEO(GSE, destdir = '.', getGPL = F, AnnotGPL = F)
    save(geo, file = paste0(GSE,'.eSet.Rdata'))
}
load(paste0(GSE,'.eSet.Rdata'))
expr <- exprs(geo[[1]])
p2s=toTable(hgu133aSYMBOL);head(p2s)
expr <- expr[p2s$probe_id,] #有的id找不到注释直接删掉
gp <- "ACTR3B ANLN BAG1 BCL2 BIRC5 BLVRA CCNB1 CCNE1 CDC20 CDC6 CDCA1 CDH3 CENPF CEP55 CXXC5 EGFR ERBB2 ESR1 EXO1 FGFR4 FOXA1 FOXC1 GPR160 GRB7 KIF2C KNTC2 KRT14 KRT17 KRT5 MAPT MDM2 MELK MIA MKI67 MLPH MMP11 MYBL2 MYC NAT1 ORC6L PGR PHGDH PTTG1 RRM2 SFRP1 SLC39A6 TMEM45B TYMS UBE2C UBE2T"
gp <- strsplit(gp, ' ')[[1]]
tmp <- dplyr::filter(p2s, p2s$symbol %in% gp)
tmp2 <- tibble::rownames_to_column(data.frame(expr),"probe_id")
tmp3 <- merge(tmp,tmp2,by="probe_id")

tmp3$median <- apply(tmp3[,3:ncol(tmp3)], 1, median)
tmp3 <- tmp3[order(tmp3$symbol, tmp3$median, decreasing = T),]
tmp3 <- tmp3[!duplicated(tmp3$symbol),]
rownames(tmp3) <- tmp3$symbol
tmp3 <- tmp3[,-c(1,2,ncol(tmp3))]

gp_expr <- log2(tmp3)
pheatmap::pheatmap(gp_expr, scale = "row") 

7.相关性计算和热图

下载数据集GSE24673的表达矩阵计算样本的相关性并且绘制热图，需要标记上样本分组信息
相关性分析：

rm(list=ls())
options(stringsAsFactors = F)

library(GEOquery)
GSE <- "GSE24673"
if(!file.exists(GSE)){
  geo <- getGEO(GSE, destdir = '.', getGPL = F, AnnotGPL = F)
  save(geo, file = paste0(GSE,'.eSet.Rdata'))
}
load(paste0(GSE,'.eSet.Rdata'))
expr <- exprs(geo[[1]])
dim(expr)
expr[1:4,1:4]
pdata <- pData(geo[[1]]) 
# 自己根据pdata第八列做一个分组信息矩阵
grp <- c('rbc','rbc','rbc',
         'rbn','rbn','rbn',
         'rbc','rbc','rbc',
         'normal','normal')

grp_df <- data.frame(group=grp)
rownames(grp_df) <- colnames(expr)
new_cor <- cor(expr)
pheatmap::pheatmap(new_cor, annotation_col = grp_df)

8.找到芯片平台对应的注释包

找到 GPL6244 platform of Affymetrix Human Gene 1.0 ST Array 对应的R的bioconductor注释包，并且安装它！
https://mp.weixin.qq.com/s/sVSsV2fWZOQmNd72Vb3SmQ
https://www.jianshu.com/p/40b560755cdf

options()$repos
options()$BioC_mirror 
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
options("repos" = c(CRAN="https://mirrors.tuna./CRAN/"))
BiocManager::install("hugene10sttranscriptcluster",ask = F,update = F)
options()$repos
options()$BioC_mirror

9.找到指定探针和对应的基因

rm(list=ls())
options(stringsAsFactors = F)

library(GEOquery)
GSE <- "GSE42872"
if(!file.exists(GSE)){
  geo <- getGEO(GSE, destdir = '.', getGPL = F, AnnotGPL = F)
  save(geo, file = paste0(GSE,'.eSet.Rdata'))
}
load(paste0(GSE,'.eSet.Rdata'))
expr <- exprs(geo[[1]])

dim(expr)
expr[1:4,1:4]

# 选出所有样本的(mean/sd/mad/)最大的探针

sort(apply(expr,1,mean),decreasing = T)[1]
sort(apply(expr,1,sd),decreasing = T)[1]
sort(apply(expr,1,mad),decreasing = T)[1]

下载数据集GSE42872的表达矩阵，并且分别挑选出所有样本的(平均表达量/sd/mad/)最大的探针，并且找到它们对应的基因。

10.limma 差异分析

下载数据集GSE42872的表达矩阵，并且根据分组使用limma做差异分析，得到差异结果矩阵

# 接第九题，得到表型信息，然后用limma做差异分析
pdata <- pData(geo[[1]]) 
grp <- unlist(lapply(pdata$title, function(x){
  strsplit(x, ' ')[[1]][4]
}))

suppressMessages(library(limma))
#limma needs：表达矩阵（expr：需要用log后的矩阵）、分组矩阵（design）、比较矩阵（contrast）
#先做一个分组矩阵～design，说明progres是哪几个样本，stable又是哪几个
design <- model.matrix(~0+factor(grp))
colnames(design) <- levels(factor(grp))
rownames(design) <- colnames(expr)
design
#再做一个比较矩阵【一般是case比control】
contrast<-makeContrasts(paste0(unique(grp),collapse = "-"),levels = design)
contrast

##step1
fit <- lmFit(expr,design)
##step2
fit2 <- contrasts.fit(fit, contrast) 
fit2 <- eBayes(fit2)  
##step3
mtx = topTable(fit2, coef=1, n=Inf)
DEG_mtx = na.omit(mtx) 
View(DEG_mtx)

# 火山图

DEG=DEG_mtx

if(T){
  logFC_cutoff <- with(DEG,mean(abs(logFC)) + 2*sd(abs(logFC)) )

  DEG$change = as.factor(ifelse(DEG$P.Value < 0.05 & abs(DEG$logFC) > logFC_cutoff,
                                ifelse(DEG$logFC > logFC_cutoff ,'UP','DOWN'),'NOT')
  )
  title <- paste0('log2FoldChange cutoff: ',round(logFC_cutoff,3),
                  '\nUp-regulated genes: ',nrow(DEG[DEG$change =='UP',]) ,
                  '\nDown-regulated genes: ',nrow(DEG[DEG$change =='DOWN',])
  )
}
library(ggplot2)
vol1 = ggplot(data=DEG, aes(x=logFC, y=-log10(P.Value), color=change)) +
  geom_point(alpha=0.4, size=1.75) +
  theme_set(theme_set(theme_bw(base_size=20)))+
  xlab("log2FoldChange") + ylab("-log10 p-value") +
  ggtitle( title ) + theme(plot.title = element_text(size=15,hjust = 0.5))+
  scale_colour_manual(values = c('blue','black','red')) # according to the levels(DEG$change)
print(vol1)

■   ■   ■

    转藏 分享 献花（0）

    来自： 健明 > 《待分类》

    举报/认领

    0条评论

    发表

    请遵守用户 评论公约

    类似文章 更多