【单细胞】烟草花冠单细胞数据实战

上面一个帖子分享了这篇文献，因为是做的烟草研究，所以看下能否重复作者的结果。

=====下载原始数据====

GEO num：GSE193464。

一共三组数据。

SRR17555533，SRR17555534，SRR17555535。分别对应ZT8，ZT12，ZT16。

然后转化成fastq格式：

fasterq-dump -O ./ --split-files -e 40 SRR17555533.sra --include-technical

fasterq-dump -O ./ --split-files -e 40 SRR17555534.sra --include-technical

fasterq-dump -O ./ --split-files -e 40 SRR17555535.sra --include-technical

然后修改名字为cellranger要求的输入格式：

Sample Name_S1_L00[Lane Number]_[Read Type]_001.fastq.gz

其中Read Type

I1：Sample index read (optional)

R1：Read 1

R2：Read 2

===准备参考基因组相关文件===

因为文中用的是NIATTr2版本的Nicotiana attenuata的基因组，所以我们先下载这个基因组：

https://plants.ensembl.org/Nicotiana_attenuata/Info/Index

主要下载基因组序列文件和注释的gff文件。

对于基因建立相应的索引文件：

利用cellranger的mkref：

其中genome代表输出的文件夹名字，fasta基因组序列文件， genes是注释文件gtf格式

/gpfs03/home/jingjing/software/cellranger-4.0.0/bin/./cellranger mkref --genome=refdata --fasta=Nicotiana_attenuata.fa --genes=Nicotiana_attenuata.gtf --nthreads=40

===得到表达矩阵===

cellranger count --id ZT8 --fastqs /gpfs03/bioinfo/Proj_SC_2022/Tobacco/001_tob_35075650/SRR17555533/ --sample SRR17555533 --transcriptome /gpfs03/home/jingjing/project/database/N.attenuata/cellranger/refdata/

这个是paper中的比对、表达结果。

这是我自己跑出来的结果，和paper给出的差不多，略有差异，不知道是不是基因组版本用的略有差异还是cellranger版本的差异。主要似乎ZT8差异大点。

但是从结果来看，个人觉得数据的quality不是很高。例如鉴定的细胞数目偏少，目前我们自己在做的情况来看，植物细胞样品现在基本能达到万的量级，1000基本是底限。但是作者的样品鉴定出来的数目还在1000左右，另外从比对率来看，可用的比对到genome的reads才30%左右，也是踩着基本线，而能比对到trancriptome的才20%左右，相当于这个数据可用的reads才20%。对于这样不理想的结果，作者也没提，更没进行进一步的处理和说明。

====聚类和细胞类型鉴定====

library(Seurat)

pbmc8h.data <- Read10X(data.dir = "ZT8/ZT8_filtered_feature_bc_matrix/")

pbmc8h <- CreateSeuratObject(counts = pbmc8h.data, project = "pbmc8h", min.cells = 3, min.features = 200)

pbmc8h

> pbmc8h

An object of class Seurat 

16515 features across 1022 samples within 1 assay 

Active assay: RNA (16515 features, 0 variable features)

pbmc12h.data <- Read10X(data.dir = "ZT12/ZT12_filtered_feature_bc_matrix/")

pbmc12h <- CreateSeuratObject(counts = pbmc12h.data, project = "pbmc12h", min.cells = 3, min.features = 200)

pbmc12h

> pbmc12h

An object of class Seurat 

16736 features across 1035 samples within 1 assay 

Active assay: RNA (16736 features, 0 variable features)

pbmc16h.data <- Read10X(data.dir = "ZT16/ZT16_filtered_feature_bc_matrix/")

pbmc16h <- CreateSeuratObject(counts = pbmc16h.data, project = "pbmc16h", min.cells = 3, min.features = 200)

pbmc16h

> pbmc16h

An object of class Seurat 

16172 features across 1718 samples within 1 assay 

Active assay: RNA (16172 features, 0 variable features)

pbmc.combine <- merge(pbmc8h, y = c(pbmc12h, pbmc16h), add.cell.ids = c("8h", "12h", "16h"), project = "pbmc36h")

pbmc.combine

> pbmc.combine

An object of class Seurat 

18094 features across 3775 samples within 1 assay 

Active assay: RNA (18094 features, 0 variable features)

head(colnames(pbmc.combine))

table(pbmc.combine$orig.ident)

可以看出经过简单filter之后，ZT8h，ZT12h，ZT16h分别获得了1022，1035和1718个cell。

下面进入单细胞数据的常规分析中。

1. 首先进行质控和过滤，但是看paper中好像没有这一步，但是我们还是准备看一下情况。

因为这个基因组没有叶绿体或者线粒体的信息，我们主要过滤一些潜在的离群值。

pbmc <- pbmc.combine

VlnPlot(pbmc, features = c("nFeature_RNA","nCount_RNA"), ncol =2)

前面的帖子也解释过这2个量的含义。#nFeature_RNA：代表的是在该细胞中共检测到的表达量大于0的基因个数，nCount_RNA：代表的是该细胞中所有基因的表达量之和。

plot <- FeatureScatter(pbmc, feature1 ="nCount_RNA", feature2 ="nFeature_RNA")

plot

从nCount和gene之间的关系图可以看到非常明显的一个相关性，当gene个数为5000时对应的umi在50000左右。以nFeature_RNA为例，可以看到数值在50000以上的细胞是非常少的，可以看做是离群值，所以在筛选时，如果一个细胞中检测到的基因个数大于5000，就可以进行过滤。

pbmc <- subset(pbmc, subset = nFeature_RNA >500& nFeature_RNA <5000)

> pbmc

An object of class Seurat 

