单细胞转录组上游分析-上游分析流程及知识点整理

背景：从生信的角度学习单细胞转录组的上游相关流程以及知识点。

10x测序平台

从cellranger的学习开始，cellranger是处理10X测序平台单细胞转录组序列的软件。
首先理解样本sample，文库library，通道lane的概念：

1.样本sample：从单一来源（组织或者细胞）提取的细胞悬液
2.文库library：单个sample制备的10x 测序文库，运行一次chromium的单个芯片通道
3.通道lane:一台测序仪上的数据会在flowcell上产出，然后这些数据又会根据flowcell上的不同lane以及lane上不同样本的index进行区分

同一个sample可以制备成多个library进行测序，从而获得更多的细胞；一个flowcell上的不同样本根据样本的index或者不同的lane进行区分；
一个sample的一个library可以在多个flowcell上进行测序，再进行merge合并。

一个sample，一个library，一个flowcell：

image.png
一个sample，一个library，多个flowcell：
image.png
一个sample，多个library，多个flowcell：
image.png
多个sample，多个library，多个flowcell：
image.png

Cell Ranger主要的流程有：拆分数据 mkfastq、细胞定量 count、定量组合 aggr、调参reanalyze，还有一些小工具比如mkref、mkgtf、upload、sitecheck、mat2csv、vdj、mkvdjref、testrun。

一、mkfastq拆分数据：

image.png

将10x sequencing后的BCLs文件转换为fastq文件。

 cellranger mkfastq --id=bcl \
                    --run=/path/to/bcl \
                    --samplesheet=samplesheet-1.2.0.csv

 cellranger mkfastq --id=bcl \
                    --run=/path/to/bcl \
                    --csv=simple-1.2.0.csv
 #'其中id指定输出目录的名称，run指的是下机的原始BCL文件目录
 #'重要的就是测序lane、样本名称、index等信息

samplesheet-1.2.0.csv（测序相关信息）：

image.png

simple-1.2.0.csv（简化后的测序相关信息）：

image.png

mkfastq拆分数据这一步后的测序文件结构：

两条reads，read1上是barcode和umi信息，read2上是转录本序列

image.png
ll查看文件：
image.png
从文件名来看pbmc_1k_v3_S1_L001_R1_001.fastq.gz中包含的信息是sample名pbmc_1k，10x测序版本v3，lane通道L001，read1序列R1

less pbmc_1k_v3_S1_L001_R1_001.fastq.gz

image.png
从上图可以看出，原始测序fastq文件结构：每四行为一个单位，第三行是碱基序列，数一下总共有28个碱基，包含16bp barcode和12bp umi；

1.10x 3`端v2版本
read1序列中16bp barcode,10bp umi
read2为98bp

image.png

2.10x 3`端v3版本
read1序列中16bp barcode,12bp umi
read2为91bp

image.png
3.10x V(D)J版本
read1序列中16bp barcode,10bp umi
image.png

二、count细胞定量

cellranger count --id=run_count_1kpbmcs \
   --fastqs=/home/data/10x_test/pbmc_1k_v3_fastqs \
   --sample=pbmc_1k_v3 \
   --transcriptome=/home/data/GRCh38_2020_A/refdata-gex-GRCh38-2020-A

探究一下每个参数：
--id=run_count_1kpbmcs：“随便”起的名字，是之后结果文件的目录名字
--sample=pbmc_1k_v3：样本名
--transcriptome=/home/data/GRCh38_2020_A/refdata-gex-GRCh38-2020-A：参考基因组（转录组）的目录

前两个都很好理解，第三个参数和序列比对等有关

image.png

这里是下载的已经构建好的基因组，如果想自己构建的话：

