单细胞分析实录(2): 使用Cell Ranger得到表达矩阵

Cell Ranger是一个“傻瓜”软件，你只需提供原始的fastq文件，它就会返回feature-barcode表达矩阵。为啥不说是gene-cell，举个例子，cell hashing数据得到的矩阵还有tag行，而列也不能肯定就是一个cell，可能考虑到这个才不叫gene-cell矩阵吧~它是10xgenomics提供的官方比对定量软件，有四个子命令，我只用过cellranger count，另外三个cellranger mkfastq、cellranger aggr、cellranger reanalyze没用过，也没啥影响。
下载：https://support.10xgenomics.com/single-cell-gene-expression/software/downloads/latest
安装：https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/installation

在讲Cell Ranger的使用之前，先来看一下10X的单细胞数据长什么样

这是一个样本5个Line的测序数据，数据量足够的话可能只有一个Line。可以看出，它们的命名格式相对规范，在收到公司的数据后，尽量不要自己更改命名。此外还要注意一个细节，就是存放这些fastq文件的目录应该用第一个下划线_前面的字符串命名，否则后续cell ranger将无法识别目录里面的文件，同时报错

[error] Unable to detect the chemistry for the following dataset.
Please validate it and/or specify the chemistry
via the --chemistry argument.

其实并不是–chemistry参数的问题。
为了更清楚地理解文件内容，我们来看一下10X单细胞的测序示意图

Read1那一段序列原本是连在磁珠上面的，有cellular barcode（一个磁珠上都一样），有UMI（各不相同），还有poly-T。Read2就是来源于细胞内的RNA。它俩连上互补配对之后，还会在Read2的另一端连上sample index序列。这段sample index序列的作用是什么呢？可以参考illumina测序中index primers的作用：

简单来说就是为了在一次测序中，测多个样本，在来源于特定样本的序列后都加上特定的index，测完之后根据对应关系拆分。一个样本对应4个index：

再看每个文件里面是什么就容易理解了，我们以一个Line为例：

less -S S20191015T1_S6_L001_I1_001.fastq.gz | head -n 8
less -S S20191015T1_S6_L001_R1_001.fastq.gz | head -n 8
less -S S20191015T1_S6_L001_R2_001.fastq.gz | head -n 8

其实这个index序列就包含在文件的第1、5、9…行，有点多余，一般不太关注它。这个文件的序列最多四种，感兴趣的小伙伴可以看看。

R1文件里面就是cellular barcode信息，多余的序列已经去掉了。10X的v2试剂碱基长度是26，v3试剂碱基长度是28

最后一个文件就是真正的转录本对应的cDNA序列
上一篇讲到cell hashing测序有转录本信息，得到的文件和上面是一样的；还有一个细胞表面蛋白信息，根据这个蛋白信息区分细胞来源，如下：

从图中可以看出，和普通转录本建库差不多，就是R2那一部分换成了HTO序列，整个片段长度也改变了。

上面两张图是我在实际处理中看到的两种cell hashing测序，第一张是TotalSeqA，第二张是TotalSeqB。TotalSeqA中，R2第一个碱基开始为HTO序列（之后是polyA序列），而TotalSeqB中，R2前10个碱基为N的任意碱基，第11个碱基为HTO序列的开始位置，HTO序列长度为16。

综上，cell hashing的测序数据有两套，一套是常规的转录本fastq，一套是蛋白信息（也可以说是样本信息）的fastq。所以处理这类数据，要跟测序公司确认清楚用的是TotalSeqA还是B，以及样本和HTO序列的对应关系。

接下来说说如何用Cell Ranger处理普通10X单细胞测序数据，以及cell hashing单细胞测序数据
普通10X

indir=/project_2019_11/data/S20191015T1
outdir=/project_2019_11/cellranger/
sample=S20191015T1
ncells=5000 #预计细胞数，这个参数对最终能得到的细胞数影响并不大，所以不用纠结
threads=20refpath=/ref/10x/human/refdata-cellranger-GRCh38-3.0.0
cellranger=/softwore/bin/cellranger
cd ${outdir}
${cellranger} count --id=${sample} \--transcriptome=${refpath} \--fastqs=${indir} \--sample=${sample} \--expect-cells=${ncells} \--localcores=${threads}