18094 features across 3197 samples within 1 assay 

Active assay: RNA (18094 features, 0 variable features)

当然也可以不过滤，paper中就没有这一步。因为这个文章中鉴定出的cell数目也不多，我们还是原始的数据进行分析。

2. 进行归一化

预处理之后，首先进行归一化：

pbmc <- NormalizeData(pbmc, normalization.method ="LogNormalize", scale.factor =10000)

默认采用LogNormalize归一化算法，首先将原始的表达量转换成相对丰度，然后乘以scale.factor的值，在取对数。

归一化之后，Seurat提取那些在细胞间变异系数较大的基因用于下游分析:

pbmc <- FindVariableFeatures(pbmc, selection.method ="vst", nfeatures =5000)

将返回5000个高变异的基因。//作者文中也没提自己用的参数，我们暂时选了默认的2000

top10 <- head(VariableFeatures(pbmc), 10)      //返回最高变异的10个基因

plot1 <- VariableFeaturePlot(pbmc)

plot1

这是2000的结果：

从图中可以看出挑选的2000个高变异的基因明显分布和误差比低变异度的高。

plot2 <- LabelPoints(plot = plot1, points = top10)

plot2     //加上前10个高变异基因的label更加清晰

all.genes <- rownames(pbmc)

pbmc <- ScaleData(pbmc, features = all.genes)  //这两步scale，保持同一个基因在不同cell之间的方差为1，均值为0，此步骤在下游分析中具有相同的权重，因此高表达的基因不会占主导地位

pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))

VizDimLoadings(pbmc, dims =1:2, reduction ="pca")

DimPlot(pbmc, reduction = "pca")

DimHeatmap(pbmc, dims = 1:15, cells = 500, balanced = TRUE)  //主要用来查看数据集中的异质性的主要来源，并且可以确定哪些PC维度可以用于下一步的下游分析。

下面我们需要确定哪些主成分所代表的基因可以进入下游分析，这里可以使用JackStraw做重抽样分析。用JackStrawPlot可视化看看哪些主成分可以进行下游分析。（注：paper中采用的参数是top50 PCAs，我不是很明白为啥选了50）

pbmc <- JackStraw(pbmc, num.replicate = 100)

pbmc <- ScoreJackStraw(pbmc, dims = 1:20)

JackStrawPlot(pbmc, dims = 1:20)

==非线性维度约化（UMAP/TSNE）========

pbmc <- FindNeighbors(pbmc, dims = 1:50)

pbmc <- FindClusters(pbmc, resolution = 0.4)

head(Idents(pbmc), 5)

pbmc <- RunUMAP(pbmc, dims = 1:50)

DimPlot(pbmc, reduction = "umap", label = TRUE)

DimPlot(pbmc, reduction = "umap", group.by="orig.ident")   

三个时间点分别显示：

pbmc <- RunTSNE(pbmc, dims = 1:50)

DimPlot(pbmc, reduction = "tsne", label = TRUE)

鉴定各个cluster的marker基因

pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)

pbmc.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_log2FC)

top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC)

DoHeatmap(pbmc, features = top10$gene) + NoLegend()

然后利用paper中给出的marker基因进行细胞类型的注释。

比如：cluster 5可以被注释为pollen

其中两个marker基因的表达箱图：

在聚类图中的分布图：

比如cluster 7：可以被注释为phloem/xylem

cluster 4为chlorenchyma cell，因为很多光合作用相关的基因特异表达。

其它cluster的类型鉴定相对来说复杂的多，比如cluster1/2/3/6:

从富集的marker我是看不出来有明显的规律，如果非要说规律的话就是脂类或者很多次生代谢合成相关的基因明显比例很高。作者在paper中提到的是：Epidermal cells are known to express wax/cutin biosynthetic genes for synthesizing cuticular wax and cutin. Several genes involved in cuticular wax biosynthesis were specifically enriched in clusters 0, 3, 4, 5, and 9: ketoacyl-CoA synthase, ECERIFERUM4, wax ester synthase/diacylglycerol acyltransferase 1, and several lipid-transfer proteins, which may participate in wax synthesis or transport, were highly expressed. We also found that clusters 0, 3, 4, 5, and 9 contained higher levels of cutin biosynthetic genes (glycerol-3-phosphate acyltransferase 6 and cutin synthase-like) compared to other clusters; their protein sequences are similar to the well-studied orthologs of Arabidopsis thaliana. These results suggest that clusters 0, 3, 4, 5, and 9 are epidermal cells.

这个可能与我们分析的聚类cluster1/2/3/6估计是类似的，但是总觉得原因不能那么solid。

最后就是最难处理的cluster 0，我个人是完全看不出来什么规律的。

作者在paper中的描述是：Clusters 1 and 2 had low levels of cuticle biosynthetic genes and relatively high levels
of photosynthetic genes. These results suggest that in the corolla,clusters 1 and 2 are parenchyma cells. 这个在我们的分析中可能指的是cluster 0。但是对于原因，我个人觉得好像不太能说的通。其实在分析单细胞数据的时候，个人觉得这一步不是很简单的一步，即便我们自己做了marker gene的库来辅助细胞类型的鉴定，依然碰到很多问题，不知道该如何处理。

最后修改我们注释的结果，可以获得最终的细胞类型鉴定结果：

new.cluster.ids <- c("Parenchyma","Epidermis","Epidermis","Epidermis","Chlorenchyma","Pollen","Epidermis","Phloem/xylem")

names(new.cluster.ids) <- levels(pbmc)

pbmc <- RenameIdents(pbmc, new.cluster.ids)

DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 1)+ NoLegend()

下面是paper中的分析：

本文使用 文章同步助手 同步