# 下载基因组
wget ftp://ftp.ensembl.org/pub/release-84/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
gunzip Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
# 下载注释
wget ftp://ftp.ensembl.org/pub/release-84/gtf/homo_sapiens/Homo_sapiens.GRCh38.84.gtf.gz
gunzip Homo_sapiens.GRCh38.84.gtf.gz
# 软件构建注释
# mkgtf <input_gtf> <output_gtf> [--attribute=KEY:VALUE...]
cellranger mkgtf Homo_sapiens.GRCh38.84.gtf Homo_sapiens.GRCh38.84.filtered.gtf \
                --attribute=gene_biotype:protein_coding \
                --attribute=gene_biotype:lincRNA \
                --attribute=gene_biotype:antisense \
                --attribute=gene_biotype:IG_LV_gene \
                --attribute=gene_biotype:IG_V_gene \
                --attribute=gene_biotype:IG_V_pseudogene \
                --attribute=gene_biotype:IG_D_gene \
                --attribute=gene_biotype:IG_J_gene \
                --attribute=gene_biotype:IG_J_pseudogene \
                --attribute=gene_biotype:IG_C_gene \
                --attribute=gene_biotype:IG_C_pseudogene \
                --attribute=gene_biotype:TR_V_gene \
                --attribute=gene_biotype:TR_V_pseudogene \
                --attribute=gene_biotype:TR_D_gene \
                --attribute=gene_biotype:TR_J_gene \
                --attribute=gene_biotype:TR_J_pseudogene \
                --attribute=gene_biotype:TR_C_gene
# 软件利用构建好的注释，去构建需要的基因组
cellranger mkref --genome=GRCh38 \
                --fasta=Homo_sapiens.GRCh38.dna.primary_assembly.fa \
                --genes=Homo_sapiens.GRCh38.84.filtered.gtf \
                --ref-version=1.2.0

学过宏基因组注释的都知道，这和基因组比对非常相似，这里不得不学习一下cellranger count的具体原理，打开软件运行日志：
？看不太懂，换个软件学习。
总之，cellranger count后就能得到一个很熟悉的文件夹，like this:

image.png

之后就是熟悉的seurat或者scanpy分析。

软件二：umitools

这也是一个可以处理10X测序平台的单细胞上游分析软件。
第一步：根据read1的信息构建barcode白名单

umi_tools whitelist --stdin /home/data/10x_test/pbmc_1k_v3_fastqs/pbmc_1k_v3_S1_L001_R1_001.fastq.gz \
        --bc-pattern=CCCCCCCCCCCCCCCCNNNNNNNNNNNN \
        --log2stderr > whitelist.txt

第二步：根据白名单给细胞打上细胞标签

umi_tools extract --bc-pattern=CCCCCCCCCCCCCCCCNNNNNNNNNNNN \
                  --stdin /home/data/10x_test/pbmc_1k_v3_fastqs/pbmc_1k_v3_S1_L001_R1_001.fastq.gz \
                  --stdout pbmc_1k_v3_S1_L001_R1_extracted.fastq.gz \
                  --read2-in /home/data/10x_test/pbmc_1k_v3_fastqs/pbmc_1k_v3_S1_L001_R2_001.fastq.gz \
                  --read2-out=pbmc_1k_v3_S1_L001_R2_extracted.fastq.gz \
                  --filter-cell-barcode \
                  --whitelist=whitelist.txt

第三步：STAR参考基因组比对

/home/yanyt/01.software/seeksoultools.1.2.0/bin/STAR --runThreadN 4 \
        --genomeDir /home/data/GRCH38/GRCh38/star \
        --readFilesIn pbmc_1k_v3_S1_L001_R2_extracted.fastq.gz \
        --readFilesCommand zcat \
        --outFilterMultimapNmax 1 \
        --outSAMtype BAM SortedByCoordinate

第四步：基因定量

featureCounts -a /home/data/GRCh38_2020_A/refdata-gex-GRCh38-2020-A/genes/genes.gtf -o gene_assigned -R BAM Aligned.sortedByCoord.out.bam -T 4
samtools sort Aligned.sortedByCoord.out.bam.featureCounts.bam -o assigned_sorted.bam
samtools index assigned_sorted.bam
umi_tools count --per-gene --gene-tag=XT --assigned-status-tag=XS --per-cell -I assigned_sorted.bam -S counts.tsv.gz

STAR是一个转录组比对的软件，如果是基因组一般使用blast或者diamond;查看了一下cellranger count那一步也是调用了STAR与参考基因组进行比对。

太难了。

