使用docker分析单细胞数据库对笔记本电脑的要求

今天来介绍一下GEO单细胞转录组下载数据以及整理，单细胞测序的原理以及数据结果都与bulk测序的方式有一定的差距，所以我们单独说一下。

桓峰基因的教程不但教您怎么使用，还会定期分析一些相关的文章，学会教程只是基础，但是如果把分析结果整合到文章里面才是目的，觉得我们这些教程还不错，并且您按照我们的教程分析出来不错的结果发了文章记得告知我们，并在文章中感谢一下我们哦！

公司英文名称：Kyoho Gene Technology (Beijing) Co.,Ltd.

桓峰基因公众号推出单细胞系列教程，有需要生信的老师可以联系我们！首选看下转录分析教程整理如下：

SCS【1】今天开启单细胞之旅，述说单细胞测序的前世今生

SCS【2】单细胞转录组 之 cellranger

SCS【3】单细胞转录组数据 GEO 下载及读取

---

前言

单细胞近几年也是越来越火热，因此数据也都上传到GEO，这样我们有些分析就可以不用测序，而是通过单细胞现有的数据进行分析了，目前GEO提供数据大致包括三类：

1. 原始数据：fastq格式原始测序数据
2. 经过cellranger标准化后的10X数据
3. Counts数据信息

当然我们都希望获得第二种数据，毕竟从原始下机数据fastq处理起来，自己的小电脑也受不起，并且不是每个人都很擅长 Linux 和 shell ，所以有了第二种，对于我们来说算是省去了很多繁琐的步骤。

数据下载

确定我们需要下载的GEO编号，我们这里使用 GSE135927， 登录网站

https://www.ncbi.nlm.nih.gov/geo/

数据集介绍

我们这里选择Cell Ranger生成的raw count。

下载GSE135927_RAW.tar后解压，我们可以看到每个样本有三个文件，包括：

1. barcodes.tsv

Barcodes通俗来讲就是每个细胞的代码，组成就是ATCG四个碱基排列组合成的不同的14个碱基组合；

AAACCTGAGGGCTCTC-1
AAACCTGAGTGAACGC-1
AAACCTGGTATGAATG-1
AAACCTGTCTGCGACG-1
AAACGGGAGGTGCAAC-1

2. features.tsv

features.tsv一般是基因的ensembl ID 和symbol

ENSG00000243485	MIR1302-2HG	Gene Expression
ENSG00000237613	FAM138A	Gene Expression
ENSG00000186092	OR4F5	Gene Expression
ENSG00000238009	AL627309.1	Gene Expression
ENSG00000239945	AL627309.3	Gene Expression
ENSG00000239906	AL627309.2	Gene Expression
ENSG00000241599	AL627309.4	Gene Expression

3. matrix.mtx

matrix.mtx就是每个细胞不同基因的表达矩阵：

%%MatrixMarket matrix coordinate integer general
%metadata_json: {"format_version": 2, "software_version": "3.0.1"}
3353813593017342
33509124
3350812
3350716
3350616
33505126
33504119
33503119
33502124
33501113

软件安装

我们使用Seurat软件包里的函数Read10X()进行数据读取。

if(!require(Seurat))
  BiocManager::install("Seurat")

数据读取

我们发现下载后的每个样本包含三个文件，准备好每个样本一个文件夹。我们先看下函数说明给出来的例子：

1. CellRanger < 3.0

如果是CellRanger < 3.0 单样本读取的时候必须是每个文件格式为 barcodes.tsv, genes.tsv, and matrix.mtx

## Not run: 
# For output from CellRanger < 3.0
data_dir <- 'path/to/data/directory'
list.files(data_dir) # Should show barcodes.tsv, genes.tsv, and matrix.mtx
expression_matrix <- Read10X(data.dir = data_dir)
seurat_object = CreateSeuratObject(counts = expression_matrix)

2. CellRanger >= 3.0

如果是CellRanger >= 3.0 单样本读取的时候必须是每个文件格式为barcodes.tsv.gz, features.tsv.gz, and matrix.mtx.gz

