上下游，合体！

生信菜鸟团

发布于 2023-08-23 09:14:23

2500

发布于 2023-08-23 09:14:23

文章被收录于专栏：生信菜鸟团

这个新专辑有以下几点希冀：

带着像我一样的单细胞小白，一步步利用我们生信技能树、生信菜鸟团、单细胞天地的资源，掌握基本的scRNAseq流程
在学习的过程中，探索出合适的学习路径，帮助大家更好地利用已有资源
对过往推文中出现的错误、更新的软件进行审查，推陈出新
在过去的基本内容上深入挖掘影响小白学习的障碍，提炼总结，拓宽深度宽度
和大家讨论我在从零开始学习过程中遇到的问题，老师们在评论区指出我的不足提出建议

而我在将自己的学习笔记排版成推文时也会遵循以下行文特点：

务必详实逐步复现，如展示原推文中没展示的过程结果，添加参考资料帮助理解
重点推陈出新，如果原推文足够详细且我没遇到其他问题，可能会直接带过这篇学习推文，只在推文中展示结果，但是仍会告诉大家我看了啥，以便梳理小白学习路径

前面我们详细地分开学习了单细胞上游、下游，这次我们就学习曾老师的一个完整的项目来实战前面学过的基本流程，并适度拓宽。其实一开始我不打算把我自己学习的这个实战贴出来，想着只在新推文开头告诉大家我拿这个练手了，学学其他的，奈何我的单细胞之路处处坑，曾老师原推文已经很详细了，我都能菜出“新花样”强行变成推文素材。

原推文：小鼠的5个样品的10x技术单细胞转录组上游定量（文末赠送全套代码）

在本专辑第一篇推文，我就提到了这个实战推文，也算是来拔掉当时立的flag

下载数据

一上来我就踩了坑

创建kingfisher环境时忘记回到base了

现在在sc-RNAseq环境下搞了个内部环境，base也变了

手动修改.bashrc配置文件，顺便把sc-RNAseq下的mamba卸载了，以除后患

所以大家使用的时候还是需要注意退到base环境下再装mamba避免麻烦

正确在base环境下装mamba后.bashrc配置文件就会变成这样：

之后就可以用mamba来管理环境了，与conda兼容很好一般命令不会发生变化，只是把conda换成mamba

正式下载数据：

可以发现kingfisher很方便，只需要一个项目号，都不需要自己去ENA获取

大家可以参考本专辑第一篇推文小鼠的5个样品的10x技术单细胞转录组上游定量（文末赠送全套代码） 折腾那么久

下载成功后还需要对文件名进行修改，修改依据也主要是单细胞建库测序的基本原理 一文搞定基本cellranger定量

cellranger定量

参考文件官网下载：

关于软件版本也需要注意:

因为新版涉及到一个默认intron mode 的问题，会保留map到内含子区域的reads，参考是的，不同版本的cellranger软件对10x技术单细胞定量结果可以相差5倍以上

如果你在循环命令中加入了nohup 你会发现所有样本定量会同时进行，这就会考验机器的资源，不太给力的机器还是不要nohup，这样循环就可以 one by one

cellranger count 是一个用于处理单细胞RNA测序数据的工具，它接受多个参数来指定不同的配置和运行选项。下面是各参数的详细介绍：

--id=<string>：为每个样本指定一个唯一的ID。该可以是任意字符串，通常用于指定输出结果的文件夹名称。
--transcriptome=<path>：指定参考转录组的路径。参考转录组是一个包含参考基因组序列和注释信息的文件夹，用于将测序数据与基因组进行比对和分析。
--fastqs=<path>：指定FASTQ文件的路径。FASTQ文件包含测序数据的原始读数和其质量信息。
--sample=<>：指定当前要处理的样本的名字。该参数通常与--id参数一致。
--nosecondary：通常在度量测序数据的过程中，会生成与参考转录组不匹配的次要比对结果。使用此参数可以禁用次要比对生成，以减少内存和磁盘间的使用。
--localcores=<>：指定本地计算机的处理核心数。
--localmem=<integer>：指定本地计算机的可用内存量，单位为GB。

然后简单的整理一下结果矩阵和报表:

全部的表达量矩阵在outputs下面的matrix文件夹，而报表在outputs下面的html文件夹

后面的代码，曾老师封装了脚本，感兴趣的同学可以去原文底下获取

我们这里简单了解每个封装起来的函数基本功能，继续分析

下游分析

step1:导入数据

