分享是一种态度
单细胞RNA-seq分析介绍 单细胞RNA-seq的设计和方法 从原始数据到计数矩阵
R_refresher
项目reviewing_R.R
data
和figures
的文件夹data
文件夹(https://github.com/hbctraining/DGE_workshop_salmon/blob/master/data/raw_counts_mouseKO.csv?raw=true)tidyverse
库read.csv()
读取所下载的文件并保存为counts
object/variable1library(tidyverse)
2counts <- read.csv("data/raw_counts_mouseKO.txt")
3class(counts)
4str(counts)
我们正在对p53野生型(WT)和敲除(KO)基因型的癌症样本进行RNA-seq。有8个样本,每个样本4个重复。编写R代码构建,如下所述。
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、Fmeta
对象的tibble 并将其命名为meta_tb
(确保不会丢失行名!) 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/