# For output from CellRanger >= 3.0 with multiple data types
data_dir <- 'path/to/data/directory'
list.files(data_dir) # Should show barcodes.tsv.gz, features.tsv.gz, and matrix.mtx.gz
data <- Read10X(data.dir = data_dir)
seurat_object = CreateSeuratObject(counts = data$`Gene Expression`)
seurat_object[['Protein']] = CreateAssayObject(counts = data$`Antibody Capture`)

3. 实际例子

我们先看下单个样本处理方法和多个样本合并。

单个样本读取

library(Seurat)
data_dir <- "./GSE135927"
dirs = list.dirs(data_dir)
dirs_sample = dirs[2:3]
dirs_sample
## [1] "./GSE135927/GSM4038043" "./GSE135927/GSM4038044"
GSM4038043 = Read10X(data.dir = dirs_sample[1])
GSM4038043[1:10, 1:3]
## 10 x 3 sparse Matrix of class "dgCMatrix"
##             AAACCTGAGGGCTCTC-1 AAACCTGAGTGAACGC-1 AAACCTGGTATGAATG-1
## MIR1302-2HG                  .                  .                  .
## FAM138A                      .                  .                  .
## OR4F5                        .                  .                  .
## AL627309.1                   .                  .                  .
## AL627309.3                   .                  .                  .
## AL627309.2                   .                  .                  .
## AL627309.4                   .                  .                  .
## AL732372.1                   .                  .                  .
## OR4F29                       .                  .                  .
## AC114498.1                   .                  .                  .
s1 = CreateSeuratObject(counts = GSM4038043, min.cells = 3, min.features = 200)
head([e-mail])
##                       orig.ident nCount_RNA nFeature_RNA
## AAACCTGAGGGCTCTC-1 SeuratProject       8506         2450
## AAACCTGAGTGAACGC-1 SeuratProject       5634         1743
## AAACCTGGTATGAATG-1 SeuratProject       8841         2309
## AAACCTGTCTGCGACG-1 SeuratProject       4484         1779
## AAACGGGAGGTGCAAC-1 SeuratProject       8221         2332
## AAACGGGTCGCTAGCG-1 SeuratProject       6471         2153

多个样本读取

GSM4038044 = Read10X(data.dir = dirs_sample[2])
s2 = CreateSeuratObject(counts = GSM4038044, min.cells = 3, min.features = 200)
m = merge(x = s1, y = s2, add.cell.ids = c("GSM4038043", "GSM4038044"), merge.data = TRUE)
as.data.frame(m@assays$RNA@counts[1:10, 1:2])
##            GSM4038043_AAACCTGAGGGCTCTC-1 GSM4038043_AAACCTGAGTGAACGC-1
## AL669831.5                             0                             0
## LINC00115                              0                             0
## FAM41C                                 0                             0
## NOC2L                                  0                             0
## KLHL17                                 0                             0
## PLEKHN1                                0                             0
## HES4                                   0                             0
## ISG15                                  2                             0
## AL645608.2                             0                             0
## AGRN                                   0                             0
head([e-mail])
##                                  orig.ident nCount_RNA nFeature_RNA
## GSM4038043_AAACCTGAGGGCTCTC-1 SeuratProject       8506         2450
## GSM4038043_AAACCTGAGTGAACGC-1 SeuratProject       5634         1743
## GSM4038043_AAACCTGGTATGAATG-1 SeuratProject       8841         2309
## GSM4038043_AAACCTGTCTGCGACG-1 SeuratProject       4484         1779
## GSM4038043_AAACGGGAGGTGCAAC-1 SeuratProject       8221         2332
## GSM4038043_AAACGGGTCGCTAGCG-1 SeuratProject       6471         2153

总结

另外还有一种，只有一个表格，下载之后可以自己整理结果，例如 GSE135928 下载之后我们看下文件格式：

