Cell Ranger是一個(gè)“傻瓜”軟件很洋，你只需提供原始的fastq文件项贺，它就會(huì)返回feature-barcode表達(dá)矩陣勇边。為啥不說(shuō)是gene-cell褂策，舉個(gè)例子横腿，cell hashing數(shù)據(jù)得到的矩陣還有tag行，而列也不能肯定就是一個(gè)cell辙培，可能考慮到這個(gè)才不叫g(shù)ene-cell矩陣吧~它是10xgenomics提供的官方比對(duì)定量軟件，有四個(gè)子命令邢锯，我只用過(guò)cellranger count扬蕊，另外三個(gè)cellranger mkfastq、cellranger aggr丹擎、cellranger reanalyze沒用過(guò)尾抑，也沒啥影響。
下載：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的使用之前蒂培，先來(lái)看一下10X的單細(xì)胞數(shù)據(jù)長(zhǎng)什么樣

這是一個(gè)樣本5個(gè)Line的測(cè)序數(shù)據(jù)再愈，數(shù)據(jù)量足夠的話可能只有一個(gè)Line』ご粒可以看出翎冲，它們的命名格式相對(duì)規(guī)范，在收到公司的數(shù)據(jù)后媳荒，盡量不要自己更改命名抗悍。此外還要注意一個(gè)細(xì)節(jié)，就是存放這些fastq文件的目錄應(yīng)該用第一個(gè)下劃線_前面的字符串命名钳枕，否則后續(xù)cell ranger將無(wú)法識(shí)別目錄里面的文件缴渊，同時(shí)報(bào)錯(cuò)

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

其實(shí)并不是--chemistry參數(shù)的問題。
為了更清楚地理解文件內(nèi)容鱼炒，我們來(lái)看一下10X單細(xì)胞的測(cè)序示意圖

Read1那一段序列原本是連在磁珠上面的衔沼，有cellular barcode（一個(gè)磁珠上都一樣），有UMI（各不相同），還有poly-T指蚁。Read2就是來(lái)源于細(xì)胞內(nèi)的RNA菩佑。它倆連上互補(bǔ)配對(duì)之后，還會(huì)在Read2的另一端連上sample index序列欣舵。這段sample index序列的作用是什么呢擎鸠？可以參考illumina測(cè)序中index primers的作用：

簡(jiǎn)單來(lái)說(shuō)就是為了在一次測(cè)序中，測(cè)多個(gè)樣本缘圈，在來(lái)源于特定樣本的序列后都加上特定的index劣光，測(cè)完之后根據(jù)對(duì)應(yīng)關(guān)系拆分。一個(gè)樣本對(duì)應(yīng)4個(gè)index：

再看每個(gè)文件里面是什么就容易理解了糟把，我們以一個(gè)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

其實(shí)這個(gè)index序列就包含在文件的第1绢涡、5、9...行遣疯，有點(diǎn)多余雄可，一般不太關(guān)注它。這個(gè)文件的序列最多四種缠犀，感興趣的小伙伴可以看看数苫。

R1文件里面就是cellular barcode信息，多余的序列已經(jīng)去掉了辨液。10X的v2試劑堿基長(zhǎng)度是26虐急，v3試劑堿基長(zhǎng)度是28

最后一個(gè)文件就是真正的轉(zhuǎn)錄本對(duì)應(yīng)的cDNA序列
上一篇講到cell hashing測(cè)序有轉(zhuǎn)錄本信息，得到的文件和上面是一樣的滔迈；還有一個(gè)細(xì)胞表面蛋白信息止吁，根據(jù)這個(gè)蛋白信息區(qū)分細(xì)胞來(lái)源，如下：

從圖中可以看出燎悍，和普通轉(zhuǎn)錄本建庫(kù)差不多敬惦，就是R2那一部分換成了HTO序列，整個(gè)片段長(zhǎng)度也改變了谈山。

上面兩張圖是我在實(shí)際處理中看到的兩種cell hashing測(cè)序俄删，第一張是TotalSeqA，第二張是TotalSeqB奏路。TotalSeqA中抗蠢，R2第一個(gè)堿基開始為HTO序列（之后是polyA序列），而TotalSeqB中思劳，R2前10個(gè)堿基為N的任意堿基迅矛，第11個(gè)堿基為HTO序列的開始位置，HTO序列長(zhǎng)度為16潜叛。

綜上秽褒，cell hashing的測(cè)序數(shù)據(jù)有兩套壶硅，一套是常規(guī)的轉(zhuǎn)錄本fastq，一套是蛋白信息（也可以說(shuō)是樣本信息）的fastq销斟。所以處理這類數(shù)據(jù)庐椒，要跟測(cè)序公司確認(rèn)清楚用的是TotalSeqA還是B，以及樣本和HTO序列的對(duì)應(yīng)關(guān)系蚂踊。

接下來(lái)說(shuō)說(shuō)如何用Cell Ranger處理普通10X單細(xì)胞測(cè)序數(shù)據(jù)约谈，以及cell hashing單細(xì)胞測(cè)序數(shù)據(jù)
普通10X

indir=/project_2019_11/data/S20191015T1
outdir=/project_2019_11/cellranger/
sample=S20191015T1
ncells=5000 #預(yù)計(jì)細(xì)胞數(shù)，這個(gè)參數(shù)對(duì)最終能得到的細(xì)胞數(shù)影響并不大犁钟，所以不用糾結(jié)
threads=20

refpath=/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
需要提前準(zhǔn)備好兩個(gè)文件夾棱诱，比如我用total_seq_A或total_seq_B存放HTO序列和樣本來(lái)源的對(duì)應(yīng)關(guān)系：

$ 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對(duì)應(yīng)哪一個(gè)樣本事先知道涝动；^(BC)可以看做正則表達(dá)式迈勋，表示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個(gè)堿基之后就是HTO序列醋粟，后面的序列隨意
lib_csv
第二個(gè)文件夾lib_csv靡菇，用來(lái)存放cell hashing兩套數(shù)據(jù)的路徑，用csv格式存儲(chǔ)米愿，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/cellranger

cd ${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}

最終的表達(dá)矩陣會(huì)輸出到

${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存儲(chǔ)的是基因信息厦凤，因?yàn)槭莄ell hashing數(shù)據(jù)，矩陣最后多了幾行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為矩陣信息，除前三行外宙搬，余下的行數(shù)等于feature乘以CB數(shù)笨腥，第二列表示CB編號(hào)拓哺，從1到10139勇垛，1重復(fù)33540次，對(duì)應(yīng)第一列的33540個(gè)feature士鸥。第三列表示UMI
下面的腳本可以將這三個(gè)文件轉(zhuǎn)換為常見的矩陣形式

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