單細胞分析實錄(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序列就包含在文件的第一、五、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

tag一、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