单细胞实战(5)：复现拟南芥单细胞文章中的数据（1）

2020-10-21 本文已影响0人周小钊

前言

目前我的课题是植物方面的单细胞测序，所以打算选择植物类的单细胞测序数据进行复现，目前选择了王佳伟老师的《A Single-Cell RNA Sequencing Profiles the Developmental Landscape of Arabidopsis Root》，希望能够得到好的结果

原始数据的下载

首先下载测序数据

prefetch SRR8485805 -O wang/
fastq-dump --split-files SRR8485805
mv SRR8485805_1.fastq data/WT_S1_L001_I1_001.fastq
mv SRR8485805_2.fastq data/WT_S1_L001_R1_001.fastq
mv SRR8485805_3.fastq data/WT_S1_L001_R2_001.fastq

下载基因组与注释文件，需要注意文献中基因组使用的是TAIR10，注释文件是Araport11。

将gff转为gtf文件

gffread Araport11.gff3 -T -o Araport11.gtf

cellranger进行比对

下载cellranger2.2版本

curl -o cellranger-2.2.0.tar.gz "https://cf.10xgenomics.com/releases/cell-exp/cellranger-2.2.0.tar.gz?Expires=1603141363&Policy=eyJTdGF0ZW1lbnQiOlt7IlJlc291cmNlIjoiaHR0cHM6Ly9jZi4xMHhnZW5vbWljcy5jb20vcmVsZWFzZXMvY2VsbC1leHAvY2VsbHJhbmdlci0yLjIuMC50YXIuZ3oiLCJDb25kaXRpb24iOnsiRGF0ZUxlc3NUaGFuIjp7IkFXUzpFcG9jaFRpbWUiOjE2MDMxNDEzNjN9fX1dfQ__&Signature=en6P4Wedmwc2aSEitfKsQp2PITVYKgRPZdzR-fEmjBl4R9yQY5QBQY05--1v8AzRD9WqfoCnddSzFvngrlwxzeCJtFyfHLa2a7ONnUT6NtrzU6RkIj1jwXpaN4NpixnCbEF-Ubj9UZX63W1rEreM0AMNdWiVneGx4bcTajl1KRWaoTNS970DSJ1wrw0g70JFQ0BAltou-qPAeZpD9Xe9EM35EdWRT6eFq~zOaCMRLTxlBjZaMItyDRH~Qecz-B5tLWcAjCKfy4o2hAWTopRRpy93LVV-x1ykxCiHpej5AuAODvUx0V73rZOkRlijcpA5d1rHV~eEdPiM1uoCOJMiSw__&Key-Pair-Id=APKAI7S6A5RYOXBWRPDA"
tar -zxvf cellranger-2.2.0.tar.gz

建立索引并比对

/datadisk02/ScRNAseq_data/cellranger-2.2.0/cellranger  mkgtf --genome=ref --fatsa=TAIR10.fa --genes=Araport11.gtf
                 
/datadisk02/ScRNAseq_data/cellranger-2.2.0/cellranger count --id=WANG --transcriptome=ref --fastqs=data/ --sample=WT --force-cells=8000

比对结果还是可以的，与原文献中差距很小

使用Seurat对数据进行分析

文献中使用到的Seurat为V3版本，要注意cellrangeV2在filtered_gene_bc_matrices生成的文件是genes、barcodes以及matrix，但Seurat识别的是features，我们需要自行对genes文件改名

cd WANG/outs/filtered_gene_bc_matrices/ref
gzip genes.tsv
gzip matrix.mtx
gzip barcodes.tsv
mv genes.tsv.gz features.tsv.gz

创建Seurat对象

library(Seurat)
library(dplyr)
library(ggplot2)
library(magrittr)
library(gtools)
library(stringr)
library(Matrix)
library(tidyverse)
library(patchwork)
setwd("D://data/ScRNAcode/wang/")
##=======================1.创建Seurat对象========================
dir <- 'filtered_gene_bc_matrices/ref/'
counts <- Read10X(dir)
wang = CreateSeuratObject(counts, project = "zxz", min.cells=3, min.features = 200)
dim(wang)
[1] 23228  8000

数据质控与标准化

##=======================2.数据质控与标准化================================
##dir.create('QC')
##提取线粒体基因
wang[["percent.mt"]] <- PercentageFeatureSet(wang, pattern='^ATMG')
violin <- VlnPlot(wang,
                  features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), 
                  pt.size = 0.1, #不需要显示点，可以设置pt.size = 0
                  ncol = 3)
