不缺好文章、idea的不要进！
大家好，我是晨曦，本次将开启一个全新的系列，依旧是单细胞，依旧是熟悉的晨曦解读，只不过这一次我们选择的课程是Single cell RNA-seq analysis workshop 这个课程是由瑞典国家生物信息学基础设施(NBIS) 负责制作的，课程的链接在下方
Single cell RNA-seq analysis course (nbisweden.github.io)
当然了，通过这门课程我们要把我们前段时间学习的单细胞再度复习以及巩固一遍，以迎接我们下一阶段的文献解读和文献复现环节。
很可惜的是，这门课程我发现的时候已经过了申请时间，但是NBIS还是十分友好的提供了全套脚本来供我们学习，这种分享知识的精神一直是晨曦颇为敬佩且极力践行的。
晨曦单细胞笔记系列传送门
1. 首次揭秘！不做实验也能发10 SCI，CNS级别空间转录组套路全解析（附超详细代码！）
2. 过关神助！99%审稿人必问，多数据集联合分析，你注意到这点了吗？
3. 太猛了！万字长文单细胞分析全流程讲解，看完就能发文章！建议收藏！（附代码）
4. 秀儿！10 生信分析最大的难点在这里！30多种方法怎么选？今天帮你解决！
5. 图好看易上手！没有比它更适合小白入手的单细胞分析了！老实讲，这操作很sao！
7. 我就不信了，生信分析你能绕开这个问题！今天一次性帮你解决！
首先介绍一下课程大纲
晨曦解读
本课程架构
1.当前scRNA-seq技术介绍（在线）
2.原始测序数据转换为基因表达矩阵的介绍（在线）
3.scRNA-seq数据的质量控制
4.数据降维和聚类分析
5.数据归一化处理
6.scRNA-seq数据的差异基因表达
7.细胞类型鉴定
8.细胞轨迹分析
9.对比Seurat、Scran和Scanpy分析流程
了解完课程结构，那我们就开始进行今天的内容吧，今天的内容为——数据质控（QC）
#本教程使用的数据是3位COVID-19患者和3位健康正常人，总共6个PBMC的10×数据，且每个样本均被二次采样，所以每个样本为1500个细胞#示例数据已经放在百度云，评论区留言获取即可#准备工作#R包加载的提示信息如果觉得乱，以后可以使用这个函数来加载R包suppressMessages(require(Seurat))suppressMessages(require(Matrix))suppressMessages(require(DoubletFinder))晨曦解读
其实有时候有些R包安装提示也很有用，所以不要盲目使用
suppressMessages函数
#准备输入文件cov.15 <- Seurat::Read10X_h5(filename = 'raw/nCoV_PBMC_15.h5', use.names = T)cov.1 <- Seurat::Read10X_h5(filename = 'raw/nCoV_PBMC_1.h5', use.names = T)cov.17 <- Seurat::Read10X_h5(filename = 'raw/nCoV_PBMC_17.h5', use.names = T)ctrl.5 <- Seurat::Read10X_h5(filename = 'raw/Normal_PBMC_5.h5', use.names = T)ctrl.13 <- Seurat::Read10X_h5(filename = 'raw/Normal_PBMC_13.h5', use.names = T)ctrl.14 <- Seurat::Read10X_h5(filename = 'raw/Normal_PBMC_14.h5', use.names = T)
晨曦解读
导入数据成功,然后我们可以通过dim函数简单看一下数据
#很遗憾我们这样还只是把h5数据导入进入R中，还得需要我们创建Seurat对象sdata.cov15 <- CreateSeuratObject(cov.15, project = 'covid_15')sdata.cov1 <- CreateSeuratObject(cov.1, project = 'covid_1')sdata.cov17 <- CreateSeuratObject(cov.17, project = 'covid_17')sdata.ctrl5 <- CreateSeuratObject(ctrl.5, project = 'ctrl_5')sdata.ctrl13 <- CreateSeuratObject(ctrl.13, project = 'ctrl_13')sdata.ctrl14 <- CreateSeuratObject(ctrl.14, project = 'ctrl_14')#创建速度很快~#然后我们需要添加分组信息来让我们合并数据集的时候知道数据集类型sdata.cov1$type = 'Covid'sdata.cov15$type = 'Covid'sdata.cov17$type = 'Covid'sdata.ctrl5$type = 'Ctrl'sdata.ctrl13$type = 'Ctrl'sdata.ctrl14$type = 'Ctrl'晨曦解读
还记得我们的分组信息添加到哪里了吗？
如果忘了，就赶快去复习吧
#合并6个样本alldata <- merge(x = sdata.cov15, y = c(sdata.cov1, sdata.cov17, sdata.ctrl5, sdata.ctrl13, sdata.ctrl14), add.cell.ids = c('covid_15', 'covid_1', 'covid_17', 'ctrl_5', 'ctrl_13', 'ctrl_14'))#整理一下分析环境rm(cov.15, cov.1, cov.17, ctrl.5, ctrl.13, ctrl.14, sdata.cov15, sdata.cov1, sdata.cov17, sdata.ctrl5, sdata.ctrl13, sdata.ctrl14)晨曦解读
然后我们来探究一下add.cell.ids参数
我分别运行下面两句代码
#代码一alldata <- merge(x = sdata.cov15, y = c(sdata.cov1, sdata.cov17, sdata.ctrl5, sdata.ctrl13, sdata.ctrl14), add.cell.ids = c('covid_15', 'covid_1', 'covid_17', 'ctrl_5', 'ctrl_13', 'ctrl_14'))#代码二shiyan <- merge(sdata.cov15, c(sdata.cov1, sdata.cov17, sdata.ctrl5, sdata.ctrl13, sdata.ctrl14))发现得到的数据集大小是一样的，说明数据集本身数量并没有改变，可能是数据集内部有一些微调，所以我探究数据集发现
代码一
代码二
晨曦解读
所以这个参数的作用是在细胞ID前添加分组信息，但是要注意，信息的顺序需要和合并的顺序保持一致
#清空一下内存gc()#简单探索一下表达矩阵和细胞信息as.data.frame(alldata@assays$RNA@counts[1:10, 1:2])head(alldata@meta.data, 10)
#质控#计算线粒体基因百分比#这里提供两种方法#方法一（自动）alldata <- PercentageFeatureSet(alldata, '^MT-', col.name = 'percent_mito')#方法二（手动）total_counts_per_cell <- colSums(alldata@assays$RNA@counts)head(total_counts_per_cell)#查询线粒体相关基因mito_genes <- rownames(alldata)[grep('^MT-', rownames(alldata))]head(mito_genes，10)#计算线粒体基因含量alldata$percent_mito <- colSums(alldata@assays$RNA@counts[mito_genes, ])/total_counts_per_cell晨曦解读
两种方法均可，方法一是Seurat的内置算法，方法二则是纯纯的数学计算
#同样的方法计算核糖体基因百分比# Way1: Doing it using Seurat function（自动）alldata <- PercentageFeatureSet(alldata, '^RP[SL]', col.name = 'percent_ribo')# Way2: Doing it manually（手动）ribo_genes <- rownames(alldata)[grep('^RP[SL]', rownames(alldata))]head(ribo_genes, 10)alldata$percent_ribo <- colSums(alldata@assays$RNA@counts[ribo_genes, ])/total_counts_per_cell晨曦解读
常规scRNA-seq计算线粒体基因比例即可，但是也可以计算核糖体基因含量，然后因为这个是PBMC，所以我们下面计算血红蛋白相关基因含量
QC这一步骤，如果想要做细会有很多方面，但是常规的就这些，按照常规处理数据即可
#计算血红蛋白基因含量# Percentage hemoglobin genes - includes all genes starting with HB except HBP.alldata <- PercentageFeatureSet(alldata, '^HB[^(P)]', col.name = 'percent_hb')alldata <- PercentageFeatureSet(alldata, 'PECAM1|PF4', col.name = 'percent_plat')晨曦解读
最后我们Seurat对象中多出这些数据
#下面把我们的QC指标进行可视化feats <- c('nFeature_RNA', 'nCount_RNA', 'percent_mito', 'percent_ribo', 'percent_hb')VlnPlot(alldata, group.by = 'orig.ident', features = feats, pt.size = 0.1, ncol = 3) NoLegend()晨曦解读
从第一个可视化可以很清楚的看到，有四个数据集的表达是有很明显的不同的，而且随着检测到核糖体表达增高，那么构成我们最后的转录景观，核糖体基因的表达比例就会增大
这里面单纯从数值上来看就不是很准确了（基因-2000左右；RNA-2万到4万左右，线粒体比例最大在25%以下哎）
#绘制点图FeatureScatter(alldata, 'nCount_RNA', 'nFeature_RNA', group.by = 'orig.ident', pt.size = 0.5)#标准过滤#保留至少检测到200个基因的细胞且需要在至少3个细胞中表达这些基因#这个过滤条件需要结合scRNA文库制备的方法#过滤至少检测到200个基因的细胞selected_c <- WhichCells(alldata, expression = nFeature_RNA > 200)#过滤至少在3个细胞中表达的基因selected_f <- rownames(alldata)[Matrix::rowSums(alldata) > 3]#提取data.filt <- subset(alldata, features = selected_f, cells = selected_c)dim(data.filt)#[1] 18147 7973
晨曦解读
这里我们还需要注释，检测到极高基因表达的细胞，还有可能是双重细胞即一个微珠内有两个或两个以上的细胞，所以我们要进行doublet预测
但是在这个教程中跳过了这一步骤，一般来说成熟的技术是很少会出现这种情况的，同时结合前面的QC可视化也没有这种情况，但是代码放在下面
# skip for now and run DoubletFinder first!# high.det.v3 <- WhichCells(data.filt, expression = nFeature_RNA > 4100)# high.det.v2 <- WhichCells(data.filt, expression = nFeature_RNA > 2000 &# orig.ident == 'v2.1k')# remove these cells data.filt <- subset(data.filt,# cells=setdiff(WhichCells(data.filt),c(high.det.v2,high.det.v3)))# check number of cellsncol(data.filt)晨曦解读
然后我们还可以查看哪些基因对表达贡献大，绘制每个基因的表达百分比
# Compute the relative expression of each gene per cell Use sparse matrix# operations, if your dataset is large, doing matrix devisions the regular way# will take a very long time.par(mar = c(4, 8, 2, 1))C <- data.filt@assays$RNA@countsC <- Matrix::t(Matrix::t(C)/Matrix::colSums(C)) * 100most_expressed <- order(apply(C, 1, median), decreasing = T)[20:1]boxplot(as.matrix(t(C[most_expressed, ])), cex = 0.1, las = 1, xlab = '% total count per cell', col = (scales::hue_pal())(20)[20:1], horizontal = TRUE)
晨曦解读
根据上图的可视化，我们可以很清楚的看到，MALAT1基因占比很高，其它占比高的基因包括线粒体基因和核糖体基因，但是线粒体和核糖体基因含量高是普遍现象，所以MALAT1基因含量高可能是技术问题或者是生物学问题，所以需要我们明确该Gene的信息，以此来帮助我们进行下游的分析。
基于我们本身的生物学知识以及前面的QC可视化，我们选择线粒体或核糖体含量作为我们后续筛选的阈值，我们选择过滤掉线粒体读数低于20%（因为大部分细胞小于这个数值），核糖体读数小于5%的细胞来作为我们后续分析的对象。
selected_mito <- WhichCells(data.filt, expression = percent_mito < 0.2)selected_ribo <- WhichCells(data.filt, expression = percent_ribo > 0.05)# and subset the object to only keep those cellsdata.filt <- subset(data.filt, cells = selected_mito)data.filt <- subset(data.filt, cells = selected_ribo)dim(data.filt)table(data.filt$orig.ident)#[1] 18147 5762晨曦解读
这里需要说一个平时不容易让人注意的点
自己计算的线粒体基因以及核糖体基因的比例与Seurat包内置算法计算的比例相差比较大。
下面是我用Seurat内置算法计算的比例，也是走到这一步发现一共就剩下34个细胞，差距比较大，所以这也给我们自己日后分析有个提示。
如果后期在质控后细胞过少，不如尝试一下这一种可能
#回到主题table(data.filt$orig.ident)#covid_1 covid_15 covid_17 ctrl_13 ctrl_14 ctrl_5 # 878 585 1042 1154 1063 1040#接下来绘制质控后的QC可视化feats <- c('nFeature_RNA', 'nCount_RNA', 'percent_mito', 'percent_ribo', 'percent_hb')VlnPlot(data.filt, group.by = 'orig.ident', features = feats, pt.size = 0.1, ncol = 3) NoLegend()晨曦解读
接下来，根据官网解释，由于线粒体和MALAT1基因的高表达被认为是技术原因，所以在后期分析之前，剔除掉这些基因是比较明智的
dim(data.filt)#[1] 18147 5762# Filter MALAT1（过滤MALAT1基因）data.filt <- data.filt[!grepl('MALAT1', rownames(data.filt)), ]# Filter Mitocondrial（过滤线粒体基因）data.filt <- data.filt[!grepl('^MT-', rownames(data.filt)), ]# Filter Ribossomal gene (optional if that is a problem on your data) data.filt（过滤核糖体基因）# <- data.filt[ ! grepl('^RP[SL]', rownames(data.filt)), ]# Filter Hemoglobin gene (optional if that is a problem on your data)（过滤血红蛋白基因）data.filt <- data.filt[!grepl('^HB[^(P)]', rownames(data.filt)), ]dim(data.filt)#[1] 18121 5762晨曦解读
一般来看不过滤核糖体基因（但是不是绝对）
接下来我们就要查看样本的性别信息，在理想情况下，应将实验设置为单性别，这样做可以避免在结论中出现性别偏见，但是，这个在实验中很难实现。
通过查看来自Y染色体（男）和XIST基因表达（女）的读数，我们可以判断每个样本的性别，这里有一个技巧，如果样本性别与我们预测的不符，那么也很有可能是样本存在混合（也就是样本处理存在问题）。
这里，我们下载公共数据的基因注释文件，请将其放置在正确的位置以使路径genes.file发挥作用。
#读取gene.files文件#这个我在百度云已经存放好了，大家直接读取就可以genes.file = 'genes.table.csv'#如果不想本地读取，也可以在github上下载，就运行下面代码if (!file.exists(genes.file)) { suppressMessages(require(biomaRt)) # initialize connection to mart, may take some time if the sites are # unresponsive. mart <- useMart('ENSEMBL_MART_ENSEMBL', dataset = 'hsapiens_gene_ensembl') # fetch chromosome info plus some other annotations genes.table <- try(biomaRt::getBM(attributes = c('ensembl_gene_id', 'external_gene_name', 'description', 'gene_biotype', 'chromosome_name', 'start_position'), mart = mart, useCache = F)) if (!dir.exists('data/results')) { dir.create('data/results') } if (is.data.frame(genes.table)) { write.csv(genes.table, file = genes.file) } if (!file.exists(genes.file)) { download.file('https://raw.githubusercontent.com/NBISweden/workshop-scRNAseq/master/labs/misc/genes.table.csv', destfile = 'data/results/genes.table.csv') genes.table = read.csv(genes.file) }} else { genes.table = read.csv(genes.file)}genes.table <- genes.table[genes.table$external_gene_name %in% rownames(data.filt), ]#上面这么长的代码其实就可以给我们把数据整理成gene.files的形式，所以还是本地读取比较好，而且以后在有类似的性别鉴定，用这个文件是一样的#这个文件其实就是染色体信息，然后让我们简单来探索一下这个数据晨曦解读
其实这个文件写的也很清楚这里就不过多解释，总而言之，现在gene table里的基因就都是我们这次分析所有的基因了（因为有了一步取子集的过程）
#然后，我们计算每个细胞中来自Y染色体的比例#第一步先得到哪些基因对应Y染色体chrY.gene = genes.table$external_gene_name[genes.table$chromosome_name == 'Y']chrY.gene#[1] 'RPS4Y1' 'AC006157.1' 'ZFY' 'ZFY-AS1' 'LINC00278' #[6] 'PCDH11Y' 'USP9Y' 'DDX3Y' 'UTY' 'TMSB4Y' #[11] 'TTTY14' 'AC010889.1' 'KDM5D' 'EIF1AY' 'RPS4Y2'#计算data.filt$pct_chrY = colSums(data.filt@assays$RNA@counts[chrY.gene, ])/colSums(data.filt@assays$RNA@counts)#绘制XIST与chrY的比例关系FeatureScatter(data.filt, feature1 = 'XIST', feature2 = 'pct_chrY')#可视化VlnPlot(data.filt, features = c('XIST', 'pct_chrY'))晨曦解读
然后我们可以从图片上很清楚的看到样本的性别
然后这里作者提出了疑问：
1.这样的性别比例会对下游分析造成影响么？
2.如果有影响如何解决？
欢迎大家在评论区讨论，下期我将在评论区说一下我的看法
晨曦解读
接下来我们对细胞周期进行评分
这里扩展一下：为什么我们要对细胞周期进行评分呢？
简单粗暴的总结：相同类型的细胞在聚类的时候会因为细胞周期的不同而被分开，这显然是我们不想看到的，所以我们要明确细胞周期相关的基因以及打分
# Before running CellCycleScoring the data need to be normalized and# logtransformed.data.filt = NormalizeData(data.filt)#使用函数对细胞周期进行评分data.filt <- CellCycleScoring(object = data.filt, g2m.features = cc.genes$g2m.genes, s.features = cc.genes$s.genes) #可视化VlnPlot(data.filt, features = c('S.Score', 'G2M.Score'), group.by = 'orig.ident', ncol = 4, pt.size = 0.1)晨曦解读
看上去不需要考虑细胞周期的影响，如果不放心，后期在高可变基因筛选完以后，还可以做一个细胞周期删除前后的聚类详细再看一下~
接下来，我们要预测双细胞。
其实双细胞就是我们的微珠内有两个细胞，但是往往我们都是没有这一步骤的，因为10×在数据质控方面做的还是可以的，但是也不妨碍我们掌握这个技术。
其中双细胞的比例与加载细胞的数量呈线性相关。
官网原文：
大多数双细胞检测工具都是通过合并细胞的计数来模拟双细胞，并将预测到的与模拟双细胞相似的细胞归为doublets。
大多数此类软件包都需要对数据集中预期双细胞数的数量/比例进行假设。这里使用的数据是经过二次抽样的，但是原始数据集中每个样本包含了大约5000个细胞，因此我们可以假设它们加载了大约9000个细胞，并且其双细胞率约为4％（其实就是因为有双细胞的存在导致我们样本和得到的数据的细胞数不匹配）。
理想情况下，我们应该对每个样本分别进行doublet预测，尤其是在不同的样本具有不同比例的细胞类型的情况下。
但是在这里，因为我们对数据进行了二次采样，因此每个样本中只有很少的细胞，并且所有样本都是经过sorted的PBMC，因此可以将它们一起运行（也就是每个样本细胞少的时候是可以的，但是如果多的话，需要每个样本单独运行，也就是提取每个组的细胞）
# remotes::install_github('chris-mcginnis-ucsf/DoubletFinder')suppressMessages(require(DoubletFinder))# 筛选高变基因data.filt = FindVariableFeatures(data.filt, verbose = F)# 进行数据归一化data.filt = ScaleData(data.filt, vars.to.regress = c('nFeature_RNA', 'percent_mito'), verbose = F)# PCA降维data.filt = RunPCA(data.filt, verbose = F, npcs = 20)# UMAP可视化data.filt = RunUMAP(data.filt, dims = 1:10, verbose = F)# Can run parameter optimization with paramSweep, but skip for now.# sweep.res <- paramSweep_v3(data.filt) # sweep.stats <- summarizeSweep(sweep.res, GT = FALSE) # bcmvn <- find.pK(sweep.stats) # barplot(bcmvn$BCmetric, names.arg = bcmvn$pK, las=2)# define the expected number of doublet cellscells.nExp <- round(ncol(data.filt) * 0.04) # expect 4% doublets# 使用doubletFinder_v3函数预测双细胞data.filt <- doubletFinder_v3(data.filt, pN = 0.25, pK = 0.09, nExp = nExp, PCs = 1:10)## [1] 'Creating 1921 artificial doublets...'## [1] 'Creating Seurat object...'## [1] 'Normalizing Seurat object...'## [1] 'Finding variable genes...'## [1] 'Scaling data...'## [1] 'Running PCA...'## [1] 'Calculating PC distance matrix...'## [1] 'Computing pANN...'## [1] 'Classifying doublets..'# name of the DF prediction can change, so extract the correct column name.DF.name = colnames(data.filt@meta.data)[grepl('DF.classification', colnames(data.filt@meta.data))]cowplot::plot_grid(ncol = 2, DimPlot(data.filt, group.by = 'orig.ident') NoAxes(), DimPlot(data.filt, group.by = DF.name) NoAxes())晨曦解读
正常情况下，双细胞中会获得比单细胞更多的基因，我们可视化一下，看我们预测的结果是否准确
VlnPlot(data.filt, features = 'nFeature_RNA', group.by = DF.name, pt.size = 0.1)晨曦解读
现在我们可以移除这些细胞
data.filt = data.filt[, data.filt@meta.data[, DF.name] == 'Singlet']dim(data.filt)#[1] 18121 5532
晨曦解读
至此，数据质控完毕~
dir.create('results', showWarnings = F)saveRDS(data.filt, 'results/seurat_covid_qc.rds')