白癜风单细胞转录组数据处理

生信技能树

发布于 2023-02-27 21:32:19

3590

发布于 2023-02-27 21:32:19

前些天我们公众号弄了一个活动，详见：春节期间单细胞转录组数据分析全免费，收到了上百个需求，本来呢我们自己就算是春节前后14天不吃不喝不眠不休也不可能完成这么多单细胞数据处理。好在我灵机一动，想起来了前面两个月培养的一百多个在线实习生，毕竟教了大家R语言，转录组，以及单细胞转录组。

所以我写了一个还算是比较自动化的单细胞转录组数据处理代码，如果是我自己来做这样的处理，可以在十几分钟就完成复现文章的第一层次降维聚类分群图，如下所示的2分组的15个样品：

降维聚类分群图

首先下载文件，在https://ngdc.cncb.ac.cn/omix/release/OMIX691 可以下载：

PRJCA006797_RAW % ls -lh *gz|cut -d" " -f 6-

   12M  1 16 09:59 OMIX691-20-01.tar.gz
   25M  1 16 09:59 OMIX691-20-02.tar.gz
   23M  1 16 09:59 OMIX691-20-03.tar.gz
   42M  1 16 10:00 OMIX691-20-04.tar.gz
   36M  1 16 10:00 OMIX691-20-05.tar.gz
  7.9M  1 16 10:00 OMIX691-20-06.tar.gz
   10M  1 16 10:00 OMIX691-20-07.tar.gz
   27M  1 16 10:00 OMIX691-20-08.tar.gz
   11M  1 16 10:00 OMIX691-20-09.tar.gz
   20M  1 16 10:00 OMIX691-20-10.tar.gz
   23M  1 16 10:00 OMIX691-20-11.tar.gz
  7.3M  1 16 10:00 OMIX691-20-12.tar.gz
   31M  1 16 10:01 OMIX691-20-13.tar.gz
   29M  1 16 10:01 OMIX691-20-14.tar.gz
   19M  1 16 10:01 OMIX691-20-15.tar.gz

每个样品是一个压缩包，需要解压：

ls *gz |while read id;do (tar zxvf $id);done

可以看到每个文件夹里面的都是一个标准的10x输出：


P09/outs/filtered_gene_bc_matrices/
P09/outs/filtered_gene_bc_matrices/GRCh38/
P09/outs/filtered_gene_bc_matrices/GRCh38/matrix.mtx
P09/outs/filtered_gene_bc_matrices/GRCh38/barcodes.tsv
P09/outs/filtered_gene_bc_matrices/GRCh38/genes.tsv
P10/outs/filtered_gene_bc_matrices/
P10/outs/filtered_gene_bc_matrices/GRCh38/
P10/outs/filtered_gene_bc_matrices/GRCh38/matrix.mtx
P10/outs/filtered_gene_bc_matrices/GRCh38/barcodes.tsv
P10/outs/filtered_gene_bc_matrices/GRCh38/genes.tsv

接下来就是读取它，然后降维聚类分群啦。值得注意的是这个时候每个文件夹里面的文件都是没有压缩状态，所以耗费磁盘空间比较大：

├── [  96]  P09
│   └── [  96]  outs
│       └── [  96]  filtered_gene_bc_matrices
│           └── [ 160]  GRCh38
│               ├── [ 77K]  barcodes.tsv
│               ├── [821K]  genes.tsv
│               └── [111M]  matrix.mtx
└── [  96]  P10
    └── [  96]  outs
        └── [  96]  filtered_gene_bc_matrices
            └── [ 160]  GRCh38
                ├── [ 55K]  barcodes.tsv
                ├── [821K]  genes.tsv
                └── [ 74M]  matrix.mtx

一般来说，我们会压缩它，但是压缩它就意味着需要把 genes.tsv 修改为 features.tsv 的格式，所以这里略。很轻松的把这么多样品的表达量矩阵读入，并且构建成为了Seurat对象：

rm(list=ls())
options(stringsAsFactors = F) 
source('scRNA_scripts/lib.R')

###### step1:导入数据 ######   
# 付费环节 800 元人民币
dir='PRJCA006797_RAW/' 
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,'outs/filtered_gene_bc_matrices/GRCh38')) ,
                          project =  gsub('_filtered_feature_bc_matrix','',pro)  ,
                          min.cells = 5,
                          min.features = 500 )
  return(sce)
})
names(sceList) 
# gsub('^GSM[0-9]*','',samples)
sce.all=merge(x=sceList[[1]],
              y=sceList[ -1 ],
              add.cell.ids =gsub('_filtered_feature_bc_matrix','',samples)     )

as.data.frame(sce.all@assays$RNA@counts[1:10, 1:2])
head(sce.all@meta.data, 10)
table(sce.all$orig.ident) 
# 分组信息，通常是隐含在文件名，样品名字里面
# phe=str_split( colnames(sce.all),'[-_]',simplify = T)
# head(phe)
# table(phe[,1])
sce.all$group= substring(colnames(sce.all),1,1)
# sce.all@meta.data$group =  ifelse(grepl('ne',sce.all$orig.ident),'ne','e')
table(sce.all@meta.data$group)
table(sce.all@meta.data$orig.ident)

上面的 scRNA_scripts 文件夹里面的 lib.R 就是加载一些单细胞转录组数据分析常见的包即可。

然后整合多个样品并且降维聚类分群

下面的scRNA_scripts 文件夹里面的 qc.R ，就是检测线粒体基因，核糖体基因等简单的质量控制：

###### 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))
print(dim(sce.all.filt))
setwd('../')

###### 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('../')

###### 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) 
sp='human'
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

可以看到，我给大家准备的基因列表主要是针对肿瘤单细胞的第一层次降维聚类分群，是：

immune (CD45+,PTPRC),
epithelial/cancer (EpCAM+,EPCAM),
stromal (CD10+,MME,fibo or CD31+,PECAM1,endo)

参考我前面介绍过 CNS图表复现08—肿瘤单细胞数据第一次分群通用规则，这3大单细胞亚群构成了肿瘤免疫微环境的复杂。绝大部分文章都是抓住免疫细胞亚群进行细分，包括淋巴系（T,B,NK细胞）和髓系（单核，树突，巨噬，粒细胞）的两大类作为第二次细分亚群。但是也有不少文章是抓住stromal 里面的fibo 和endo进行细分，并且编造生物学故事的。

如果仅仅是使用我们的前面的 scRNA_scripts/check-all-markers.R ，会发现很多细胞亚群其实没办法给名字。而这个文章里面的各个单细胞亚群给出来了基因列表，所以需要我们自己肉眼看看这些基因：