barcode	is_cell	contig_id	high_confidence	length	chain	v_gene	d_gene	j_gene	c_gene	full_length	productive	cdr3	cdr3_nt	reads	umis	raw_clonotype_id	raw_consensus_id
AAACGGGTCGTCGTTC-1TRUE	AAACGGGTCGTCGTTC-1_contig_1	TRUE582	TRA	TRAV5	None	TRAJ8	TRAC	TRUETRUE	CAESSGTGFQKLVF	TGTGCAGAGAGTTCCGGCACAGGCTTTCAGAAACTTGTATTT	19453	clonotype17	clonotype17_consensus_2
AAACGGGTCGTCGTTC-1TRUE	AAACGGGTCGTCGTTC-1_contig_2	TRUE710	TRB	TRBV4-3	TRBD2	TRBJ2-7	TRBC2	TRUETRUE	CASSQDSGPSYEQYF	TGCGCCAGCAGCCAAGATAGTGGGCCATCCTACGAGCAGTACTTC	468810	clonotype17	clonotype17_consensus_1
AAACGGGTCTCGAGTA-1TRUE	AAACGGGTCTCGAGTA-1_contig_1	TRUE689	TRB	TRBV5-1	TRBD1	TRBJ2-2	TRBC2	TRUETRUE	CASSYRGVNTGELFF	TGCGCCAGCAGCTACAGGGGCGTGAACACCGGGGAGCTGTTTTTT	619013	clonotype18	clonotype18_consensus_2
AAACGGGTCTCGAGTA-1TRUE	AAACGGGTCTCGAGTA-1_contig_3	TRUE343	Multi	None	None	TRBJ2-3	TRBC2	FALSE	None	None	None	7374	clonotype18	None
AAACGGGTCTCGAGTA-1TRUE	AAACGGGTCTCGAGTA-1_contig_4	TRUE694	TRA	TRAV8-6	None	TRAJ52	TRAC	TRUETRUE	CAVIKANAGGTSYGKLTF	TGTGCTGTGATCAAAGCTAATGCTGGTGGTACTAGCTATGGAAAGCTGACATTT	54838	clonotype18	clonotype18_consensus_1
AAAGATGAGGTAAACT-1TRUE	AAAGATGAGGTAAACT-1_contig_1	TRUE605	TRA	TRAV12-3	None	TRAJ22	TRAC	TRUETRUE	CARGSARQLTF	TGTGCACGGGGTTCTGCAAGGCAACTGACCTTT	77239	clonotype19	clonotype19_consensus_1
AAAGATGAGGTAAACT-1TRUE	AAAGATGAGGTAAACT-1_contig_3	TRUE543	TRB	TRBV6-3	TRBD1	TRBJ2-5	TRBC2	TRUETRUE	CASSYSFSGQGGETQYF	TGTGCCAGCAGTTACTCATTTTCGGGACAAGGCGGAGAGACCCAGTACTTC	1098420	clonotype19	clonotype19_consensus_2
AAAGATGCACGACGAA-1TRUE	AAAGATGCACGACGAA-1_contig_1	TRUE517	TRB	TRBV4-1	TRBD2	TRBJ1-1	TRBC1	TRUETRUE	CASSQDGLMNTEAFF	TGCGCCAGCAGCCAAGATGGACTTATGAACACTGAAGCTTTCTTT	34097	clonotype20	clonotype20_consensus_2
AAAGATGCACGACGAA-1TRUE	AAAGATGCACGACGAA-1_contig_2	TRUE543	TRA	TRAV23DV6	None	TRAJ37	TRAC	TRUETRUE	CAASRVNTGKLIF	TGTGCAGCAAGCAGAGTCAACACAGGCAAACTAATCTTT	20283	clonotype20	clonotype20_consensus_1

这种格式起始也非常好用，利用TPM计算原理，有length，reads，barcode，raw_consensus_id 等信息计算出来获得一个表达矩阵即可。

桓峰基因，铸造成功的您！

未来桓峰基因公众号将不间断的推出单细胞系列生信分析教程，

敬请期待！！

References:

1. Singh J, Chen ELY, Xing Y, Stefanski HE et al. Generation and function of progenitor T cells from StemRegenin-1-expanded CD34+ human hematopoietic progenitor cells. Blood Adv 2019 Oct 22;3(20):2934-2948. PMID: 31648315