seurat提取表达矩阵_（单细胞-SingleCell）Seurat流程文献复现

109 阅读 0 评论 72 点赞

我是靠谱客的博主微笑钥匙，最近开发中收集的这篇文章主要介绍seurat提取表达矩阵_（单细胞-SingleCell）Seurat流程文献复现，觉得挺不错的，现在分享给大家，希望可以做个参考。

概述

单细胞项目：来自于30个病人的49个组织样品，跨越3个治疗阶段

Therapy-Induced Evolution of Human Lung Cancer Revealed by Single-Cell RNA Sequencing

这篇教程我将分为四个阶段完整的阐述单细胞的主流的下游分析流程

数据预处理（数据准备阶段）
seurat 基础分析
免疫细胞识别
inferCNV 的实现

先来第一部分数据预处理

这里的数据预处理很自由，其中一部分是必须按照严格的标准进行比如计算外源基因，线粒体，核糖体等等

# 导入包
library(Seurat)
library(dplyr)
library(Matrix)
library(magrittr)
library(stringr)
library(tidyverse)
options(stringsAsFactors = F)
options(as.is = T)
setwd('/Users/yuansh/Desktop/单细胞/scell_lung_adenocarcinoma-master')

文章的话一共是有 2 个表达矩阵需要处理

我们这里的话要很明确的知道表达谱的预处理到底预处理什么东西

看一下 count 矩阵是否是 hista 比对结果，如果是需要剔除 hista 比对的注释信息。（有的也是已经处理完了，但是也得检查）
确定一下 cell 名的组成，构建好 metadata（行名）和 count（列名）矩阵能够对应
过滤掉外源 RNA （ERCC）
合并成为 seurat 对象
计算线粒体（MT），核糖体的 RNA（^RP[SL]）比例并筛选数据

我们先处理第一个表达矩阵

# 这个是治疗前的第一个表达矩阵
# 导入数据
raw.data = read.csv('csv_files/S01_datafinal.csv', header=T, row.names = 1)
metadata = read.csv("csv_files/S01_metacells.csv", row.names=1, header=T)
dim(raw.data)
dim(metadata)

判断一下是否是 hista 比对的结果并且没有处理，如果是的话后面 5 行是没用的
然后我们发现并不是
继续观察发现,raw.data 是 count 矩阵,metadata 是信息矩阵
count 矩阵的列名是 meta 矩阵的 well 和 plate 合并的
稍微对数据进行处理一下,然后保存，确保 count 矩阵的列名和 metadata 的行名一样
这样的话下次读取就很快

> tail(raw.data[,1:4])
A10_1001000329 A10_1001000362 A10_1001000366 A10_1001000367
ZXDC                0              0              0              0
ZYG11A              0              0              0              0
ZYG11B              0              0              0              0
ZYX               292              0              0              0
ZZEF1              66              0              0              0
ZZZ3                0              0              0              0

> raw.data[1:4,1:4]
A10_1001000329 A10_1001000362 A10_1001000366 A10_1001000367
A1BG                  3              0              0              0
A1BG-AS1              0              0              0              0
A1CF                  0              0              0              0
A2M                  68              0              0              0

> metadata[1:4,1:4]
well      plate        cell_id sample_name
1  A10 1001000329 A10_1001000329      LT_S07
2  A10 1001000362 A10_1001000362      LT_S12
3  A10 1001000366 A10_1001000366      LT_S13
4  A10 1001000367 A10_1001000367      LT_S13
paste0(metadata$well[1],'_',metadata$plate[1]) # 在进行大幅度修改的时候要判断代码有没有写错
colnames(raw.data)[1]
# 修改一下 metadata 的行名
rownames(metadata) = paste0(metadata$well,'_',metadata$plate)
# 确定了行名metadata 的行名和 count 矩阵 raw.data 的列名相等
identical(colnames(raw.data),rownames(metadata))

接下来就是基因的处理，按照上面说的步骤：
1. 剔除外源基因，然后构建 seurat 对象 2. 计算线粒体基因（有没有都算就行了） 3. 计算核糖体基因