dir='outputs/matrix/' 
samples=list.files( dir )
samples

# samples = head(samples,10) 
sceList = lapply(samples,function(pro){ 
  # pro=samples[1] 
  print(pro)  
  sce =CreateSeuratObject(counts =  Read10X(file.path(dir,pro )) ,
                          project = gsub('^GSM[0-9]*_','',pro)  ,# pro, #
                          min.cells = 5,
                          min.features = 500 )
  return(sce)
})
names(sceList) 

sce.all=merge(x=sceList[[1]],
              y=sceList[ -1 ],
              add.cell.ids =  gsub('^GSM[0-9]*_','',samples)    )

> as.data.frame(sce.all@assays$RNA@counts[1:10, 1:2])
              SRR19904954_AAACCCAAGCCTGAAG-1 SRR19904954_AAACCCAAGTGCCAGA-1
Xkr4                                       0                              0
Gm37381                                    0                              0
Mrpl15                                     2                              1
Lypla1                                     0                              0
Tcea1                                      1                              0
Rgs20                                      0                              0
Atp6v1h                                    0                              0
Rb1cc1                                     0                              2
4732440D04Rik                              0                              0
Alkal1                                     0                              0
> head(sce.all@meta.data, 10)
                                orig.ident nCount_RNA nFeature_RNA
SRR19904954_AAACCCAAGCCTGAAG-1 SRR19904954       5505         2169
SRR19904954_AAACCCAAGTGCCAGA-1 SRR19904954       8050         2456
SRR19904954_AAACCCACAACCCTCT-1 SRR19904954       5531         2251
SRR19904954_AAACCCACAATAGGAT-1 SRR19904954       4531         1743
SRR19904954_AAACCCACAGGAATCG-1 SRR19904954      11716         3289
SRR19904954_AAACCCAGTGTCATGT-1 SRR19904954      10420         2804
SRR19904954_AAACCCATCCATCGTC-1 SRR19904954       6525         2213
SRR19904954_AAACCCATCTCTGACC-1 SRR19904954       4245         1333
SRR19904954_AAACCCATCTGATGGT-1 SRR19904954       3414         1600
SRR19904954_AAACGAAAGTCAGGGT-1 SRR19904954       5719         1996
> table(sce.all$orig.ident) 

SRR19904954 SRR19904955 SRR19904956 SRR19904957 SRR19904958 
       9084        9924        9240        6333        9804

> # 分组信息，通常是隐含在文件名，样品名字里面
> # phe=str_split( colnames(sce.all),'[-_]',simplify = T)
> # head(phe)
> # table(phe[,1])
> # sce.all$group= phe[,1]
> 
> #sce.all$group=toupper( substring(colnames(sce.all),1,5))
> # sce.all@meta.data$group =  ifelse(grepl('ne',sce.all$orig.ident),'ne','e')
> # sce.all$group=  gsub( '^[0-9]*_' ,'',sce.all$orig.ident) 
> # sce.all$group = sce.all$orig.ident
> #table(sce.all@meta.data$group)
> table(sce.all@meta.data$orig.ident)

SRR19904954 SRR19904955 SRR19904956 SRR19904957 SRR19904958 
       9084        9924        9240        6333        9804

step2:QC质控

###### step2:QC质控 ######
dir.create("./1-QC")
setwd("./1-QC")
# 如果过滤的太狠，就需要去修改这个过滤代码
source('../scRNA_scripts/qc.R')
sce.all.filt = basic_qc(sce.all)
print(dim(sce.all))  # 20513 44385
print(dim(sce.all.filt))  # 20513 43952
setwd('../')


sp='mouse'

在qc.R脚本中，创建了一个basic_qc函数

计算线粒体基因比例
计算核糖体基因比例
计算红血细胞基因比例
可视化细胞的上述比例情况

过滤（脚本中默认不执行）

根据上述指标，过滤低质量细胞/基因

过滤指标1:最少表达基因数的细胞&最少表达细胞数的基因

一般来说，在CreateSeuratObject的时候已经是进行了这个过滤操作

如果后期看到了自己的单细胞降维聚类分群结果很诡异，就可以回过头来看质量控制环节

先走默认流程即可

selected_c <- WhichCells(input_sce, expression = nFeature_RNA > 500)
selected_f <- rownames(input_sce)[Matrix::rowSums(input_sce@assays$RNA@counts > 0 ) > 3]
input_sce.filt <- subset(input_sce, features = selected_f, cells = selected_c)
dim(input_sce) 
dim(input_sce.filt)

