单细胞测序—GSE218208(流程简化)

untar("GSE218208_RAW.tar")
rm(list = ls())
a = data.table::fread("GSM6736629_10x-PBMC-1_ds0.1974_CountMatrix.tsv.gz",data.table = F)
#fread无法直接设置行名
a[1:4,1:4]
##                 alias:gene AAACCCAAGTAGGGTC AAACCCACACCATTCC AAACCCATCTACACTT
## 1   TSPAN6:ENSG00000000003                0                0                0
## 2     DPM1:ENSG00000000419                0                1                0
## 3    SCYL3:ENSG00000000457                0                0                0
## 4 C1orf112:ENSG00000000460                0                0                0
library(tidyverse)
#只保留基因名称，并去重
a$`alias:gene` = str_split(a$`alias:gene`,":",simplify = T)[,1]
#str_split_i(a$`alias:gene`,":",i = 1)
a = distinct(a,`alias:gene`,.keep_all = T)
#把基因的名字设置为行名
a = column_to_rownames(a,var = "alias:gene")
a[1:4,1:4]
##          AAACCCAAGTAGGGTC AAACCCACACCATTCC AAACCCATCTACACTT AAACGAAAGCACGTCC
## TSPAN6                  0                0                0                0
## DPM1                    0                1                0                0
## SCYL3                   0                0                0                0
## C1orf112                0                0                0                0
library(Seurat)
#CreateSeuratObject允许直接接受一个矩阵
pbmc <- CreateSeuratObject(counts = a, 
                           project = "a", 
                           min.cells = 3, 
                           min.features = 200)

pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
head(pbmc@meta.data, 3)
##                  orig.ident nCount_RNA nFeature_RNA percent.mt
## AAACCCAAGTAGGGTC          a      10768         3213   7.030089
## AAACCCACACCATTCC          a       4102         1676   5.046319
## AAACCCATCTACACTT          a       4694         1740   6.305922
VlnPlot(pbmc, 
        features = c("nFeature_RNA",
                     "nCount_RNA", 
                     "percent.mt"), 
        ncol = 3,pt.size = 0.5)

pbmc = subset(pbmc,nFeature_RNA < 4200 &
                nCount_RNA < 18000 &
                percent.mt < 18)

f = "obj.Rdata"
if(!file.exists(f)){
  pbmc = pbmc %>% 
  NormalizeData() %>%  
  FindVariableFeatures() %>%  
  ScaleData(features = rownames(.)) %>%  
  RunPCA(features = VariableFeatures(.))  %>%
  FindNeighbors(dims = 1:15) %>% 
  FindClusters(resolution = 0.5) %>% 
  RunUMAP(dims = 1:15) %>% 
  RunTSNE(dims = 1:15)
  save(pbmc,file = f)
}

load(f)
ElbowPlot(pbmc)

p1 <- DimPlot(pbmc, reduction = "umap",label = T)+NoLegend();p1

library(dplyr)
f = "markers.Rdata"
if(!file.exists(f)){
  pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE,min.pct = 0.25)
  save(pbmc.markers,file = f)
}
load(f)
mks = pbmc.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_log2FC)
#如果一个基因同时是两簇的marker基因的话下面的代码会报错，因此需要去重
g = unique(mks$gene)
#可视化
DoHeatmap(pbmc, features = g) + NoLegend()+
  scale_fill_gradientn(colors = c("#2fa1dd", "white", "#f87669"))
DotPlot(pbmc, features = g,cols = "RdYlBu") +
  RotatedAxis()
#查看感兴趣的基因
VlnPlot(pbmc, features = g[1:3])
FeaturePlot(pbmc, features = g[1:4])

a = read.delim("../supp/markers.txt",header = F)
gt = split(a[,2],a[,1])

DotPlot(pbmc, features = gt,cols = "RdYlBu") +
  RotatedAxis()

#打出以下内容，复制到新建的anno.tx文档里
#照着DotPlot图用unique(a$V1)的输出结果补全细胞类型
writeLines(paste0(0:11,","))
## 0,
## 1,
## 2,
## 3,
## 4,
## 5,
## 6,
## 7,
## 8,
## 9,
## 10,
## 11,
celltype = read.table("anno.txt",sep = ",") #自己照着DotPlot图填的
celltype
##    V1           V2
## 1   0  Naive CD4 T
## 2   1   CD14+ Mono
## 3   2            B
## 4   3        CD8 T
## 5   4           NK
## 6   5  Naive CD4 T
## 7   6  Naive CD4 T
## 8   7 FCGR3A+ Mono
## 9   8   CD14+ Mono
## 10  9     Platelet
## 11 10           DC
## 12 11           DC
new.cluster.ids <- celltype$V2
names(new.cluster.ids) <- levels(pbmc)
seu.obj <- RenameIdents(pbmc, new.cluster.ids)
save(seu.obj,file = "seu.obj.Rdata")
p2 <- DimPlot(seu.obj, 
        reduction = "umap", 
        label = TRUE, 
        pt.size = 0.5) + NoLegend()
p2

library(celldex)
library(SingleR)
ls("package:celldex")
## [1] "BlueprintEncodeData"              "DatabaseImmuneCellExpressionData"
## [3] "HumanPrimaryCellAtlasData"        "ImmGenData"                      
## [5] "MonacoImmuneData"                 "MouseRNAseqData"                 
## [7] "NovershternHematopoieticData"
f = "../supp/single_ref/ref_BlueprintEncode.RData"
if(!file.exists(f)){
  ref <- celldex::BlueprintEncodeData()
  save(ref,file = f)
}
#ref兼容不同的变量名称
ref <- get(load(f))
library(BiocParallel)
scRNA = pbmc
test = scRNA@assays$RNA@layers$data
rownames(test) = Features(scRNA)
colnames(test) = Cells(scRNA)
pred.scRNA <- SingleR(test = test, 
                      ref = ref,
                      labels = ref$label.main, 
                      clusters = scRNA@active.ident)
# 查看注释成的哪些类型
pred.scRNA$pruned.labels
##  [1] "CD4+ T-cells" "Monocytes"    "B-cells"      "CD8+ T-cells" "NK cells"    
##  [6] "CD4+ T-cells" "CD4+ T-cells" "Monocytes"    "Monocytes"    "Monocytes"   
## [11] "DC"           "Monocytes"
#查看注释的准确性，每一列最黄的是注释结果
plotScoreHeatmap(pred.scRNA, clusters=pred.scRNA@rownames, fontsize.row = 9,show_colnames = T)
new.cluster.ids <- pred.scRNA$pruned.labels
names(new.cluster.ids) <- levels(scRNA)
levels(scRNA)
##  [1] "0"  "1"  "2"  "3"  "4"  "5"  "6"  "7"  "8"  "9"  "10" "11"
scRNA <- RenameIdents(scRNA,new.cluster.ids)
levels(scRNA)
## [1] "CD4+ T-cells" "Monocytes"    "B-cells"      "CD8+ T-cells" "NK cells"    
## [6] "DC"
p3 <- DimPlot(scRNA, reduction = "umap",label = T,pt.size = 0.5) + NoLegend()
#对比下手动注释与自动注释
p2+p3