urat提取表达矩阵_（单细胞-SingleCell）Seurat流程⽂献复
现
单细胞项⽬：来⾃于30个病⼈的49个组织样品，跨越3个治疗阶段
Therapy-Induced Evolution of Human Lung Cancer Revealed by Single-Cell RNA Sequencing
熊猫历险记这篇教程我将分为四个阶段完整的阐述单细胞的主流的下游分析流程
1. 数据预处理（数据准备阶段）
2. urat 基础分析
3. 免疫细胞识别
4. inferCNV 的实现
先来第⼀部分数据预处理
这⾥的数据预处理很⾃由，其中⼀部分是必须按照严格的标准进⾏⽐如计算外源基因，线粒体，核糖体等等
# 导⼊包
library(Seurat)
library(dplyr)
library(Matrix)
library(magrittr)
library(stringr)
library(tidyver)
options(stringsAsFactors = F)
options(as.is = T)
twd('/Urs/yuansh/Desktop/单细胞/scell_lung_adenocarcinoma-master')
⽂章的话⼀共是有 2 个表达矩阵需要处理
我们这⾥的话要很明确的知道表达谱的预处理到底预处理什么东西
1. 看⼀下 count 矩阵是否是 hista ⽐对结果，如果是需要剔除 hista ⽐对的注释信息。（有的也是已经处理完了，但是也得检查）
2. 确定⼀下 cell 名的组成，构建好 metadata（⾏名） 和 count（列名） 矩阵能够对应
3. 过滤掉外源 RNA （ERCC）
4. 合并成为 urat 对象
5. 计算线粒体（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. 剔除外源基因，然后构建 urat 对象
2. 计算线粒体基因（有没有都算就⾏了）
3. 计算核糖体基因
# 这⾥要过滤 ERCC(External RNA Control Consortium)
# ERCC 是外源 RNA 主要的作⽤就是⽤来质控的
erccs = grep('^ERCC-', x= rownames(x=raw.data),value = T) # value = T 获取名字
# 计算 ercc ⽐例
< = Matrix::colSums(raw.data[erccs, ])/Matrix::colSums(raw.data)
)
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
构建 urat 对象
# 构建对象
# lls 某⼀个基因⾄少在多少个基因中表达
# min.features 某个细胞⾄少表达多少个基因
和的多音字组词min.features = 0 # 150 # 同上
sce = CreateSeuratObject(counts = raw.data,metadata = lls =lls,min.features =min.features)
sce
我们先看⼀下 sce 对象中的 meta.data，这个是构建 urat 对象⾃动⽣成的，缺少了很多信息，所以要把 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, , col.name = "")
# Head to check
head(sce@meta.data)
然后就是数据清洗
数据清洗原则： 1. 不符合质量的：这⾥⼀般指细胞中基因数太少或者某些基因只在两三个细胞中表达 2. 过滤线粒体 DNA占⽐过⾼的 3. 过滤核糖体占⽐过⾼的 在过滤前先计算
table(grepl("^MT-",rownames(sce)))
table(grepl("^RP[SL][[:digit:]]",rownames(sce)))
sce[[""]] = 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[[""]] )
summary(sce[[""]] )
summary(sce[["percent.rp"]] )
⽂章中是啥都没去，就正常的过滤了⼀下基因和表达过低的细胞
# 过滤筛选
# ⼀般根据实际情况进⾏筛选，这次的筛选主要还是根据⽂章⾥⾯的流程
# 这个过滤蛮⾃由的，所以随便搞⼀下就⾏
myData1 = subt(x=sce, subt = 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 67598 44778 47959 500993
> raw.data[1:4,1:4]
A10_B000561_S10.homo A10_B000563_S10.homo A10_B000568_S10.homo A10_B001544_S10.homo
A1BG 0 0 0 0
A1BG-AS1 0 0 0 0
A1CF 0 0 0 0
A2M 0 0 0 0
> metadata[1:4,1:4]
plate sample_name patient_id sample_type
1 B001567 AZ_01 AZ003 tumor
2 B001544 AZ_01 AZ00
3 tumor
3 B001546 AZ_01 AZ003 tumor
4 B001481 AZ_01 AZ003 tumor
然后我们进⼀步的发现，这套数据好像是没有 well 信息的 所以⼜要进⼀步处理⼀下，确保count 矩阵的列名和 metadata 矩阵的⾏名可以对应
gsub('_S[[:digit:]]*.homo','',colnames(raw.data))[1:5]
colnames(raw.data) = gsub('_S[[:digit:]]*.homo','',colnames(raw.data))
raw.data = raw.data[-grep('__',rownames(raw.data)),]
tail(raw.data[,1:4])
metadata$well[1] # 发现这套数据中的 metadata 中没有 well 信息
metadata$plate[1] # 好在有 plate信息
所以稍微进⾏处理⼀下
# 将 count 矩阵的名字拆开⾃⼰构建 well 信息矩阵
library(stringr)
wellInfo = str_split(colnames(raw.data),'_',simplify = T)%>% as.data.frame()
colnames(wellInfo) = c("well", "plate")
rownames(wellInfo) = paste0(wellInfo[,1],'_',wellInfo[,2])
# 修改⼀下 metadata 的⾏名
wellInfo[1:4,]
metadata = left_join(wellInfo, metadata, by = "plate")
metadata[1:4,]
wellInfo[1:4,]
rownames(metadata) = rownames(wellInfo)
# 判断⼀下是否⼀致
identical(colnames(raw.data),rownames(metadata))
> metadata[1:5,1:5]
well plate sample_name patient_id sample_type
A10_B000561 A10 B000561 AZ_04 AZ008_NAT NAT
A10_B000563 A10 B000563 AZ_04 AZ008_NAT NAT
A10_B000568 A10 B000568 AZ_05 AZ008 tumor
A10_B001544 A10 B001544 AZ_01 AZ003 tumor
A10_B001546 A10 B001546 AZ_01 AZ003 tumor
接下来的步骤和上⾯完全⼀样
# ERCC 是外源 RNA 主要的作⽤就是⽤来质控的
erccs = grep('^ERCC-', x= rownames(x=raw.data),value = T) # value = T 获取名字
# 计算 ercc ⽐例
< = Matrix::colSums(raw.data[erccs, ])/Matrix::colSums(raw.data)
) # 查看前 5 个结果
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 3507
# 构建对象
# 构建对象
# lls 某⼀个基因⾄少在多少个基因中表达
# min.features 某个细胞⾄少表达多少个基因
min.features = 0 # 150 #
sce = CreateSeuratObject(counts = raw.data,metadata = lls =lls,min.features =min.features) # 添加metadata 和外源基因信息
sce = AddMetaData(object = sce,metadata = metadata)
sce = AddMetaData(object = sce, , col.name = "")
# 2.数据清洗
# 数据清洗原则
# 1. 不符合质量的
# 2. 过滤线粒体 DNA
# 3. 过滤外原DNA
# 在过滤前先计算
table(grepl("^MT-",rownames(sce)))
table(grepl("^RP[SL][[:digit:]]",rownames(sce)))
sce[[""]] = 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[[""]] )
summary(sce[[""]] )
summary(sce[["percent.rp"]] )
# 过滤筛选
# ⼀般根据实际情况进⾏筛选，这次的筛选主要还是根据⽂章⾥⾯的流程
学期家长寄语myData2 = subt(x=sce, subt = nCount_RNA > 50000 & nFeature_RNA > 500)
myData2
# 26485 features across 2070 samples within 1 assay
save(myData2,file='myData2.Rdata')
第⼀部分结束，准备进⼊第⼆部分
第⼆部分 urat 标准流程
⼩提琴图观察数据分布情况
rm(list = ls())
# ⼩提琴图可视化
load("myData1.Rdata")
load("myData2.Rdata")
allData = merge(x = myData1, y = myData2)
table(allData@meta.data$sample_name) %>% length()
VlnPlot(allData, features = 'nFeature_RNA', group.by = 'sample_name')
