当前位置:网站首页>/silicosis/geo/GSE184854_scRNA-seq_mouse_lung_ccr2/GSE184854_RAW/GSM5598265_matrix_inflection_demult
/silicosis/geo/GSE184854_scRNA-seq_mouse_lung_ccr2/GSE184854_RAW/GSM5598265_matrix_inflection_demult
2022-07-02 03:06:00 【youngleeyoung】
step1
library(Seurat)
#https://www.jianshu.com/p/5b26d7bc37b7
getwd()
path="G:/silicosis/geo/GSE184854_scRNA-seq_mouse_lung_ccr2/GSE184854_RAW/"
setwd(path)
getwd()
#new_counts=read.table(file = "G:/silicosis/geo/GSE184854_scRNA-seq_mouse_lung_ccr2/GSE184854_RAW/GSM5598265_matrix_inflection_demulti_SilicaWT.txt/GSM5598265_matrix_inflection_demulti_SilicaWT.txt")
head(new_counts)
#或者可以直接 create色urat object(counts=counts)
mydata <- CreateSeuratObject(counts = read.table(file = "G:/silicosis/geo/GSE184854_scRNA-seq_mouse_lung_ccr2/GSE184854_RAW/GSM5598265_matrix_inflection_demulti_SilicaWT.txt/GSM5598265_matrix_inflection_demulti_SilicaWT.txt"),
min.cells = 10, project = "mydata_scRNAseq")
dim(mydata)
Idents(mydata)
mydata$orig_stim=Idents(mydata)
table(Idents(mydata))
colnames(mydata)
library(stringr)
teststring=colnames(mydata)[1:50]
str_split(teststring,"_",simplify = T)
str_split(teststring,"_",simplify = T)[,2]
table(str_split(colnames(mydata),"_",simplify = T)[,2])
table(str_split(colnames(mydata),"_",simplify = T)[,1])
mydata$percent.mt=PercentageFeatureSet(mydata,pattern = "^Mt")
Idents(mydata)
grepl(colnames(mydata),pattern = 'A0301')
colnames(mydata)[grepl(colnames(mydata),pattern = 'A0301')]
table(Idents(mydata))
mydata=RenameIdents(mydata,'A0301'='day_21_1','A0302'='day_21_2','A0303'='day_21_3',
'A0304'='day_7_1','A0305'='day_7_2','A0306'='day_7_3',
'A0307'='day_3_1','A0308'='day_3_2','A0309'='day_3_3',
'doublet'='doublet')
mydata$mystim=Idents(mydata)
table(mydata$mystim)
dim(mydata)
Idents(mydata)=mydata$orig_stim
mydata=RenameIdents(mydata,'A0301'='Day_21','A0302'='Day_21','A0303'='Day_21',
'A0304'='Day_7','A0305'='Day_7','A0306'='Day_7',
'A0307'='Day_3','A0308'='Day_3','A0309'='Day_3',
'doublet'='Doublet')
mydata$my_stim_3=Idents(mydata)
save(mydata,file = "G:/silicosis/geo/GSE184854_scRNA-seq_mouse_lung_ccr2/GSE184854_RAW/mydata.rds")
step2
library(dplyr)
library(cowplot)
library(Seurat)
library(harmony)
getwd()
path="G:/silicosis/geo/GSE184854_scRNA-seq_mouse_lung_ccr2/GSE184854_RAW_all_merged/step_2_harmony/"
dir.create(path)
setwd(path)
getwd()
#load("G:/silicosis/geo/GSE184854_scRNA-seq_mouse_lung_ccr2/GSE184854_RAW/mydata.rds")
load("G:/silicosis/geo/GSE184854_scRNA-seq_mouse_lung_ccr2/GSE184854_RAW_all_merged/merged.rds")
table(Idents(All.merge))
All=subset(All.merge,idents = c('A0301', 'A0302', 'A0303',
'A0304', 'A0305' ,'A0306',
'A0307' , 'A0308', 'A0309'))
dim(All)
table(Idents(All))
#All$percent.mt=PercentageFeatureSet(All,pattern = "^mt-")
pdf("1_contorl_QC.pdf")
VlnPlot(All, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
dev.off()
All <- subset(All, subset = nFeature_RNA > 200 & percent.mt < 20)## 15499,16133;
All = All%>%Seurat::NormalizeData(verbose = FALSE) %>%
FindVariableFeatures(selection.method = "vst", nfeatures = 2000) %>%
ScaleData(verbose = FALSE)
All = RunPCA(All, npcs = 50, verbose = FALSE)
pdf("2_ElbowPlot.pdf")
ElbowPlot(All, ndims = 50)
dev.off()
table(Idents(All))
All$stim=Idents(All)
#All@meta.data$stim <- c(rep("case", length(grep("1$", All@assays$RNA@counts@Dimnames[[2]]))), rep("ctrl", length(grep("2$", All@assays$RNA@counts@Dimnames[[2]])))) ## 8186,7947;
pdf("2_pre_harmony_harmony_plot.pdf")
options(repr.plot.height = 5, repr.plot.width = 12)
p1 <- DimPlot(object = All, reduction = "pca", pt.size = .1, group.by = "stim")
p2 <- VlnPlot(object = All, features = "PC_1", group.by = "stim", pt.size = .1)
plot_grid(p1, p2)
dev.off()
##########################run harmony
All <- All %>% RunHarmony("stim", plot_convergence = TRUE)
harmony_embeddings <- Embeddings(All, 'harmony')
pdf("1_contorl_QC_.pdf")
VlnPlot(All, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
dev.off()
pdf("2_after_harmony_harmony_plot.pdf")
options(repr.plot.height = 5, repr.plot.width = 12)
p3 <- DimPlot(object = All, reduction = "harmony", pt.size = .1, group.by = "stim")
p4 <- VlnPlot(object = All, features = "harmony_1", group.by = "stim", pt.size = .1)
plot_grid(p3, p4)
dev.off()
#######################cluster
All <- All %>%
RunUMAP(reduction = "harmony", dims = 1:30) %>%
RunTSNE(reduction = "harmony", dims = 1:30) %>%
FindNeighbors(reduction = "harmony", dims = 1:30)
All<-All%>% FindClusters(resolution = 3) %>% identity()
options(repr.plot.height = 4, repr.plot.width = 10)
pdf("3_after_harmony_umap_two_group.pdf")
DimPlot(All, reduction = "umap", group.by = "stim", pt.size = .1)
dev.off()
pdf("3_after_harmony_cluster_UMAP.pdf")
DimPlot(All, reduction = "umap", label = TRUE, pt.size = .1)
dev.off()
pdf("3_umap_samples_split.pdf")
DimPlot(All, reduction = "umap", pt.size = .1, split.by = "stim", label = T)
dev.off()
pdf("3_after_harmony_tsne_two_group.pdf")
DimPlot(All, reduction = "tsne", group.by = "stim", pt.size = .1)
dev.off()
pdf("3_after_harmony_cluster_tSNE.pdf")
DimPlot(All, reduction = "tsne", label = TRUE, pt.size = .1)
dev.off()
pdf("3_tSNE_samples_split.pdf")
DimPlot(All, reduction = "tsne", pt.size = .1, split.by = "stim", label = T)
dev.off()
#######################################################################
stat = as.matrix(table(Idents(All), All$stim))
write.table(stat, "cluster_stat.txt", sep = "\t", quote = F, col.names = T, row.names = T)
################################################################
All<-All%>% FindClusters(resolution = 2) %>% identity()
Disease.markers <- FindAllMarkers(All, min.pct = 0.25, logfc.threshold = 0.35, only.pos = T)
head(Disease.markers)
top100markers <- Disease.markers %>% group_by(cluster) %>% slice_max(avg_log2FC,n=100)
openxlsx::write.xlsx(top100markers,'top100markers_for_all_res2.xlsx')
#write.table(top20markers, "top20_markers.txt", sep = "\t", quote = F, col.names = T, row.names = F)
#save(All, file = "G:/silicosis/geo/GSE184854_scRNA-seq_mouse_lung_ccr2/GSE184854_RAW_all_merged/step_2_harmony/all_harmony.rds")
load("G:/silicosis/geo/GSE184854_scRNA-seq_mouse_lung_ccr2/GSE184854_RAW_all_merged/step_2_harmony/all_harmony.rds")
getwd()
DimPlot(All,label = T)
FeaturePlot(All,features = c('Cd68'))
FeaturePlot(All,features = c('Lyz2'))
FeaturePlot(All,features = c('Mrc1'))
FeaturePlot(All,features = c('Adgre1'))
FeaturePlot(All,features = c('Axl'))
FeaturePlot(All,features = c('Csf1r'))
FeaturePlot(All,features = c('Ly6c2'))
FeaturePlot(All,features = c('Il7r'))
FeaturePlot(All,features = c('Car4'))
FeaturePlot(All,features = c('Siglecf'))
FeaturePlot(All,features = c('Mki67'))
FeaturePlot(All,features = c('Jchain','Tnfrsf17'),split.by = 'mystim')
DotPlot(All,features = c('Jchain','Tnfrsf17'),split.by = 'mystim')+RotatedAxis()+ggplot2::coord_flip()
colnames(All)
grepl(colnames(All),pattern = 'WT')
table(grepl(colnames(All),pattern = 'WT'))
table(grepl(colnames(All),pattern = 'CCR'))
All$orig_stim=Idents(All)
if(1==1){
All$mystim_validate=ifelse(grepl(colnames(All),pattern = 'WT') &
(grepl(colnames(All),pattern = 'A0301')|grepl(colnames(All),pattern = 'A0302')|grepl(colnames(All),pattern = 'A0303')),
paste0(colnames(All),'_Day_21_WT'),
ifelse(grepl(colnames(All),pattern = 'WT') &
(grepl(colnames(All),pattern = 'A0304')|grepl(colnames(All),pattern = 'A0305')|grepl(colnames(All),pattern = 'A0306'))
,paste0(colnames(All),'_Day_7_WT'),
ifelse(grepl(colnames(All),pattern = 'WT') &
(grepl(colnames(All),pattern = 'A0307')|grepl(colnames(All),pattern = 'A0308')|grepl(colnames(All),pattern = 'A0309'))
,paste0(colnames(All),'_Day_3_WT'),
ifelse(grepl(colnames(All),pattern = 'CCR') &
(grepl(colnames(All),pattern = 'A0301')|grepl(colnames(All),pattern = 'A0302')|grepl(colnames(All),pattern = 'A0303'))
,paste0(colnames(All),'_Day_21_CCR'),
ifelse(grepl(colnames(All),pattern = 'CCR') &
(grepl(colnames(All),pattern = 'A0304')|grepl(colnames(All),pattern = 'A0305')|grepl(colnames(All),pattern = 'A0306'))
,paste0(colnames(All),'_Day_7_CCR'),
ifelse(grepl(colnames(All),pattern = 'CCR') &
(grepl(colnames(All),pattern = 'A0307'))
,paste0(colnames(All),'_Day_3_CCR'),paste0(colnames(All),'_Day_3_CCR'))
)
))
)
)
All$mystim=ifelse(grepl(colnames(All),pattern = 'WT') &
(grepl(colnames(All),pattern = 'A0301')|grepl(colnames(All),pattern = 'A0302')|grepl(colnames(All),pattern = 'A0303')),
paste0('_Day_21_WT'),
ifelse(grepl(colnames(All),pattern = 'WT') &
(grepl(colnames(All),pattern = 'A0304')|grepl(colnames(All),pattern = 'A0305')|grepl(colnames(All),pattern = 'A0306'))
,paste0('_Day_7_WT'),
ifelse(grepl(colnames(All),pattern = 'WT') &
(grepl(colnames(All),pattern = 'A0307')|grepl(colnames(All),pattern = 'A0308')|grepl(colnames(All),pattern = 'A0309'))
,paste0('_Day_3_WT'),
ifelse(grepl(colnames(All),pattern = 'CCR') &
(grepl(colnames(All),pattern = 'A0301')|grepl(colnames(All),pattern = 'A0302')|grepl(colnames(All),pattern = 'A0303'))
,paste0('_Day_21_CCR'),
ifelse(grepl(colnames(All),pattern = 'CCR') &
(grepl(colnames(All),pattern = 'A0304')|grepl(colnames(All),pattern = 'A0305')|grepl(colnames(All),pattern = 'A0306'))
,paste0('_Day_7_CCR'),
ifelse(grepl(colnames(All),pattern = 'CCR') &
(grepl(colnames(All),pattern = 'A0307'))
,paste0('_Day_3_CCR'),paste0('_Day_3_CCR'))
)
))
)
)
}
table(All$mystim)
table(All$orig_stim)
table(All$orig.ident)
table(grep(colnames(All),pattern = 'WT') & (grep(colnames(All),pattern = 'A0301')|
grep(colnames(All),pattern = 'A0302')|
grep(colnames(All),pattern = 'A0303'))
)
All$mystim=ifelse(grepl(colnames(All),pattern = 'WT'),'WT','CCR')
table(All$mystim)
table(ifelse(grepl(colnames(All),pattern = 'WT'),'WT','CCR'))
length(grep(colnames(All),pattern = 'WT'))
table(All$orig.ident)
Idents(All)=All$mystim
table(grepl(colnames(All),pattern = 'WT') )
grepl(c('af','faf','cd'),pattern = 'a')
step3
getwd()
path="G:/silicosis/geo/GSE184854_scRNA-seq_mouse_lung_ccr2/step_3_mono_macrophage/"
dir.create(path)
setwd(path)
getwd()
load("G:/silicosis/geo/GSE184854_scRNA-seq_mouse_lung_ccr2/step_2_harmony/silicosis_mouse_harmony.rds")
table(Idents(All))
#csf1r and cd68
macrophage_and_monocyte=subset(All,idents =c('1','7','11','12','15','22',
'31','35') )
table(Idents(macrophage_and_monocyte))
DimPlot(macrophage_and_monocyte)
macrophage_and_monocyte$orig_idents_from_All=Idents(macrophage_and_monocyte)
table(Idents(macrophage_and_monocyte),macrophage_and_monocyte$stim)
subset_data=macrophage_and_monocyte
subset_data = subset_data %>%
Seurat::NormalizeData(verbose = FALSE) %>%
FindVariableFeatures(selection.method = "vst", nfeatures = 2000) %>%
ScaleData(verbose = FALSE) %>%
RunPCA(npcs = 50, verbose = FALSE)
ElbowPlot(subset_data, ndims = 50)
VlnPlot(subset_data, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
##########################run harmony
#BiocManager::install('harmony')
library('harmony')
subset_data <- subset_data %>% RunHarmony("stim", plot_convergence = TRUE)
harmony_embeddings <- Embeddings(subset_data, 'harmony')
#######################cluster
dims = 1:30
subset_data <- subset_data %>%
RunUMAP(reduction = "harmony", dims = dims) %>%
RunTSNE(reduction = "harmony", dims = dims) %>%
FindNeighbors(reduction = "harmony", dims = dims)
DimPlot(subset_data)
save(subset_data, file="sepsis_Endothelial.rds")
file=getwd()
r = 0.2
load(paste(file, "sepsis_Endothelial.rds", sep="/"))
dir.create(paste(file,r,sep="/"))
setwd(paste(file,r,sep="/"))
subset_data <- FindClusters(subset_data, resolution = r)
save(subset_data, file=paste0("sepsis_Endothelial_r", r, ".rds"))
subset_data <- ScaleData(subset_data, verbose = FALSE, features = rownames(subset_data))
pdf("1_Endothelial_cluster_TSNE.pdf")
p = DimPlot(subset_data, reduction = "tsne", label = TRUE, pt.size = 1,label.size = 6)
print(p)
dev.off()
pdf("1_Endothelial_cluster_UMAP.pdf")
p = DimPlot(subset_data, reduction = "umap", label = TRUE, pt.size = 1,label.size = 6)
print(p)
dev.off()
stat = cbind(table(Idents(subset_data), subset_data$stim), All = rowSums(table(Idents(subset_data), subset_data$stim)))
stat = rbind(stat, All = colSums(stat))
stat<-as.data.frame(stat)
stat
openxlsx::write.xlsx(stat, "2_Endothelial_cluster_stat.xlsx", col.names=T, row.names=T)
######################## find marker
markers <- FindAllMarkers(subset_data, min.pct = 0.25, logfc.threshold = 0.25, only.pos=T)
openxlsx::write.xlsx(markers,"3_Endothelial_cluster_markers.xlsx", col.names=T,row.names=F)
FeaturePlot(subset_data,features = c('C1qa',
'C1qb',
'C1qc'))
FeaturePlot(subset_data,features = c('Mmp9',
'Spp1',
'Marcks',
'Il1rn'))
FeaturePlot(subset_data,features = c('Mmp9',
'Spp1',
'Marcks',
'Il1rn'),split.by = 'stim')
DotPlot(subset_data,features = c('Mmp9',
'Spp1',
'Marcks',
'Il1rn',
'C1qa','C1qb','C1qc'))
r=0.3
load(paste(file, "sepsis_Endothelial.rds", sep="/"))
dir.create(paste(file,r,sep="/"))
setwd(paste(file,r,sep="/"))
subset_data <- FindClusters(subset_data, resolution = r)
save(subset_data, file=paste0("sepsis_Endothelial_r", r, ".rds"))
subset_data <- ScaleData(subset_data, verbose = FALSE, features = rownames(subset_data))
pdf("1_Endothelial_cluster_TSNE.pdf")
p = DimPlot(subset_data, reduction = "tsne", label = TRUE, pt.size = 1,label.size = 6)
print(p)
dev.off()
pdf("1_Endothelial_cluster_UMAP.pdf")
p = DimPlot(subset_data, reduction = "umap", label = TRUE, pt.size = 1,label.size = 6)
print(p)
dev.off()
stat = cbind(table(Idents(subset_data), subset_data$stim), All = rowSums(table(Idents(subset_data), subset_data$stim)))
stat = rbind(stat, All = colSums(stat))
stat<-as.data.frame(stat)
stat
openxlsx::write.xlsx(stat, "2_Endothelial_cluster_stat.xlsx", col.names=T, row.names=T)
######################## find marker
markers <- FindAllMarkers(subset_data, min.pct = 0.25, logfc.threshold = 0.25, only.pos=T)
openxlsx::write.xlsx(markers,"3_Endothelial_cluster_markers.xlsx", col.names=T,row.names=F)
DimPlot(subset_data,label = T)
step4
getwd()
path="G:/silicosis/geo/GSE184854_scRNA-seq_mouse_lung_ccr2/step_4_macrophage_from_0.3"
dir.create(path)
setwd(path)
getwd()
load("G:/silicosis/geo/GSE184854_scRNA-seq_mouse_lung_ccr2/step_3_mono_macrophage/0.3/sepsis_Endothelial_r0.3.rds")
DimPlot(subset_data,label = T)
subset_data=subset(subset_data,idents = c('1','2','4'))
DimPlot(subset_data,label = T)
subset_data$cluster124=Idents(subset_data)
#合并
subset_data = subset_data %>%
Seurat::NormalizeData(verbose = FALSE) %>%
FindVariableFeatures(selection.method = "vst", nfeatures = 2000) %>%
ScaleData(verbose = FALSE) %>%
RunPCA(npcs = 50, verbose = FALSE)
ElbowPlot(subset_data, ndims = 50)
VlnPlot(subset_data, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
##########################run harmony
#BiocManager::install('harmony')
library('harmony')
subset_data <- subset_data %>% RunHarmony("stim", plot_convergence = TRUE)
harmony_embeddings <- Embeddings(subset_data, 'harmony')
#######################cluster
dims = 1:30
subset_data <- subset_data %>%
RunUMAP(reduction = "harmony", dims = dims) %>%
RunTSNE(reduction = "harmony", dims = dims) %>%
FindNeighbors(reduction = "harmony", dims = dims)
DimPlot(subset_data)
FeaturePlot(subset_data,features = c('C1qa',
'C1qb',
'C1qc'))
FeaturePlot(subset_data,features = c('Mmp9',
'Spp1',
'Marcks',
'Il1rn'))
FeaturePlot(subset_data,features = c('Mmp9',
'Spp1',
'Marcks',
'Il1rn'),split.by = 'stim')
DotPlot(subset_data,features = c('Mmp9',
'Spp1',
'Marcks',
'Il1rn',
'C1qa','C1qb','C1qc',
'Csf1','Csf1r',
'Car4', 'Ctsk', 'Chil3', 'S100a1','Wfdc21',
'Ear1', 'Fabp1',
'Itgam', 'Cd36','Gpnmb',
'Litaf', 'Jund', 'Bhlhe40', 'Bhlhe41','Klf9'))
subset_data=RenameIdents(subset_data,'4'='IM','2'='AM','1'='Mo-AM')
getwd()
#save(subset_data,file = "G:/silicosis/geo/GSE184854_scRNA-seq_mouse_lung_ccr2/step_4_macrophage_from_0.3/im_am_from_res0.3.rds")
load("G:/silicosis/geo/GSE184854_scRNA-seq_mouse_lung_ccr2/step_4_macrophage_from_0.3/im_am_from_res0.3.rds")
library(Seurat)
DotPlot(subset_data,features = c('Chi3l1', 'Marcks', 'Il1rn', 'Pla2g7', 'Mmp9', 'Spp1',
'Mmp12','Mmp14'))
DotPlot(All.merge,features = c('Spp1','Mmp14','Il1rn','Marcks','Pla2g7'))
VlnPlot(subset_data,features = c('Chi3l1', 'Marcks', 'Il1rn', 'Pla2g7', 'Mmp9', 'Spp1',
'Mmp12','Mmp14'))
#'Apoe', 'Mafb',
DotPlot(subset_data,features = c('C1qa','C1qb','C1qc',
'Ear1','Fabp1', 'Ctsk','Car4', 'Chil3', 'S100a1' , 'Wfdc21',
'Itgam', 'Litaf', 'Gpnmb','Apoe','Mafb',
'Mmp14','Spp1'))+RotatedAxis()+ggplot2:::coord_flip()
DotPlot(All.merge,features = c('Il34','Csf1'))
VlnPlot(All.merge,features = 'Il34')
VlnPlot(All.merge,features = 'Csf1')
VlnPlot(All.merge,features = c('Pdgfa','App','Tgfb1'))
DotPlot(All.merge,features = c('Il34','Csf1','Fgf2'))
边栏推荐
- Mathematical calculation in real mode addressing
- After marriage
- 设置状态栏颜色
- 表单自定义校验规则
- Principle of computer composition - interview questions for postgraduate entrance examination (review outline, key points and reference)
- 2022-2028 global deep sea generator controller industry research and trend analysis report
- Verilog 时序控制
- Unit · elementary C # learning notes
- buu_ re_ crackMe
- 多线程查询,效率翻倍
猜你喜欢
Après le mariage
Render header usage of El table
結婚後
MongoDB非关系型数据库
批量检测url是否存在cdn—高准确率
2022-2028 global nano abrasive industry research and trend analysis report
Formatting logic of SAP ui5 currency amount display
数据传输中的成帧
el-table的render-header用法
[staff] diacritical mark (ascending sign | descending sign B | double ascending sign x | double descending sign BB)
随机推荐
Ten minutes will take you in-depth understanding of multithreading - multithreaded teamwork: synchronous control
Gradle 笔记
【无标题】
QT实现界面跳转
[untitled]
Start a business
2022-2028 global soft capsule manufacturing machine industry research and trend analysis report
Mongodb base de données non relationnelle
Mathematical calculation in real mode addressing
寻找重复数[抽象二分/快慢指针/二进制枚举]
Render header usage of El table
Redis set command line operation (intersection, union and difference, random reading, etc.)
Connected block template and variants (4 questions in total)
Which brand of sports headset is better? Bluetooth headset suitable for sports
The number one malware in January 2022: lokibot returned to the list, and emotet returned to the top
[road of system analyst] collection of wrong topics in enterprise informatization chapter
[staff] pitch representation (bass clef | C1 36 note pitch representation | C2 48 note pitch representation | C3 60 note pitch representation)
Unit · elementary C # learning notes
PMP personal sprint preparation experience
Baohong industry | four basic knowledge necessary for personal finance