# 这里要过滤 ERCC(External RNA Control Consortium)
# ERCC 是外源 RNA 主要的作用就是用来质控的
erccs = grep('^ERCC-', x= rownames(x=raw.data),value = T) # value = T 获取名字
# 计算 ercc 比例
percent.ercc = Matrix::colSums(raw.data[erccs, ])/Matrix::colSums(raw.data)
fivenum(percent.ercc) 
ercc.index = grep(pattern = "^ERCC-", x = rownames(x = raw.data), value = FALSE) # 获取 index
# 删除 Count 矩阵中的 ERCC
raw.data = raw.data[-ercc.index,]
dim(raw.data) 
# [1] 26485 27489

构建 seurat 对象

# 构建对象
# min.cells 某一个基因至少在多少个基因中表达
# min.features 某个细胞至少表达多少个基因
min.cells = 0 # 3 因为文章没有筛选，所以这里设为0
min.features = 0 # 150 # 同上
sce = CreateSeuratObject(counts = raw.data,metadata = metadata,min.cells =min.cells,min.features =min.features)
sce

我们先看一下 sce 对象中的 meta.data,这个是构建 seurat 对象自动生成的，缺少了很多信息，所以要把 metdata 添加上去

> sce@meta.data[1:5,1:5]
orig.ident nCount_RNA nFeature_RNA sample_name patient_id
A10_1001000329 SeuratProject     681122         2996      LT_S07      TH103
A10_1001000362 SeuratProject         30           17      LT_S12      TH172
A10_1001000366 SeuratProject         29           19      LT_S13      TH169
A10_1001000367 SeuratProject        286            3      LT_S13      TH169
A10_1001000372 SeuratProject        599          240      LT_S14      TH153
sce = AddMetaData(object = sce, metadata = metadata)
sce = AddMetaData(object = sce, percent.ercc, col.name = "percent.ercc")
# Head to check
head(sce@meta.data)

然后就是数据清洗
数据清洗原则： 1. 不符合质量的：这里一般指细胞中基因数太少或者某些基因只在两三个细胞中表达 2. 过滤线粒体 DNA占比过高的 3. 过滤核糖体占比过高的在过滤前先计算

table(grepl("^MT-",rownames(sce)))
table(grepl("^RP[SL][[:digit:]]",rownames(sce)))
sce[["percent.mt"]]  = PercentageFeatureSet(sce, pattern = "^MT-")
sce[["percent.rp"]]  = PercentageFeatureSet(sce, pattern = "^RP[SL][[:digit:]]")
dim(sce)
summary(sce[["nCount_RNA"]])
summary(sce[["nFeature_RNA"]])
summary(sce[["percent.mt"]]  )
summary(sce[["percent.ercc"]] )
summary(sce[["percent.rp"]]  )

文章中是啥都没去，就正常的过滤了一下基因和表达过低的细胞

# 过滤筛选
# 一般根据实际情况进行筛选，这次的筛选主要还是根据文章里面的流程
# 这个过滤蛮自由的，所以随便搞一下就行
myData1 = subset(x=sce, subset = nCount_RNA > 50000 & nFeature_RNA > 500)
myData1
save(myData1,file='myData1.Rdata')

第一阶段的处理流程就是这样

接着处理第二个表达矩阵

# 这个是治疗前的第一个表达矩阵
# 导入数据
raw.data = read.csv('Data_input/csv_files/neo-osi_rawdata.csv', header=T, row.names = 1)
metadata = read.csv("Data_input/csv_files/neo-osi_metadata.csv", row.names=1, header=T)
raw.data[1:4,1:4]
metadata[1:4,1:4]

判断一下是否是 hista 比对的结果并且没有处理，如果是的话后面 5 行是没用的

然后发现果然是

继续观察发现,第二套数据的命名规则和第一套有点不太一样,多了_S10.homo,因此要把他去掉

稍微对数据进行处理一下,然后保存，确保 count 矩阵的列名和 metadata 的行名一样这样的话下次读取就很快

> tail(raw.data[,1:4])
A10_B000561_S10.homo A10_B000563_S10.homo A10_B000568_S10.homo A10_B001544_S10.homo
ZZZ3                                      0                    0                    0                    0
__no_feature                          31084                92098               131276              1197071
__ambiguous                             353                 2231                  570                21247
__too_low_aQual                           0                    0                    0                    0
__not_aligned                             0                    0                    0                    0
__alignment_not_unique