根据我的宏基因组数据和单细胞测序数据的不同，补充学习了生信测序中常用的几种文件格式：
一、fasta文件
宏基因组分析时使用的文件：

image.png

后缀可以是：fa/fasta/fna
结构：两行为一个单位，首行以>开头，存储序列的描述信息；第二行是序列本身，包括碱基（DNA或者RNA）或者氨基酸（蛋白质）的字符序列
二、fastq文件
单细胞上游测序时看到的文件：

image.png
后缀可以是：fastq/fq
结构：每四行为一个单位；首行@开头，存储序列的名称或者标识符等信息；第二行是序列本身，包括碱基（DNA或者RNA）序列；第三行+；第四行存储碱基质量信息
三、sam/bam文件
STAR比对后的结果文件
image.png

后缀：sam/bam
结构：文件头以@开头，存储测序等信息：

image.png
比对记录行：
image.png
四、暂时未接触的vcf文件
存储基因组变异数据的文本格式，通常用于描述DNA或RNA测序数据中的单核苷酸变异和结构变异
结构：文件头以##开头，存储文件的属性和来源信息
image.png
主体信息：
image.png
在比较文件格式的过程中，发现转录组fa文件中的碱基也是ATCG，根据微薄的高中生物知识表明这是DNA序列碱基，因此单细胞转录组测序得到的文件存储的是cDNA（？）。
那么存在一个疑问，为啥不能使用宏基因组分析中的blast来比对，而是使用STAR呢？
image.png

而且还有一个小问题，基因组测序没有说表达量而是基因的丰度，转录组则出现了基因的表达量，通过知乎提问，解答如下：

image.png

隐隐约约好像说了点什么，这得回到高中生物那个中心法则，DNA表达了才是RNA，因此转录组从一定程度上刻画基因的表达量（个人理解）。

我又去学习了一下转录组的测序流程：

首先转录组也有质控，也就是fastp、fastqc等；
然后就是比对，因为测序得到的原始文件中没有基因名字，只有一些碱基片段，因此需要一个参考基因组，根据序列相似性，将测序的碱基片段比对到基因区间上去，并根据gtf/gff文件赐予它们名字；
有了基因信息后，又返回去"比对"，计算每个样本中基因的丰度/表达量。
而单细胞转录组不同的地方在于它是单细胞，多了一个对细胞barcode的判别。

此时，我又接触了一篇文献：
Mapping single-cell atlases throughout Metazoa unravels cell type evolution - PubMed (nih.gov)

image.png

可厉害了，可以说是涉及了单细胞转录组上下游以及算法。
在细胞类型映射细胞类型时候，不仅考虑细胞与细胞之间的关系，还要考虑基因与基因之间的关系。
细胞与细胞的关系很简单，可以根据基因表达矩阵计算细胞与细胞的相似性作为距离矩阵；
而基因与基因的关系呢？这里是使用tblastx比对来得到基因与基因的相似性的：
需要准备单细胞的上游参考基因组STAR比对后的fasta文件，将两个单细胞上游测序fasta文件做tblastx比对，将基因与基因之间的相似性作为基因的距离矩阵；
此后，通过某种算法，即考虑细胞与细胞之间基因表达量特征，又通过基因的距离矩阵校正基因的权重，实现数据集间细胞类型的映射关系。
这种算法通常有2个用处：
第一种是比较2个已经注释好细胞类型的数据集之间细胞类型的演化关系：

image.png

第二种就是给未知细胞类型的数据集打上细胞类型标签，方法是将未知细胞类型的数据集的metadata中的cluster赋值为unknow，再进行映射，得到映射.csv文件后通过左连接left_join重新将标签赋值给未知数据集的metadata。

总结：以上基于初学者的理解会存在错误。

参考：
单细胞实战(三) Cell Ranger使用初探 (qq.com)
umi_tools for dealing with UMIs and cell barcodes - 简书 (jianshu.com)
10x Genomics文库结构知多少【3' v2 & V(D)J】 - 简书 (jianshu.com)
生物信息学——常见的四种文件格式(fasta,fastq,sam,vcf)_fastq文件-CSDN博客