前言
目前我的課題是植物方面的單細胞測序，所以打算選擇植物類的單細胞測序數據進行復現膜宋，目前選擇了王佳偉老師的《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 mkref --genome=ref --fasta=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
創(chuàng)建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.創(chuàng)建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)
結果是有區(qū)別的埋凯，我的聚類比原文獻中要多一個点楼，而且數字不對應，所以我要用文獻中列出的某些基因的小提琴圖確定我的聚類
根據文獻對應自己數據聚類
原文獻中有所有聚類的特異基因白对，所以我根據列出的基因去匹配我的聚類結果
##==============================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)
修改之后的聚類結果
一些基因的小提琴圖對應效果
結語
對于這次的數據重復盟步，基本符合預期結果，和文章的結果有點差距躏结，需要自己進一步研究問題出在哪里，下一次將繼續(xù)這篇文獻的數據復現狰域，主要是偽時間分析媳拴，目前的數據與代碼我已上傳github
轉載請注明：周小釗的博客>>>單細胞實戰(zhàn)(5)：復現文章中的聚類圖（1）