ggsave("QC/vlnplot-before-qc.pdf", plot = violin, width = 15, height = 6) 
ggsave("QC/vlnplot-before-qc.png", plot = violin, width = 15, height = 6) 
plot1 <- FeatureScatter(wang, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(wang, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
pearplot <- CombinePlots(plots = list(plot1, plot2), nrow=1, legend="none") 
ggsave("QC/pearplot-before-qc.pdf", plot = pearplot, width = 12, height = 5) 
ggsave("QC/pearplot-before-qc.png", plot = pearplot, width = 12, height = 5)
##设置质控标准
wang<-subset(wang,subset=nFeature_RNA>500 & nFeature_RNA<5000 &percent.mt<0.5)
dim(wang)
[1] 23228  7626
## 绘制质量控制后的图
violin <-VlnPlot(wang,
                 features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), 
                 pt.size = 0.1, 
                 ncol = 3)
ggsave("QC/vlnplot-after-qc.pdf", plot = violin, width = 15, height = 6) 
ggsave("QC/vlnplot-after-qc.png", plot = violin, width = 15, height = 6)
## 基因表达量标准化
## 它的作用是让测序数据量不同的细胞的基因表达量具有可比性。计算公式如下：
## 标准化后基因表达量 = log1p（10000*基因counts/细胞总counts）
wang <- NormalizeData(wang, normalization.method = "LogNormalize", scale.factor = 10000)

质控后细胞数目为7626，基因数为23228，原文献中两者的数据分别是7695与23161

数据降维与聚类

##=======================3.数据降维与聚类==================================
## 寻找高变基因
## dir.create("cluster")
wang <- FindVariableFeatures(wang,mean.cutoff=c(0.0125,3),dispersion.cutoff =c(1.5,Inf) )
top10 <- head(VariableFeatures(wang), 10)
plot1 <- VariableFeaturePlot(wang) 
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE, size=2.5) 
plot <- CombinePlots(plots = list(plot1, plot2),legend="bottom")
## 横坐标是某基因在所有细胞中的平均表达值，纵坐标是此基因的方差。
## 红色的点是被选中的高变基因，黑色的点是未被选中的基因，变异程度最高的10个基因在如图中标注了基因名称。
ggsave("cluster/VariableFeatures.pdf", plot = plot, width = 8, height = 6) 
ggsave("cluster/VariableFeatures.png", plot = plot, width = 8, height = 6)
## 数据缩放
scale.genes <-  rownames(wang)
wang <- ScaleData(wang, features = scale.genes)
## PCA降维并提取主成分
wang <- RunPCA(wang, features = VariableFeatures(wang),npcs = 100) 
plot1 <- DimPlot(wang, reduction = "pca") 
plot2 <- ElbowPlot(wang, ndims=40, reduction="pca") 
plotc <- plot1+plot2
ggsave("cluster/pca.pdf", plot = plotc, width = 8, height = 4) 
ggsave("cluster/pca.png", plot = plotc, width = 8, height = 4)
## 细胞聚类
## 此步利用 细胞-PC值 矩阵计算细胞之间的距离，
## 然后利用距离矩阵来聚类。其中有两个参数需要人工选择，
## 第一个是FindNeighbors()函数中的dims参数，需要指定哪些pc轴用于分析，选择依据是之前介绍的cluster/pca.png文件中的右图。
## 第二个是FindClusters()函数中的resolution参数，需要指定0.1-1.0之间的一个数值，用于决定clusters的相对数量，数值越大cluters越多。
wang <- FindNeighbors(object = wang, dims = 1:100)
wang <- FindClusters(object = wang, resolution = 1.0)
table(wang@meta.data$seurat_clusters)
## 非线性降维
## tsne
wang <- RunTSNE(wang, dims =1:40)
embed_tsne <- Embeddings(wang, 'tsne')
write.csv(embed_tsne,'cluster/embed_tsne_new.csv')
plot1 = DimPlot(wang, reduction = "tsne" ,label = "T", pt.size = 1,label.size = 4)
ggsave("cluster/tSNE_cluster.pdf", plot = plot1, width = 8, height = 7)
ggsave("cluster/tSNE_cluster.png", plot = plot1, width = 8, height = 7)

## UMAP'
wang <- RunUMAP(wang,n.neighbors = 30,metric = 'correlation',min.dist = 0.3,dims = 1:40)
embed_umap <- Embeddings(wang, 'umap')
write.csv(embed_umap,'cluster/embed_umap_new.csv')
plot2 = DimPlot(wang, reduction = "umap",label = "T", pt.size = 1,label.size = 4) 
ggsave("cluster/UMAP_cluster_new.pdf", plot = plot2, width = 8, height = 7)
ggsave("cluster/UMAP_cluster_new.png", plot = plot2, width = 8, height = 7)

结果是有区别的，我的聚类比原文献中要多一个，而且数字不对应，所以我要用文献中列出的某些基因的小提琴图确定我的聚类

根据文献对应自己数据聚类

原文献中有所有聚类的特异基因，所以我根据列出的基因去匹配我的聚类结果

##==============================5.修改聚类标号=====================
##修改聚类号重新做图
new.cluster.ids<-c("2",'1','4','5','13','3','12','21','8','6','11',
                   '9','7','10','6','15','22','14','17','19','16',
                   '20','18','23','24')
names(new.cluster.ids) <- levels(wang)
wang <- RenameIdents(wang, new.cluster.ids)
Idents(wang)<-factor(Idents(wang),levels=mixedsort(levels(Idents(wang))))
wang <- RunTSNE(wang, dims =1:40)
embed_tsne <- Embeddings(wang, 'tsne')
write.csv(embed_tsne,'cluster/embed_tsne-new.csv')
plot1 = DimPlot(wang, reduction = "tsne" ,label = "T", pt.size = 1,label.size = 4)
ggsave("cluster/tSNE_cluster-new.pdf", plot = plot1, width = 8, height = 7)
ggsave("cluster/tSNE_cluster-new.png", plot = plot1, width = 8, height = 7)

## UMAP
wang <- RunUMAP(wang,n.neighbors = 30,metric = 'correlation',min.dist = 0.3,dims = 1:40)
embed_umap <- Embeddings(wang, 'umap')
write.csv(embed_umap,'cluster/embed_umap-new.csv')
plot2 = DimPlot(wang, reduction = "umap",label = "T", pt.size = 1,label.size = 4) 
ggsave("cluster/UMAP_cluster.pdf", plot = plot2, width = 8, height = 7)
ggsave("cluster/UMAP_cluster.png", plot = plot2, width = 8, height = 7)