一文搞定基本cellranger定量
引言
在上次推文中，我们正式开启了“初学者暑期搞定单细胞”这个专辑，并学习了如何下载单细胞测序上游fastq文件，这其中我们着重强调了scRNAseq测序原理的重要性，这会导致我们输入输出文件的不同，影响接下来的分析，这次我们就来系统地学习scRNAseq测序以及10X技术配套软件cellranger的使用。
参考原推文：
scRNAseq测序技术基础知识：
测序流程：现在主流的主要10X Genomics Chromium（较多细胞），SAMRT-seq2（较多基因）和Fluidigm C1等。当然还有其他的如：CELL-seq、Drop-seq、mas-seq和Wafergen ICELL8等。
最上游的实验阶段，我们简单地理解成两部分：单细胞分离和单细胞定量
- 单细胞分离 （Isolation of Single Cells）
单细胞分离方法：包括有限稀释法（Limiting dilution technique）显微操作法（micromanipulation）、流式细胞分选(FACS)、激光捕获显微分离(LCM)、微孔（microwell）、微流控（microfluidics）和微滴（droplet）。
微滴技术 是将单个细胞包裹在µl级别的液滴中，液滴被搭载到建库所用的酶上，每个微滴包含一个独特的条码(barcode)，由那个被包装好的细胞产生的所有reads都被贴上了该条码，也是为了之后对于不同细胞reads的分辨
- 单细胞定量
包括两种类型：全长以及基于标签(tag)。前者对每个转录本都试图获得一致的read覆盖度，后者只捕获5‘或者3’端的RNA。定量方法的选择也影响了后续分析的方法选择。
理论上，全长的方法应该得到转录本的平均覆盖度，但是实际上，覆盖度经常是有偏差的。基于标签的方法能够利用特异性分子标记(Unique Molecular Identifiers, UMIs)提高定量准确度，但是呢，这种限制了转录组一端的方法有降低了转录本的可拼接性，让以后的isoform识别变得困难。
详细看看10X技术：
可以发现和bulk测序中illumina双端很像，其实本质上是一样的，只不过其中一端（R1）变成了barcode+UMI
我们获取scRNAseq的fastq文件至少是两个，除了barcode+UMI和reads外，有的还有index
index和barcode有什么区别，为什么用两个fq文件进行区分？ i7 sample index是加到Illumina测序接头上的，保证多个测序文库可以在同一个flow-cell上或者同一个lane上进行混合测序(multiplexed)。当然可以自己指定index，但更多情况下会使用10X公司提供的index序列(bundled index sets)，针对不同项目使用的index也是不同的。不过共性就是：96孔板的每个孔中都加入了4种不同的index oligos混合(详见：https://kb.10xgenomics.com/hc/en-us/articles/218168503-What-oligos-are-in-my-sample-index-)。 它的作用就是在CellRanger的
mkfastq功能中体现出来的，它自动识别样本index名称（例如：SA-GA-A1），将具有相同4种oligo的fq文件组合在一起，表示同一个样本
我们使用barcode来鉴别细胞，UMI来控制PCR偏差和鉴别基因（转录本）
我根据我个人的理解，绘制了一个草图帮助理解”UMI来控制PCR偏差“
每个UMI对应一个mRNA，mRNA可以是不同\相同基因的、不同\相同转录本的，但都对应不同UMI（一个barcode上每种UMI都是唯一的）
cellranger基本流程：
修改文件名称
在上一期 今年暑假一起学单细胞吧（附上游数据下载tips） 我们获得了sra通过--split-files得到的3个fq文件，其实这后面还需要修改一下名字
R1就是barcode+UMI序列
使用原推文代码 单细胞实战(三) Cell Ranger使用初探
# 以P2586-4为例
mkdir -p $wkd/qc
cd $wkd/qc
find $wkd/raw/P2586-4 -name '*R1*.gz'>P2586-4-id-1.txt
find $wkd/raw/P2586-4 -name '*R2*.gz'>P2586-4-id-2.txt
cat P2586-4-id-1.txt P2586-4-id-2.txt >P2586-4-id-all.txt
cat P2586-4-id-all.txt| xargs fastqc -t 20 -o ./
根据自己当前工作目录进行修改
fastqc质控结果解读：
https://www.bioinformatics.babraham.ac.uk/projects/fastqc/Help/3%20Analysis%20Modules/
背景补充
原推文单细胞实战(三) Cell Ranger使用初探 介绍了许多10X多种不同的测序情况，并且介绍了如何用cellranger来处理这些不同的情况
主要根据sample、library、flowcell的数量来定义分析的复杂程度(由浅入深)
先学最简单的，一个sample 一个library 一个flowcell
目前应该是最简单的，因为我看sample和lane都是唯一的，S1_L001
原推文的cellranger版本已经较老了，这里我去官网下载的最新版：
https://support.10xgenomics.com/single-cell-gene-expression/software/downloads/latest?
参考基因组：
cellranger定量
在原推文中，介绍了一种cellranger mkfastq拆分BCLs（每个flowcell 的Illumina sequencer's base call files）数据得到fastq文件的方法
我们这里直接使用前面sra拆分得到的fastq文件
联系前面提到
许多10X多种不同的测序情况，并且介绍了如何用cellranger来处理这些不同的情况 主要根据sample、library、flowcell的数量来定义分析的复杂程度(由浅入深)
原推文提到，这些不同情况也有不同的fq文件位置需要注意，这里我们就不深入探究了
原推文v2版本代码：
我使用最新版定量代码：
ref=../soft/refdata-gex-mm10-2020-A/
cellranger=../soft/cellranger-7.1.0/cellranger
ls ../raw/*gz|cut -d"_" -f 1 |sort -u|cut -d"/" -f 3 | cut -d "_" -f 1 | uniq | while read id;do
nohup $cellranger count --id=$id \
--transcriptome=$ref \
--fastqs=../raw \
--sample=$id \
--nosecondary \
--localcores=4 \
--localmem=30 &
done
使用资源限制：Cell Ranger默认在本地运行(或者使用
--jobmode=local指定)，它会占用90%的空余内存以及所有空余的CPU。如果要进行资源限制，可以使用—localmem或者--localcores如果使用的是共享服务器就需要注意这个问题
输出文件比较：
cellranger7.1.0输出文件目录：
可以发现默认输出文件格式基本一致
从上到下依次来看：
- web_summary.html：官方说明 summary HTML file
- metrics_summary.csv：CSV格式数据摘要
- possorted_genome_bam.bam：比对文件
- possorted_genome_bam.bam.bai：索引文件
- filtered_gene_bc_matrices：是重要的一个目录 ，下面又包含了 barcodes.tsv.gz、features.tsv.gz、matrix.mtx.gz，是下游Seurat、Scater、Monocle等分析的输入文件
- filtered_feature_bc_matrix.h5：过滤掉的barcode信息HDF5 format
- raw_feature_bc_matrix：原始barcode信息
- raw_feature_bc_matrix.h5：原始barcode信息HDF5 format
- analysis：数据分析目录，下面又包含聚类clustering（有graph-based & k-means）、差异分析diffexp、主成分线性降维分析pca、非线性降维tsne
- molecule_info.h5：下面进行aggregate使用的文件
- cloupe.cloupe：官方可视化工具Loupe Cell Browser 输入文件
此外，原推文还提到了一些内置软件和算法、如何自主构建参考信息以及多个文库的整合 aggr
很多时候，我们需要根据自己的需要，自定义一套参考信息
当处理多个生物学样本或者一个样本存在多个重复/文库时，最好的操作就是先分别对每个文库进行单独的count定量，然后将定量结果利用
aggr组合起来