根据代码的结果，我们可以看到：

WhichCells() 函数用于选取符合条件 nFeature_RNA > 500 的细胞
rownames(input_sce)[Matrix::rowSums(input_sce@assays$RNA@counts > 0 ) > 3] 这一行代码用于选取具有至少3个大于0计数的特征的基因
subset(input_sce, features = selected_f, cells = selected_c) 这一行代码用于根据选定的细胞和基因子集来创建一个新的 Seurat 对象 input_sce.filt。它使用了 subset() 函数来选择 input_sce 中的指定细胞和基因。

过滤指标2:线粒体/核糖体基因比例(根据violin图)
可视化过滤后的情况

step3:harmony整合多个单细胞样品

###### step3: harmony整合多个单细胞样品 ######
dir.create("2-harmony")
getwd()
setwd("2-harmony")
source('../scRNA_scripts/harmony.R')
# 默认 ScaleData 没有添加"nCount_RNA", "nFeature_RNA"
# 默认的代码都是可以自由修改的，需要对全部的代码有基本的认识。
sce.all.int = run_harmony(sce.all.filt)
setwd('../')

在harmony.R脚本中，创建了一个run_harmony函数

run_harmony对一个包含单细胞RNA测序数据的对象进行一系列处理和分析。主要步骤如下：

第一步是对输入数据进行标准化处理，使用的方法是LogNormalize。
接下来找到变异特征，这一步骤使用了FindVariableFeatures函数。
进行数据归一化，使用了ScaleData函数。
运行PCA分析，提取出变异特征的主成分，这一步骤使用了RunPCA函数。
使用RunHarmony函数对数据进行整合，通过指定"orig.ident"来整合不同样本的细胞。
运行UMAP降维，并将降维后的结果保存在seuratObj对象中。
根据UMAP降维结果，在整合后的数据上使用FindNeighbors函数建立近邻关系。
创建input_sce.all作为保存整合和降维后数据的备份。
使用不同的分辨率参数(0.01, 0.05, 0.1, 0.2, 0.3, 0.5, 0.8, 1)进行聚类分析，观察不同分辨率下的聚类效果。
分析聚类结果，输出各个聚类簇的细胞数目。
使用DimPlot函数绘制不同分辨率下的UMAP结果，以及不同分辨率下聚类结果的树状图，并保存为PDF文件。

可以发现也是我们之前学习的基本流程 初探单细胞下游

12.输出活跃标识(active.ident)的频数统计表格。

13.将整合和降维后的数据对象保存为RDS文件。

14.返回整合和降维后的数据对象。

step4: 降维聚类分群和看标记基因库

###### step4: 降维聚类分群和看标记基因库 ######
# 原则上分辨率是需要自己肉眼判断，取决于个人经验
# 为了省力，我们直接看  0.1和0.8即可
table(Idents(sce.all.int))
table(sce.all.int$seurat_clusters)
table(sce.all.int$RNA_snn_res.0.1) 
table(sce.all.int$RNA_snn_res.0.8) 

getwd()
dir.create('check-by-0.1')
setwd('check-by-0.1')
sel.clust = "RNA_snn_res.0.1"
sce.all.int <- SetIdent(sce.all.int, value = sel.clust)
table(sce.all.int@active.ident) 

source('../scRNA_scripts/check-all-markers.R')
setwd('../') 
getwd()

dir.create('check-by-0.8')
setwd('check-by-0.8')
sel.clust = "RNA_snn_res.0.8"
sce.all.int <- SetIdent(sce.all.int, value = sel.clust)
table(sce.all.int@active.ident) 
source('../scRNA_scripts/check-all-markers.R')
setwd('../') 
getwd()

last_markers_to_check

在check-all-markders.R脚本中，

①首先根据先验知识定义了许多markers对应的基因列表

# ...
markers = c('gastric_cancer_markers','lung_epi_markers',
            'Tcells_markers',
            'stromal_markers', 
            'last_markers' )
markers_list <- c(
  'GSE150580_breast_cancer_markers_list' ,
  'SCP1661_meyloids_markers_list' ,
  'myeloids_markers_list1' ,
  'myeloids_markers_list2' ,
  'CD4_markers_list' ,
  'CD8_markers_list1' ,
  'CD8_markers_list2' ,
  'cd4_and_cd8T_markers_list'   ,
  'Bcels_markers_list' ,
  'Hepatic_stellate_markers_list' 
)

UMAP降维可视化

