-
Notifications
You must be signed in to change notification settings - Fork 0
/
single_cell_cluster.r
75 lines (64 loc) · 3.23 KB
/
single_cell_cluster.r
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
*Scripts : single_cell_cluster.r*
`rm(list = ls())
setwd("F:/R_wenjian")
library(Seurat)
library(data.table)
library(stringr)
library(dplyr)
library(rtracklayer)
library(ggplot2)
counts_d2 <- as.data.frame(fread("hebing_count.txt"))
gtf <- import("gencode.vM21.chr_patch_hapl_scaff.annotation.gtf.gz")
protein <- unique(gtf$gene_name[gtf$transcript_type %in% "protein_coding"])
mito <- unique(gtf$gene_name[as.character(seqnames(gtf)) %in% "chrM"])
genes_name <- unique(c(protein, mito))
rownames(counts_d2) <- counts_d2$Geneid
counts_d2 <- counts_d2[, c(7:93)]
#colnames(counts_d2) <- str_split_fixed(colnames(counts_d2), "/", 2)[, 1]
colnames(counts_d2)<- rownames(metadata2)
counts_d2<- counts_d2[genes_name, ]
klf5_d2 <- CreateSeuratObject(counts = counts_d2, project = "klf5_day2", min.cells = 3, min.features = 200)
klf5_d2[["percent.mt"]] <- PercentageFeatureSet(klf5_d2, pattern = "^mt-")
VlnPlot(klf5_d2, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
klf5_d2 <- subset(klf5_d2, subset = nFeature_RNA > 1000 & nCount_RNA >200 & percent.mt < 10)
klf5_d2 <- NormalizeData(klf5_d2, normalization.method = "LogNormalize", scale.factor = 10000)
klf5_d2 <- FindVariableFeatures(klf5_d2, selection.method = "vst", nfeatures = 3000)
#top15 <- head(VariableFeatures(klf5_d2), 15)
#plot1 <- VariableFeaturePlot(klf5_d2)
#LabelPoints(plot = plot1, points = top15, repel = TRUE)
all.genes <- rownames(klf5_d2)
klf5_d2 <- ScaleData(klf5_d2, features = all.genes)
klf5_d2 <- RunPCA(klf5_d2, features = VariableFeatures(object = klf5_d2))
klf5_d2 <- JackStraw(klf5_d2, num.replicate = 100)
klf5_d2 <- ScoreJackStraw(klf5_d2, dims = 1:20)
klf5_d2 <- FindNeighbors(klf5_d2, dims = 1:10)
klf5_d2 <- FindClusters(klf5_d2, resolution = 1)
klf5_d2 <- RunTSNE(klf5_d2, dims.use = 1:10, perplexity = 5)
TSNEPlot(object = klf5_d2,pt.size = 2)
#cluster1.markers <- FindMarkers(object = klf5_d2, ident.1 = 3, min.pct = 0.25)
#VlnPlot(object = klf5_d2, features = c("Sox9", "Krt19"))
klf5_d2.markers <- FindAllMarkers(klf5_d2, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
write.table(klf5_d2.markers,file = "fenqun_marker3.tsv",row.names = F, sep="\t", quote=F)
###picture
klf5_d2.markers %>% group_by(cluster) %>% top_n(n = 3, wt = avg_logFC)
VlnPlot(object = klf5_d2, features = c("Ttn", "D10Wsu102e","Sfrp5",
"Cldn2", "Dmbt1", "Ly6d",
"Lurap1l", "Apoe", "Anxa5",
"Vcam1", "Slc34a2", "Hba-a1"), ncol = 3)
FeaturePlot(object = klf5_d2, features = c("Ttn", "D10Wsu102e","Sfrp5",
"Cldn2", "Dmbt1", "Ly6d",
"Lurap1l", "Apoe", "Anxa5",
"Vcam1", "Slc34a2", "Hba-a1"), ncol = 3)
top50 <- klf5_d2.markers %>% group_by(cluster) %>% top_n(n = 50, wt = avg_logFC)
DoHeatmap(klf5_d2, features = top50$gene, size = 2)`
#### change heatmap's color
DoHeatmap(clu_qubatch_hoxb_sig_2, features = top15$gene, size = 6)+
theme(
axis.text.y = element_text(size = 10,
family = "Times",face = "bold",colour = "black"))+
scale_fill_gradient2(
low = "blue",
mid = "white",
high ="red",
name = "Expression"
)