10X小鼠单细胞转录组和空间转录组预处理
1.数据下载
GEO 数据下载：进入EBI数据库网站：The European Bioinformatics Institute < EMBL-EBI，找到自己感兴趣的数据的下载连接
image.png
2.linux下载命令 下载超快！！(但是由于单细胞数据比较大，最好用迅雷或者更稳健的下载器，不然下载不完全就有可能停)
1.nohup ascp -QT -l 20m -P33001 -i /home/WHL/anaconda2/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:vol1/srr/SRR144/097/SRR14458797 ./ &
2.nohup wget -c -t 0 -O path/SRR1482463.sra https://sra-downloadb.be-md.ncbi.nlm.nih.gov/sos1/sra-pub-run-5/SRR14458797/SRR14458797 &
#-c -t 配合使用可以防止下载数据的过程中链接中断的问题，-O则可以指定下载路径和文件名。
3.原始下机数据预处理
3.1. .sra转fastq.gz文件
nohup fastq-dump --gzip --split-files /home/WHL/singalcell_test/result/SRR14458797.sra &
#运行后得到两个fastq文件：SRR14458797_1.fastq.gz、SRR14458797_2.fastq.gz、SRR14458797_3.fastq.gz
3.2.Cell Ranger： 10X公司对单细胞 RNA 测序输出结果进行比对、定量、聚类及基因表达分析的分析流程
Cellranger原理介绍（上） - 简书 (jianshu.com)
单细胞实战(二) cell ranger使用前注意事项 (qq.com)
##小鼠参考基因组下载：
curl -O https://cf.10xgenomics.com/supp/cell-exp/refdata-gex-mm10-2020-A.tar.gz
cellranger count --id=sample345 \
--transcriptome=/home/WHL/singalcell_test/mm10refseq/refdata-gex-mm10-2020-A \
--fastqs=/home/wanghongli/singalcell_test/result \
--sample=mysample \
--expect-cells=1000 \
--nosecondary
# id指定输出文件存放目录名
# transcriptome指定与CellRanger兼容的参考基因组
# fastqs指定mkfastq或者自定义的测序文件
# sample要和fastq文件的前缀中的sample保持一致，作为软件识别的标志
# expect-cells指定复现的细胞数量，这个要和实验设计结合起来
# nosecondary 只获得表达矩阵，不进行后续的降维、聚类和可视化分析(因为后期会自行用R包去做)
结果解读：Cellranger原理介绍（下） - 简书 (jianshu.com)
image.png
image.png
3.3.多个文库的整合 aggr(重复样本合并)
当处理多个生物学样本或者一个样本存在多个重复/文库时，最好的操作就是先分别对每个文库进行单独的count定量，然后将定量结果利用aggr组合起来
3.3.1.第一步 得到count结果
例如现在分别进行3个定量流程
$ cellranger count --id=LV123 ...
... wait for pipeline to finish ...
$ cellranger count --id=LB456 ...
... wait for pipeline to finish ...
$ cellranger count --id=LP789 ...
... wait for pipeline to finish ...
3.3.2.第二步 构建Aggregation CSV
# AGG123_libraries.csv
library_id,molecule_h5
LV123,/opt/runs/LV123/outs/molecule_info.h5
LB456,/opt/runs/LB456/outs/molecule_info.h5
LP789,/opt/runs/LP789/outs/molecule_info.h5
# 其中
molecule_h5：文件molecule_info.h5 file的路径
3.3.3.第三步 运行aggr
cellranger aggr --id=AGG123 \
--csv=AGG123_libraries.csv \
--normalize=mapped
# 结果输出到AGG123这个目录中
新版本的cellranger可以实现细胞过滤无需用 DropletUtils进行过滤
亲测：
1.cellranger raw_feature_bc_matrix数据绘制的图
未过滤的（图1）：
image.png
过滤完之后的（图1.2）：
image.png
2.cellranger filtered_feature_bc_matrix数据绘制的图
未执行过滤的（图2.1）：
image.png
如果过滤的话会报如下错误：
image.png
对比1.2和2.1会发现几乎一样
3.4 DropletUtils简介
OSCA单细胞数据分析笔记-14、Empty/Doublet droplet - 简书 (jianshu.com)
3.4.1 DropletUtils介绍
R包DropletUtils针对10X Genomics平台，根据观察到的每个液滴的表达谱与周围溶液的表达谱来区分空液滴（empty droplets，只含溶液中RNA）和含细胞的液滴。
DropletUtils主要功能如下：
1.读入10X Genomics平台的UMI count matrix
2.读入CellRanger生成的molecule information file (molecule_info.h5)
3.降采样：downsampling the UMI count matrix or the raw reads
4.识别空液滴和doublets
5.去除Illumina 4000测序仪的barcode swapping效应
4 用Seurat进行数据标准化、降维、聚类、
安装Seurat包
install.packages("Seurat", repos="https://mran.microsoft.com/snapshot/2019-02-01/")
install.packages(c('dplyr','patchwork'))
install.packages('magrittr')
install.packages('hdf5r')
library(Seurat)
library(dplyr)
library(magrittr)
library(patchwork)
library(hdf5r)
读取目录下包括（barcodes.tsv.gz、features.tsv.gz、matrix.mtx.gz）读取数据主要使用Read10X这个函数
pbmc98f.data <- Read10X(data.dir = "D:\\scRNA-seq\\filtered_feature_bc_matrix")用原始数据（non-normalized data）
创建Seurat对象
使用CreateSeuratObject函数创建Seurat对象，这里要求细胞中基因数目不能小于200，且基因至少在三个细胞中有表达，否则过滤掉
pbmc98_f <- CreateSeuratObject(counts = pbmc98f.data, project = "pbmc98f", min.cells = 3, min.features = 200)
pbmc98_f
# An object of class Seurat
# 20579 features across 4901 samples within 1 assay
# Active assay: RNA (20579 features, 0 variable features)
# 总共得到4901个细胞和20579个基因
# nCount_RNA 每个细胞的UMI数目
# nFeature_RNA 每个细胞所检测到的基因数目
Seurat对象就是一个S4类，里面装着单细胞数据集，如count矩阵、细胞矩阵、聚类信息都存储在这个容器中。
我们可以通过查看这个对象来得到我们想要的信息。如我们想看一下前10个细胞的前10个基因的表达矩阵，可以这样使用：
pbmc98_f@assays$RNA@counts[1:10,1:10]
image.png
提取线粒体基因集
pbmc98_f[["percent.mt"]] <- PercentageFeatureSet(pbmc98_f, pattern = "^MT-")
qc指标可视化
VlnPlot(pbmc98_f, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
image.png
计算每个特征之间的相关性
plot1 <- FeatureScatter(pbmc98_f, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot2 <- FeatureScatter(pbmc98_f, feature1 = "nCount_RNA", feature2 = "percent.mt")
image.png
从上图的可视化结果进行细胞选择。我们过滤具有超过2500或少于200个基因计数的细胞，同时过滤掉线粒体比例超过5%的细胞。
pbmc98_f <- subset(pbmc98_f, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
从数据集中移除不需要的细胞后，下一步是标准化数据。默认情况下，我们使用“lognormalize”全局缩放的归一化方法，通过总表达值对每个细胞的基因表达值归一化，并将其乘以缩放因子(默认为10,000)，最后对结果进行对数变换，标准化数据存储在pbmc98_f @assays$RNA@data中。
pbmc98_f <- NormalizeData(pbmc98_f, normalization.method = "LogNormalize", scale.factor = 10000)
pbmc98_f@assays$RNA@counts[1:5,1:10]
计算数据集中表现出高细胞间差异的基因子集(即，它们在一些细胞中高表达，而在另一些细胞中低表达)。使用均值与方差之间的关系，来挑选高变基因，在下游分析中关注这些基因有助于突出单细胞数据集中的生物信号。默认情况下，每个数据集选择2000个高变异基因用于下游分析。
pbmc98_f <- FindVariableFeatures(pbmc98_f, selection.method = "vst", nfeatures = 2000)
示例：提取前10的高变异基因，并将前10的高变异基因在图中标注出来
top10 <- head(VariableFeatures(pbmc98_f), 10)
plot1 <- VariableFeaturePlot(pbmc98_f)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
image.png
PCA
使用线性变换(“scaling”)处理数据，过程：改变每个基因的表达，使细胞间的平均表达为0；缩放每个基因的表达，使细胞间的差异为1，这样高表达基因就不会影响后续分析；这步是pca等降维技术之前的标准预处理步骤，缩放的数据储存在pbmc98_f@assays$RNA@scale.data中。
#默认只对2千个高变基因进行数据归一化
pbmc98_f <- ScaleData(pbmc98_f)
# 也可以使用全部的基因进行归一化，但是耗时会比较久
all.genes <- rownames(pbmc98_f)
pbmc98_f <- ScaleData(pbmc98_f, features = all.genes)
接下来，我们对缩放后的数据执行pca降维，默认降至50维。PCA(Principal Component Analysis)，即主成分分析方法，是一种使用最广泛的数据降维算法。PCA的主要思想是将n维特征映射到k维上，在尽可能保留原始数据信息的同时降低数据维度来加速数据分析。过程就是从原始高维的空间中按顺序地找一组相互正交的坐标轴系统，新的坐标轴的选择与数据本身是密切相关的。其中，第一个新坐标轴选择是原始数据中方差最大的方向，第二个新坐标轴选取是与第一个坐标轴正交的平面中使得方差最大的，第三个轴是与第1,2个轴正交的平面中方差最大的。依次类推，可以得到n个这样的坐标轴。大部分方差都包含在前面k个坐标轴中，后面的坐标轴所含的方差几乎为0。于是，我们可以忽略余下的坐标轴，只保留前面k个含有绝大部分方差的坐标轴，实现对数据降维。
pbmc98_f <- RunPCA(pbmc98_f, features = VariableFeatures(object = pbmc98_f))
Seurat的VariableFeatures中包含scran检测到的可变基因
#选择前5个维度进行查看
print(pbmc98_f[["pca"]], dims = 1:5, nfeatures = 5)
使用DimPlot可视化函数查看样品中的细胞在前pca中的分布情况：
DimPlot(pbmc98_f, reduction = "pca")
#可视化 前4个pca相关的top基因
VizDimLoadings(pbmc98_f, dims = 1:4, reduction = "pca")
image.png
image.png
使用DimHeatmap可视化函数查看样品中前500个细胞在前6个PCA中的热图：
DimHeatmap(pbmc98_f, dims = 1:6, cells = 500, balanced = TRUE)
image.png
每个pc本质上代表一个“元特征”，它将相关特征集中的信息组合在一起。因此，越在顶部的主成分越可能代表数据集。然而，我们应该选择多少个主成分才认为我们选择的数据包含了绝大部分的原始数据信息呢？我们可以使用JackStraw函数来选择PC数目：
#将数据的1%打乱重新运行pca，构建特征得分的“零分布”，这个过程重复100次
pbmc <- JackStraw(pbmc98_f, num.replicate = 100)
pbmc1 <- ScoreJackStraw(pbmc, dims = 1:10)
plot1<-JackStrawPlot(pbmc1, dims = 1:10)
plot2<-ElbowPlot(pbmc1)
CombinePlots(plots = list(plot1, plot2))
image.png
Seurat3使用基于图聚类算法（graph-based clustering approach）进行细胞聚类分析。这部分聚类按我的理解是在pca空间中分布着各个细胞，首先在这个空间中的计算离每个细胞最近的k个细胞的欧氏距离来构造一个个knn图，并基于任意两个细胞在其局部邻域中的共享重叠来细化它们之间的边缘权重，尝试将该图划分为高度互连的社区，最后应用模块化优化技术（默认Louvain 算法）将细胞迭代分组在一起。该步骤使用FindNeighbors和FindClusters两个函数完成：
#计算最邻近距离
pbmc2 <- FindNeighbors(pbmc1, dims = 1:15)
#聚类，包含设置下游聚类的“间隔尺度”的分辨率参数resolution ，增加值会导致更多的聚类。
pbmc3 <- FindClusters(pbmc2, resolution = 0.5)
# 查看前5个细胞的聚类id
head(Idents(pbmc3), 5)
seurat提供了几种非线性降维技术，例如tsne和umap，来可视化和探索这些数据集。我们将上一步PC维度中的数据进一步降至二维，实现数据的可视化。之前pca是线性降维，能很好处理高维度数据，tsne和umap是非线性降维方式，对低维度数据处理速度快，但难以处理高维度数据。相比tsne，umap能更好地反映原始数据结构，有着更好的连续性。
pbmc3 <- RunUMAP(pbmc3, dims = 1:15)
pbmc3 <- RunTSNE(pbmc3, dims = 1:15)
DimPlot(pbmc3, reduction = "umap")
DimPlot(pbmc3, reduction = "tsne")
image.png
image.png
seurat可以通过差异表达找到每个聚类的markgene，差异分析可以有多种形式，如找到所有聚类的markene（如cluster1中所有的markgene是指cluster1相对于其余所有cluster是差异的）、两个cluster之间的差异分析、某个cluster中两个样品之间差异分析等。
#找到cluster1中的markgene
#其中ident.1参数设置待分析的细胞类别，min.pct表示该基因表达数目占该类细胞总数的比例
cluster1.markers <- FindMarkers(pbmc3, ident.1 = 1, min.pct = 0.9)
#找到cluster5和cluster0、3之间的markgene
cluster5.markers <- FindMarkers(pbmc3, ident.1 = 5, ident.2 = c(0, 3), min.pct = 0.25)
#找到每个cluster相比于其余cluster的markgene，只报道阳性的markgene
pbmc.markers <- FindAllMarkers(pbmc3, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
如果你对每个cluster已经鉴定出细胞类型，就可以在可视化中将细胞类型标注在图中
new.cluster.ids <- c("Naive CD4 T", "Memory CD4 T", "CD14+ Mono", "B", "CD8 T","lll")
names(new.cluster.ids) <- levels(pbmc3)
pbmc3 <- RenameIdents(pbmc3, new.cluster.ids)
DimPlot(pbmc3, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
image.png
查看特定基因在聚类图中的表达量分布情况
FeaturePlot(pbmc3, features = c("Dusp5", "Rnf31"))
image.png
绘制每个cluster按logfc排名前10的markgene的热图
top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 50, wt = avg_log2FC)
DoHeatmap(pbmc3, features = top10$gene) + NoLegend()
image.png
查看特定基因在cluster中的表达山脊图
features <- c("Tuba1b", "Rps2")
RidgePlot(pbmc3, features = features, ncol = 2)
image.png
查看特定基因在cluster中的分布点图，点的大小代表表达该基因的细胞比例，颜色代表平均表达水平
features <- c("Tuba1b", "Rps2")
DotPlot(pbmc3, features = features) + RotatedAxis()
image.png
利用scMCA进行细胞注释
安装scMCA R包（我github在线安装scMCA包失败 下载到本安装）
install.packages('devtools')
install.packages("rJava")
install.packages("usethis")
install.packages("shinythemes")
library(devtools)
library(rJava)
library(usethis)
library(shinythemes)
devtools::install_local("C:\\Users\\Administrator\\Desktop\\scMCA-master.zip")
library(scMCA)
提取表达矩阵GetAssayData()
从Assay中提取
d <- as.matrix(GetAssayData(object = pbmc3[["RNA"]], slot = "data"))
#前10个基因和前20个细胞
##②从Seurat对象中提取
f <- as.matrix(GetAssayData(pbmc3[["RNA"]], slot = "data")[1:10,1:20])
#提取lll细胞的表达矩阵，用于其他分析,lll是之前注释的cluster
nk.raw.data <- as.matrix(GetAssayData(pbmc3[["RNA"]], slot = "data")[, WhichCells(pbmc3, ident = "lll")])
mca_result <- scMCA(scdata = d, numbers_plot = 1)
scMCA_vis(mca_result)
image.png
10X空间转录组数据处理
示例数据下载(需要fastq文件和tif或jpg文件)
nohup wget https://s3-us-west-2.amazonaws.com/10x.files/samples/spatial-exp/1.1.0/V1_Adult_Mouse_Brain/V1_Adult_Mouse_Brain_fastqs.tar &
wget https://cf.10xgenomics.com/samples/spatial-exp/1.1.0/V1_Adult_Mouse_Brain/V1_Adult_Mouse_Brain_image.tif
空间转录组原始数据处理
spaceranger下载：
curl -o spaceranger-1.0.0.tar.gz "http://cf.10xgenomics.com/releases/spatial-exp/spaceranger-1.0.0.tar.gz?Expires=1575402715&Policy=eyJTdGF0ZW1lbnQiOlt7IlJlc291cmNlIjoiaHR0cDovL2NmLjEweGdlbm9taWNzLmNvbS9yZWxlYXNlcy9zcGF0aWFsLWV4cC9zcGFjZXJhbmdlci0xLjAuMC50YXIuZ3oiLCJDb25kaXRpb24iOnsiRGF0ZUxlc3NUaGFuIjp7IkFXUzpFcG9jaFRpbWUiOjE1NzU0MDI3MTV9fX1dfQ__&Signature=fACB1rzbHv1rwUicNqL8SheRe6FkFOKxow5cTXcZPfOPBOTBEplElFMnOi4Xv4A2X3kydX45B-JnIaRj7I6a2doGEMTyqv84BnM5LxHAVBtWrXJyQqXbKKtgl9Dxe4BDnM9rPKhs6o2UbmWWAHX8Xu4J3~vgP3yXbhovuyl6OqCxu5p82oxTeOfN0bONqZdZ33svlAXJhatUTdpse2YCSRJZzov69NSHF6gE5DXl6iu5RWU7AgnjFgCuEFkQMwyn-FoYi2~i0s2fOFK0RCVI07~YKNDsjz3eXgOoHjWGPtWw5DAbPpTB2~32xkGzYeIYeZjH6m5JEgNGuvfWEyj~Aw__&Key-Pair-Id=APKAI7S6A5RYOXBWRPDA"
##解压安装设置环境变量
tar -xzvf spaceranger-1.0.0.tar.gz
vi ~/.bashrc # 填入 export PATH=/opt/spaceranger-1.0.0:$PATH
source ~/.bashrc
spaceranger count --id=V1_Adult_Mouse_Brain \
--transcriptome=/home/whl/singalcell_test/mm10refseq/refdata-gex-mm10-2020-A \
--fastqs=/home/whl/spatia_Trans/test2/V1_Adult_Mouse_Brain_fastqs \
--sample=V1_Adult_Mouse_Brain \
--image=/home/whl/spatia_Trans/test2/V1_Adult_Mouse_Brain_image.ti\
--slide=V19L01-041 \
--area=C1 \
--localcores=32 \
--localmem=128
##
--fastqs fastq_path文件夹的路径
--transcriptome Space Ranger 兼容转录组参考的路径
--sample 提供给 的示例工作表中指定的示例名称。
--image .jpg或.tiff格式的明场组织H&E图像。
--slide Visium 幻灯片序列号。必需，除非 --unknown-slide 通过。
--area Visium 捕获区域标识符。必需，除非 --unknown-slide 通过。Visium 的选项包括 A1、B1、C1、D1。
--slidefile 幻灯片版式文件，指示捕获点和基准点位置。
--localcores 限制为使用指定数量的核心来执行管道阶段。默认情况下，将使用系统上所有可用的核心。
--localmem 限制使用指定的内存量（以 GB 为单位）来执行管道阶段。默认情况下，将使用系统上 90% 的可用内存。
Space Ranger软件（10x Genomics）对读数进行解复用并映射到参考基因组，计数矩阵被加载到Seurat和STutility，用于所有后续的数据过滤、归一化、过滤、降维和可视化。使用Seurat的SCTransform功能中实现的方差稳定转换方法。我们使用STutility（nFactors = 20）对归一化的表达矩阵进行了非负矩阵分解。
用Stereoscope空间反卷积
Stereoscope是一个概率模型，用于使用注释的scRNA-seq数据作为输入估计细胞类型比例。
