基因表达差异分析前的准备工作

分享是一种态度

回顾

学习目标

• 了解R言语使用的各种数据类型和数据结构

• 在R中使用函数并了解如何获取有关参数的帮助

• 使用dplyr包中的管道(%>%)

• 了解ggplot2用于绘图的语法

配置

1. 创建一个新的项目目录

• 创建一个名为R_refresher项目

• 创建一个名为reviewing_R.R

• 项目目录中创建data和figures的文件夹

• 将counts文件下载到data文件夹(https://github.com/hbctraining/DGE_workshop_salmon/blob/master/data/raw_counts_mouseKO.csv?raw=true)

1. 加载库并读入数据，同时并思考以下问题

• 加载tidyverse 库

• 使用read.csv()读取所下载的文件并保存为counts object/variable

• 函数的语法是什么？

• 我们如何获得帮助 ？

• 什么是数据结构？

• R中有哪些主要的数据结构？

• 列的数据类型是什么？

• R中提供哪些数据类型？

1library(tidyverse)
2counts <- read.csv("data/raw_counts_mouseKO.txt")
3class(counts)
4str(counts)

创建向量/因子和数据框

我们正在对p53野生型(WT)和敲除(KO)基因型的癌症样本进行RNA-seq。有8个样本，每个样本4个重复。编写R代码构建，如下所述。

• 为每列创建vectors/factors(提示：您可以键入每个vectors/factors，如果您希望更快速的创建，可以尝试使用rep()函数)

• 将它们放到一个数据框中，这个数据框命名为meta

• 使用rowames()函数给数据框定义行名(提示：您可以键入行名作为向量，如果您希望该过程进行得更快，可以尝试使用paste0()函数)。
  创建好的数据框中应该包含sex、stage、genotype和myc ：
  

  
  

1sex <- rep(c("M","F"),4)
2stage <- rep(c(1,2),4)
3genotype <- c(rep("KO",4),rep("WT",4))
4myc <- c(23,4,45,90,34,35,9, 10)
5meta <- data.frame(sex=sex,
6                   stage=stage,
7                   genotype=genotype,
8                   myc=myc)
9rownames(meta) <- c(paste0(rep("KO",4),1:4),paste0(rep("WT",4),1:4))

探索数据

既然我们已经创建了元数据数据框，在执行任何分析之前获取一些关于数据的描述性统计数据通常是一个好习惯。

• 总结meta对象的内容，表示多少种数据类型？

• 检查meta数据框中的行名称是否与counts(内容和顺序)中的列名称相同

• 将现有 stage列转换为因子数据类型

1str(meta)
2all(rownames(meta) %in% colnames(counts))
3all(rownames(meta) == colnames(counts))
4
5meta$stage <- factor(meta$stage)
6str(meta)

提取数据

使用上一个问题中创建的meta数据框，执行以下练习(问题之间不是相互依赖)：

• 使用[]仅返回genotype和sex列

• 使用[]返回样本1、7和8的genotype值

• 用于filter()返回基因型为WT的样本的所有数据

• 使用filter()/ select()仅返回myc> 50的那些样本的stage和genotype列

• 在数据框的开头添加一个名为pre_treatment的列，其值为T、F、T、F、T、F、T、F

• 使用%>%创建meta对象的tibble 并将其命名为meta_tb（确保不会丢失行名！）

• 将列的名称更改为：“ A”，“ B”，“ C”，“ D”，“ E”

 1meta[,c(2,3)]
 2#or
 3meta[,c("stage","genotype")]
 4
 5meta[c(1,7,8),"genotype"]
 6
 7filter(meta, genotype == "WT")
 8
 9meta %>% filter(myc > 50) %>% select(stage, genotype)
10#or
11select(filter(meta, myc > 50), stage, genotype)
12
13pretreatment <- c(T, F, T, F, T, F, T, F)
14meta <- cbind(pretreatment, meta)
15###思考一下这样为什么会有问题
16
17meta_tb <- meta %>% rownames_to_column(var="sampleIDs") %>%as_tibble()
18
19colnames(meta_tb)[2:6] <- LETTERS[1:5]

可视化数据

通常，当我们使用各种图形进行可视化探索时，更容易看到数据的模式或性质。让我们使用ggplot2来探索基于基因型的Myc基因表达的差异。

• 使用theme_minimal()为KO和WT样本绘制Myc表达式的箱线图，并为绘图指定新的轴名和居中的标题。

1ggplot(meta) +
2  geom_boxplot(aes(x = genotype, y = myc)) +
3  theme_minimal() +
4  ggtitle("Myc expression") +
5  ylab("Myc level") +
6  xlab("Genotype") +
7  theme(plot.title = element_text(hjust=0.5, size = rel(2)))

为下游分析做准备

许多不同的统计工具或分析包都希望作为输入的所有数据都在列表结构中。让我们创建一个包含count和metadata的数据列表，为后续分析做准备。

• 使用meta和count对象创建名为project1的列表，并从两个数据框之一中提取所有样本名称创建一个新向量。

1project1 <- list(meta, counts, rownames(meta))
2project1

---

注：以上内容来自哈佛大学生物信息中心（HBC）的教学团队的生物信息学培训课程。原文链接：https://hbctraining.github.io/scRNA-seq/schedule/

往期回顾

---

---

---

如果你对单细胞转录组研究感兴趣，但又不知道如何入门，也许你可以关注一下下面的课程

• 生信爆款入门-全球听（买一得五）（第4期） 你的生物信息入门课

• 数据挖掘第2期（两天变三周，实力加量）医学生/医生首选技能提高课

• 生信技能树的2019年终总结 你的生物信息成长宝藏

看完记得顺手点个“在看”哦！

生物 | 单细胞 | 转录组丨资料每天都精彩

长按扫码可关注