cell hashing

total_seq_A
需要提前准备好两个文件夹，比如我用total_seq_A或total_seq_B存放HTO序列和样本来源的对应关系：

$ ls
feature.reference1.csv
$ cat feature.reference1.csv
id,name,read,pattern,sequence,feature_type
tag1,tag1,R2,^(BC),GTCAACTCTTTAGCG,Antibody Capture
tag2,tag2,R2,^(BC),TGATGGCCTATTGGG,Antibody Capture

tag1、tag2对应哪一个样本事先知道；^(BC)可以看做正则表达式，表示R2序列以barcode(也就是HTO序列)开始
total_seq_B

$ ls
feature.reference.csv
$ cat feature.reference.csv
id,name,read,pattern,sequence,feature_type
tag6,tag6,R2,5PNNNNNNNNNN(BC)NNNNNNNNN,GGTTGCCAGATGTCA,Antibody Capture
tag7,tag7,R2,5PNNNNNNNNNN(BC)NNNNNNNNN,TGTCTTTCCTGCCAG,Antibody Capture

5PNNNNNNNNNN(BC)NNNNNNNNN表示从5端开始，10个碱基之后就是HTO序列，后面的序列随意
lib_csv
第二个文件夹lib_csv，用来存放cell hashing两套数据的路径，用csv格式存储，sample这一列为文件夹名称

$ cat S20200612P1320200702N.libraries.csv
fastqs,sample,library_type
/project_2019_11/data/fastq/,S20200612P1320200702N,Gene Expression
/project_2019_11/data/antibody_barcode/,S20200612P13F20200702N,Antibody Capture

最终脚本如下

lib_dir=/script/cellranger/1/lib_csv/
#need to be changed based on your seq-tech: total_seq_A or total_seq_B
feature_ref_dir=/script/cellranger/1/total_seq_A/
outdir=/project_2019_11/cellranger/
sample=S20191017P11
ncells=5000
threads=20
refpath=/ref/10x/human/refdata-cellranger-GRCh38-3.0.0
cellranger=/softwore/bin/cellrangercd ${outdir}
${cellranger} count --libraries=${lib_dir}${sample}.libraries.csv \--r1-length=28 \--feature-ref=${feature_ref_dir}feature.reference1.csv \--transcriptome=${refpath} \--localcores=${threads} \--expect-cells=${ncells} \--id=${sample}

最终的表达矩阵会输出到

${outdir}${sample_id}/outs/filtered_feature_bc_matrix

$ cd S20200619P11120200716NC/outs/filtered_feature_bc_matrix/
$ ls
barcodes.tsv.gz  features.tsv.gz  matrix.mtx.gz$ less -S features.tsv.gz
ENSG00000243485 MIR1302-2HG Gene Expression
ENSG00000237613 FAM138A Gene Expression
......
ENSG00000277475 AC213203.1  Gene Expression
ENSG00000268674 FAM231C Gene Expression
tag7    tag7    Antibody Capture
tag8    tag8    Antibody Capture

features.tsv.gz存储的是基因信息，因为是cell hashing数据，矩阵最后多了几行tag信息，共33540行

$ less -S barcodes.tsv.gz | head -n 4
AAACCCAAGACTTAAG-1
AAACCCAAGCTACTGT-1
AAACCCAAGGACTGGT-1
AAACCCAAGGCCTGCT-1

barcodes.tsv.gz存放的是最后得到的cellular barcode，共10139行

$ less -S matrix.mtx.gz | head -n 8
%%MatrixMarket matrix coordinate integer general
%metadata_json: {"format_version": 2, "software_version": "3.1.0"}
33540 10139 15746600
65 1 1
103 1 1
155 1 2
179 1 2
191 1 1