②接着，根据不同物种（人或鼠）进行基因分析和绘图

如果物种是人类（sp=='human'）：
- 利用 lapply函数对 markers列表中的每一个元素进行循环操作。
- 在每次循环中，获取每个marker对应的基因列表，并将标记基因转换为大写字母形式并存储在"genes_to_check"中。
- 使用 DotPlot函数绘制基因的表达图，其中设置的特征为"genes_to_check"，并保存为PDF文件。
- 对生成的图形进行坐标翻转，并设置x轴标签旋转角度为45度。
- 根据"genes_to_check"的数量计算图形的高度"h"，并使用"ggsave"函数保存生成的图形为PDF文件，文件名以"check_for_"和标记基因名字为前缀。

这里以分辨率为0.1 check-by-0.1 的结果为例：

以这两个基因列表为例

可以发现绘图也可以接受列表元素，绘出来的图Features之间也分开了

如果物种是小鼠（sp=='mouse'），过程与人类类似
如果物种既不是人类也不是小鼠，则输出提示信息，表示只接受人类或小鼠的数据。

最后一个基因列表中非冗余基因的表达情况 DotPlot(sce.all.int , features = last_markers_to_check )，与前面的UMAP降维结果一块比较查看

③可视化细胞的上述比例情况

④使用COSG（Cell Marker Gene Identification）方法对给定数据集中的细胞组进行标记基因的识别

我们前面初探单细胞下游识别高变基因使用的是Seurat包自带的FindVariableFeatures函数，现在已经有了许多其它方法来探索单细胞数据集的高变基因，如COSG包

marker_cosg <- cosg(
  sce.all.int,
  groups='all',
  assay='RNA',
  slot='data',
  mu=1,
  n_genes_user=100)

使用COSG（Cell Marker Gene Identification）方法对给定数据集中的细胞组进行标记基因的识别

COSG是一种准确且快速的方法，用于在给定数据集中识别细胞组的标记基因。它可以帮助研究人员识别具有不同功能或表型的细胞组，并进一步了解它们的生物学特征。

参数说明如下：

sce.all.int是一个SingleCellExperiment对象，其中包含了RNA-seq数据。
groups='all'用于定义要分析的细胞组。在这里，我们使用了数据集中的所有细胞组。
assay='RNA'指定要分析的assay类型为RNA。
slot='data'表示从SingleCellExperiment对象的"data" slot中获取数据。
mu=1是cosg函数中的一个参数，用于调整标记基因识别的准确性和灵敏度。默认值为1。
n_genes_user=100是cosg函数中的一个参数，指定在计算标记基因时要考虑的top基因数量。

⑤可视化COSG找到的top标记基因不同分组的表达情况

分别用 DoHeatmap和 DotPlot可视化各分组top10和top3

需要注意的是，使用 SetIdent函数可以为特定细胞或整个数据集中的所有细胞设置标识类别信息。这对于对细胞进行分组或分类非常有用，可以进一步进行聚类、差异表达分析和其他下游分析

我在分析代码时拿到分辨率为0.1的结果后，发现与当前 marker_cosg变量等其他包含分组信息变量，对不上，其实是因为在代码执行的过程中，我们先运行分辨率0.1后运行的0.8，并且使用相同的变量名获取

step5: 确定单细胞亚群生物学名字

###### step5: 确定单细胞亚群生物学名字 ######
# 一般来说，为了节省工作量，我们选择0.1的分辨率进行命名
# 因为命名这个步骤是纯人工 操作
# 除非0.1确实分群太粗犷了，我们就选择0.8 
source('scRNA_scripts/lib.R')
# sce.all.int = readRDS('2-harmony/sce.all_int.rds')
# sp='mouse'
colnames(sce.all.int@meta.data) 
table(sce.all.int$RNA_snn_res.0.8)

# keratinocyte：SFN
#melanocyte1 "MLANA" 
selected_genes=c('MUC4', 'PI3', 'SIX3', # nose
                 'SCGB1A1', 'TFF3',   'KRT1',   'KRT10',  
                 'KRT5', 'TP63',   'DLK2',
                 'MKI67', 'TOP2A', 'CDC20' ,
                 'DMD','PUS7','PGAP1','IFI44L',"MLANA" ,'SFN',
                 'MUC5AC' , 'MUC5B',# secretory cells
                 'FOXJ1', 'TPPP3',  'SNTN',  # multiciliated cells 
                 'DEUP1', 'FOXN4',  'CDC20B' # deuterosomal cells
)
# p1 <- DotPlot(sce.all.int, features = unique(str_to_upper(selected_genes)),
#               assay='RNA' ,group.by = 'RNA_snn_res.0.8'  )  + coord_flip() # +th
# 
# p1

