jamesbang

V1

2022/10/05阅读:34主题:雁栖湖

🤩 scRNA-seq | 吐血整理的单细胞入门教程(数据格式和处理)(四)

1写在前面

我们还是在正式进行代码操作前想几个小问题:👇

  • 如何将单细胞数据导入R中?
  • 不同类型的数据/信息(如细胞信息、基因信息等)是如何存储和操作的?
  • 如何获得细胞基因的基本信息并对数据进行相应的过滤

2用到的包

目前常用的scRNA-seq分析包,包括SeuratScanpy(python)、ScaterMonocle2Monocle3等。🤒

rm(list = ls())
library(tidyverse)
library(SingleCellExperiment)
library(DropletUtils)

3数据格式

scRNA-seq的标准格式为SingleCellExperiment。 😘
Note! Seurat包有其自己的格式,即Seurat格式,可能因为Seurat太火了吧,越来越多的包都开始兼容Seurat格式的文件了。😂

我们拿到的数据通常是一个feature-by-sample的表达矩阵。
scRNA-seq分析中,我们一般需要从counts矩阵开始分析,代表每个cellfeaturereads/UMIfeature可以是geneisoformsexons。🤒

singlecellexperiment
singlecellexperiment

4SingleCellExperiment对象

4.1 读入数据

counts <- read.table("./molecules.txt", sep = "\t")
annotation <- read.table("./annotation.txt", sep = "\t", header = TRUE)
counts
counts
annotation
annotation

4.2 创建SingleCellExperiment对象

# 注意assays必须是matrix
tung <- SingleCellExperiment(
assays = list(counts = as.matrix(counts)),
colData = annotation
)

tung

Note! 这个SingleCellExperiment包含:

  • 19027个基因(行)和864个细胞(列);
  • 一个名为countsassays;
  • 预览部分基因名(rownames)和细胞名(colnames);
  • 基因metadata为空 (rowData);
  • 细胞metabdata(colData names)

510X 文件的读取

像上面那样创建SingleCellExperiment对象,几乎适用于任何情况,我们只需要一个counts矩阵即可。 🤗 但有时候我们获取的文件是cellranger(用于10X Chromium数据)的输出文件,这个时候我们可以用DropletUtils包中的read10xCounts函数。😀

sce <- read10xCounts("./pbmc_1k_raw")

sce <- read10xCounts("./pbmc_1k_filtered")
sce

6查看数据

6.1 查看counts矩阵

这里需要说一下的是,如果你的矩阵就叫counts,那么counts(tung) = assay(tung, "counts")

counts(tung)
assay(tung, "counts")[1:4,1:4]

6.2 查看cell信息

colData(tung)

6.3 查看batch信息

table(tung$batch)

🧐 补充!

Function Description
rowData(sce) gene信息
colData(sce) cell信息
assay(sce, "counts") 查看counts矩阵
reducedDim(sce, "PCA") 查看PCA降维table
sce$colname 查看colDatacolname, 等价于colData(sce)$colname
sce[<rows>, <columns>] 查看特定

7修改SingleCellExperiment对象

举个栗子 🌰
我们把counts矩阵进行log2(x+1)转换,命名为logcounts

assay(tung, "logcounts") <- log2(counts(tung) + 1)
# 查看一下
logcounts(tung)[1:4, 1:4]

🧐 补充!

Code Description
assay(sce, "name") <- matrix 添加新矩阵
rowData(sce) <- data_frame 替换rowData
colData(sce) <- data_frame 替换colData
colData(sce)$column_name <- values colData中添加/替换一个新的列
rowData(sce)$column_name <- values rowData中添加/替换一个新的行
reducedDim(sce, "name") <- matrix 添加降维矩阵

8基本数据处理

8.1 举个栗子🌰

计算我们的数据集中每个细胞)的平均counts
然后我们将它加到column metadata中作为新的一。😘

colData(tung)$mean_counts <- colMeans(counts(tung))

这个时候我们可以查看一下,多了一列mean_counts。🤔

colData(tung)[1:4,]

8.2 补充一下

这里我们再补充一下基本函数。有兴趣的小伙伴都试试吧。🫶

  • row (feature) summaries
rowSums(counts(tung))  # sum
rowMeans(counts(tung)) # mean
rowSds(counts(tung)) # standard deviation
rowVars(counts(tung)) # variance
rowIQRs(counts(tung)) # inter-quartile range
rowMads(counts(tung)) # mean absolute deviation

  • column (sample) summaries
colSums(counts(tung))  # sum
colMeans(counts(tung)) # mean
colSds(counts(tung)) # standard deviation
colVars(counts(tung)) # variance
colIQRs(counts(tung)) # inter-quartile range
colMads(counts(tung)) # mean absolute deviation

9CPM的计算与查看

9.1 CPM的计算

这里我们直接使用sweep函数来搞定countsCPM,在这之前我们先计算一下total counts并新增一列。

colData(tung)$total_counts <- colSums(counts(tung))
assay(tung, "cpm") <- sweep(counts(tung),2,tung$total_counts/1e6,'/')

# check一下
colSums(cpm(tung))[1:10]

9.2 CPM矩阵的查看

cpm(tung)[1:4, 1:4]
assay(tung, "cpm")[1:4, 1:4]

10常用过滤方法

拿到矩阵后,我们一般需要过滤一些低表达或异常表达值。
举个栗子 🌰,假设细胞至少25000counts基因在至少一半的细胞中拥有超过5个counts

10.1 过滤细胞

cell_filter <- colSums(counts(tung)) >= 25000

# check一下
table(cell_filter)

10.2 过滤基因

gene_filter <- rowSums(counts(tung) > 5) > ncol(tung)/2

# check一下
table(gene_filter)

10.3 最终过滤

tung_filtered <- tung[gene_filter, cell_filter]

tung_filtered

10.4 补充一下

colSums(counts(sce)) > x # 对细胞进行过滤, Total counts per cell greater than x。

colSums(counts(sce) > x) > y # 对细胞进行过滤,Cells with at least y genes having counts greater than x.

rowSums(counts(sce)) > x # 对基因进行滤过,Total counts per gene greater than x.

rowSums(counts(sce) > x) > y # 对基因进行过滤,Genes with at least y cells having counts greater than x.

最后祝大家早日不卷!~

需要示例数据的小伙伴,在公众号回复scRNAseq获取吧!

点个在看吧各位~ ✐.ɴɪᴄᴇ ᴅᴀʏ 〰

分类:

后端

标签:

后端

作者介绍

jamesbang
V1

wx🔍: Grassssss 卷起来了