matrix.mtx.gz为矩阵信息，除前三行外，余下的行数等于feature乘以CB数，第二列表示CB编号，从1到10139，1重复33540次，对应第一列的33540个feature。第三列表示UMI
下面的脚本可以将这三个文件转换为常见的矩阵形式

path1=/softwore/biosoft/cellranger-3.1.0/cellranger
path2=/project_2019_11/cellranger/i=S20191211P71
${path1} mat2csv ${path2}${i}/outs/filtered_feature_bc_matrix ${path2}Feature_Barcode_Matrices/${i}.mat.count.csv
sed 's/,/\t/g' ${path2}Feature_Barcode_Matrices/${i}.mat.count.csv  > ${path2}Feature_Barcode_Matrices/${i}.mat.count.txt
sed -i 's/^\t//g' ${path2}Feature_Barcode_Matrices/${i}.mat.count.txt
rm -f ${path2}Feature_Barcode_Matrices/${i}.mat.count.csv

单细胞分析实录(2): 使用Cell Ranger得到表达矩阵相关推荐

单细胞分析实录(7): 差异表达分析/细胞类型注释
前面已经讲解了: 单细胞分析实录(1): 认识Cell Hashing 单细胞分析实录(2): 使用Cell Ranger得到表达矩阵单细胞分析实录(3): Cell Hashing数据拆分单细胞 ...
seurat提取表达矩阵_单细胞分析实录(5): Seurat标准流程
前面我们已经学习了单细胞转录组分析的:使用Cell Ranger得到表达矩阵和doublet检测,今天我们开始Seurat标准流程的学习.这一部分的内容,网上有很多帖子,基本上都是把Seurat官网P ...
单细胞分析实录(14): 细胞类型注释的另一种思路 — CellID
之前写过一篇细胞类型注释的推文:单细胞分析实录(7): 差异表达分析/细胞类型注释,里面介绍了细胞类型注释的两种思路.在谈到自动注释时,我们基本都会用SingleR这个软件,类似的软件还有很多,原理上 ...
单细胞分析实录(8): 展示marker基因的4种图形（一）
今天的内容讲讲单细胞文章中经常出现的展示细胞marker的图:tsne/umap图.热图.堆叠小提琴图.气泡图,每个图我都会用两种方法绘制. 使用的数据来自文献:Single-cell transcr ...
单细胞分析实录(15): 基于monocle2的拟时序分析
关于什么是"拟时序分析",可以参考本期推送的另一篇推文.这一篇直接演示代码 monocle2这个软件用得太多了,很多文章都是monocle2的图.因为只使用表达矩阵作为输入,相比于 ...
单细胞分析实录(9): 展示marker基因的4种图形（二）
在上一篇中,我已经讲解了展示marker基因的前两种图形,分别是tsne/umap图.热图,感兴趣的读者可以回顾一下.这一节我们继续学习堆叠小提琴图和气泡图. 3. 堆叠小提琴图展示marker基因 ...
初步分析CCLE和GDSC的数据——RNA表达矩阵
CCLE是Cancer Cell Line Encyclopedia的简称, 是与肿瘤有关的细胞系的百科全书项目. GDSC是Genomics of Drug Sensitivity in Cance ...
单细胞分析之质控（四）
学习目标知道如何导入和读取数据,并了解数据的质控,能够对数据进行质控和分析. 1. 质控准备在基因表达定量后,需要将这些数据导入到 R 中,以生成用于执行 QC(质控).下面将讨论定量数据的格式, ...
seurat提取表达矩阵_Hemberg-lab单细胞转录组数据分析
单细胞RNA-seq简介混合RNA-seq2000年末的重大技术突破,取代微阵列表达芯片被广泛使用通过混合大量细胞获取足够RNA用于建库测序,来定量每个基因的平均表达水平用于比较转录组,例如比较 ...

单细胞分析实录(2): 使用Cell Ranger得到表达矩阵

cell hashing

单细胞分析实录(2): 使用Cell Ranger得到表达矩阵相关推荐

最新文章

热门文章