# gplots::balloonplot(table(sce.all.int$RNA_snn_res.0.8,sce.all.int$orig.ident))

# 还是以0.1分辨率为例

# > table(sce.all.int$RNA_snn_res.0.1)
# 
# 0     1     2     3     4     5     6     7     8     9    10 
# 15770 11663  6554  4366  2430   925   536   517   429   409   353 

gplots::balloonplot(table(sce.all.int$RNA_snn_res.0.1,sce.all.int$orig.ident))

输入一个列联表

自定义细胞亚群并根据定义分组再次进行check

# 前面的出图需要仔细看
# 理论上可以把细胞周期回归掉
# 也可以把文库大小回归掉

# 付费环节 800 元人民币
# 纯手工操作，非常耗时耗力
# 这个时候可以自由选择不同的分辨率
# 一般来说，为了偷懒就选择0.1的分辨率，但是不建议
# 如果是自己的项目，建议起码选择0.8分辨率以上。
# 会有不同的结果，但是影响不大。

if(F){
  sce.all.int
  celltype=data.frame(ClusterID=0:10 ,
                      celltype= 0:10) 
  #定义细胞亚群  
  celltype[celltype$ClusterID %in% c( 0,5 ),2]='Bcells' 
  celltype[celltype$ClusterID %in% c( 4 ),2]='plasma'   
  celltype[celltype$ClusterID %in% c( 6,8,10 ),2]='myeloids'   
  celltype[celltype$ClusterID %in% c( 1,9 ),2]='cd4_conv_Tcells' 
  celltype[celltype$ClusterID %in% c( 2 ),2]='cd8_Tcells'   
  celltype[celltype$ClusterID %in% c( 3 ),2]='cd4_Treg_Tcells' 
  celltype[celltype$ClusterID %in% c( 7 ),2]='NK'   
  
  head(celltype)
  celltype
  table(celltype$celltype)
 #        Bcells cd4_conv_Tcells cd4_Treg_Tcells      cd8_Tcells        myeloids 
 #             2               2               1               1               3 
 #            NK          plasma 
 #             1               1 
  sce.all.int@meta.data$celltype = "NA"
  
  for(i in 1:nrow(celltype)){
    sce.all.int@meta.data[which(sce.all.int@meta.data$RNA_snn_res.0.1 == celltype$ClusterID[i]),'celltype'] <- celltype$celltype[i]}
  Idents(sce.all.int)=sce.all.int$celltype
  
  sel.clust = "celltype"
  sce.all.int <- SetIdent(sce.all.int, value = sel.clust)
  table(sce.all.int@active.ident) 
#         Bcells      cd8_Tcells cd4_conv_Tcells cd4_Treg_Tcells          plasma 
#          16695            6554           12072            4366            2430 
#             NK        myeloids 
#            517            1318 
  
  dir.create('check-by-celltype')
  setwd('check-by-celltype')
  source('../scRNA_scripts/check-all-markers.R')
  setwd('../') 
  getwd()
}

phe=sce.all.int@meta.data
save(phe,file = 'phe.Rdata')

还是以gastric_cancer_markers、GSE150580_breast_cancer_markers_list这两个基因集为例

发现脚本画的gastric_cancer_markers dotplot导出为pdf挤到了一起

这里单独画一下这个基因集

需要注意的是

> DotPlot(sce.all.int , features = str_to_upper(gastric_cancer_markers) )  + 
+     coord_flip() + 
+     theme(axis.text.x=element_text(angle=45,hjust = 1))
Error in FetchData.Seurat(object = object, vars = features, cells = cells) : 
  None of the requested variables were found (10 out of 18 shown): PTPRC, MUC2, ITLN1, FABP1, APOA1, CEACAM5, CEACAM6, EPCAM, KRT18, MUC1
> 
> DotPlot(sce.all.int , features = str_to_title(gastric_cancer_markers) )  + 
+     coord_flip() + 
+     theme(axis.text.x=element_text(angle=45,hjust = 1))
Warning message:
In FetchData.Seurat(object = object, vars = features, cells = cells) :
  The following requested variables were not found: Muc2, Ceacam5, Ceacam6, Tff2, Pga4, Pga3, Muc5ac, Chgb

基因名格式应该保持一致