單細胞分析實錄(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=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
需要提前準備好兩個資料夾，比如我用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/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}

最終的表達矩陣會輸出到

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