转录组和单细胞下游基于R的数据分析-01

生信技能树

发布于 2024-03-25 15:32:12

27300

代码可运行

文章被收录于专栏：生信技能树生信技能树

运行总次数：0

代码可运行

单细胞转录组数据情况

数据链接是：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE212199

数据情况是：有两个分组Mock以及WNV

#samples
GSM6513450 Female Mock F1
GSM6513451 Male Mock M1
GSM6513452 Female WNV F2
GSM6513453 Male WNV M2

#数据详情
GSM6513450_F1_barcodes.tsv.gz 5.0 Kb
GSM6513450_F1_features.tsv.gz 272.8 Kb
GSM6513450_F1_matrix.mtx.gz 5.2 Mb
GSM6513451_M1_barcodes.tsv.gz 7.5 Kb
GSM6513451_M1_features.tsv.gz 272.8 Kb
GSM6513451_M1_matrix.mtx.gz 8.2 Mb
GSM6513452_F2_barcodes.tsv.gz 6.9 Kb
GSM6513452_F2_features.tsv.gz 272.8 Kb
GSM6513452_F2_matrix.mtx.gz 7.1 Mb
GSM6513453_M2_barcodes.tsv.gz 7.6 Kb
GSM6513453_M2_features.tsv.gz 272.8 Kb
GSM6513453_M2_matrix.mtx.gz 8.2 Mb

数据下载整理以及读取

所有的代码及运行脚本均来源于生信技能树！

提供的是10X格式的标准三个文件，选择下载数据之后需要对数据进行整理，将三个文件分别整理到对应的文件夹中。

#整理文件
fs=list.files('./','features')
fs

samples1= gsub('.tsv.gz','',gsub('features.','',fs))
samples1

samples2 = samples1

lapply(1:length(samples2), function(i){
  x=samples2[i]
  y=fs[i]
  dir.create(x,recursive = T)
  file.copy(from= y ,
            to=file.path(x,  'features.tsv.gz' )) 
  file.copy(from= gsub('features','matrix',gsub('tsv','mtx',y)),
            to= file.path(x, 'matrix.mtx.gz' ) ) 
  file.copy(from= gsub('features','barcodes',y),
            to= file.path(x, 'barcodes.tsv.gz' )) 
  
})

整理前

整理后

加载需要的R包，然后使用Read10X()函数将数据读取进来，然后创建seurta对象，即可进行后续的降维聚类分群。

#指定数据存放位置
samples=list.files("./GSE212199_RAW/outputs/")
samples
dir <- file.path('./GSE212199_RAW/outputs/',samples)
names(dir) <- samples
#读取数据创建Seurat对象
counts <- Read10X(data.dir = dir)

sce.all = CreateSeuratObject(counts,
                             min.cells = 5,
                             min.features = 300 )

dim(sce.all)   #查看基因数和细胞总数
as.data.frame(sce.all@assays$RNA$counts[1:10, 1:2])
table(sce.all@meta.data$orig.ident)  #查看每个样本的细胞数
head(sce.all@meta.data)

# 整理分组信息，通常是隐含在文件名，样品名字里面
sce.all$group=toupper( substring(colnames(sce.all),12,13))

table(sce.all@meta.data$group)
table(sce.all@meta.data$orig.ident)

质控

所有的代码及运行脚本均来源于生信技能树！

质控——去除掉低质量的基因与细胞

#step2:QC质控(运行命令)
dir.create("./1-QC")
setwd("./1-QC")
# 如果过滤的太狠，就需要去修改这个过滤代码
source('../scRNA_scripts/qc.R')
sce.all.filt = basic_qc(sce.all)
print(dim(sce.all))
print(dim(sce.all.filt))
setwd('../')
#指明物种
sp='mouse'

basic_qc函数

过滤结果:确定一下过滤前后的数量差异，如果相差太多要修改质控代码！

> print(dim(sce.all))
[1] 15968  4775
> print(dim(sce.all.filt))
[1] 15968  4568

质控前

质控后

harmony整合多个单细胞样品

所有的代码及运行脚本均来源于生信技能树！

因为GSE212199有四个样品，两个分组，咱们就需要使用Harmony整合一下

###### step3: harmony整合多个单细胞样品 ######
dir.create("2-harmony")
getwd()
setwd("2-harmony")
source('../scRNA_scripts/harmony.R')
# 默认 ScaleData 没有添加"nCount_RNA", "nFeature_RNA"
# 默认的
sce.all.int = run_harmony(sce.all.filt)
setwd('../')

harmony整合

harmony整合结果：可以根据结果去选择后续的resolution

Tree_diff_resolution

resolution_low

resolution_high

降维聚类分群

所有的代码及运行脚本均来源于生信技能树！

降维的resolution一般我们是选择0.1以及0.8，但是这次根据文章里面的结果图，所以还选择了0.2的分辨率

###### step4: 降维聚类分群和看标记基因库 ######
#原则上分辨率是需要自己肉眼判断，取决于个人经验。为了省力，我们直接看0.1和0.8即可
table(Idents(sce.all.int))
table(sce.all.int$seurat_clusters)
table(sce.all.int$RNA_snn_res.0.1)
table(sce.all.int$RNA_snn_res.0.2) 
table(sce.all.int$RNA_snn_res.0.8) 

getwd()
dir.create('check-by-0.1')
setwd('check-by-0.1')
sel.clust = "RNA_snn_res.0.1"
sce.all.int <- SetIdent(sce.all.int, value = sel.clust)
table(sce.all.int@active.ident) 

source('../scRNA_scripts/check-all-markers.R')
setwd('../') 
getwd()

dir.create('check-by-0.2')
setwd('check-by-0.2')
sel.clust = "RNA_snn_res.0.2"
sce.all.int <- SetIdent(sce.all.int, value = sel.clust)
table(sce.all.int@active.ident) 
source('../scRNA_scripts/check-all-markers.R')
setwd('../') 
getwd()

dir.create('check-by-0.8')
setwd('check-by-0.8')
sel.clust = "RNA_snn_res.0.8"
sce.all.int <- SetIdent(sce.all.int, value = sel.clust)
table(sce.all.int@active.ident) 
source('../scRNA_scripts/check-all-markers.R')
setwd('../') 
getwd()

last_markers_to_check

那就直接选择0.2进行后续的命名吧！

亚群命名

所有的代码及运行脚本均来源于生信技能树！

###### step5: 确定单细胞亚群生物学名字 ######
# 一般来说，为了节省工作量，我们选择0.1的分辨率进行命名
# 因为命名这个步骤是纯人工 操作
# 除非0.1确实分群太粗狂了，我们就选择0.8 
source('scRNA_scripts/lib.R')
#sce.all.int = readRDS('2-harmony/sce.all_int.rds')
colnames(sce.all.int@meta.data)
DimPlot(sce.all.int,split.by = "orig.ident"  )
DimPlot(sce.all.int,group.by = "group"  )
FeaturePlot(sce.all.int,'Apoe',split.by  = "group" )
FeaturePlot(sce.all.int,'Cd74',split.by  = "group" )


astrocytes = c("AQP4", "ADGRV1", "GPC5", "RYR3") 
endothelial = c("CLDN5", "ABCB1", "EBF1") 
excitatory = c("CAMK2A", "CBLN2", "LDB2") 
inhibitory = c("GAD1", "LHFPL3", "PCDH15") 
microglia = c("C3", "LRMDA", "DOCK8") 
oligodendrocytes = c("MBP", "PLP1", "ST18") 
OPC='Tnr,Igsf21,Neu4,Gpr17'
Ependymal='Cfap126,Fam183b,Tmem212,pifo,Tekt1,Dnah12'
pericyte=c(  'DCN', 'LUM',  'GSN' ,'FGF7','MME', 'ACTA2','RGS5')

gene_list = list(
  Astro = astrocytes,
  Endo = endothelial,
  Excit = excitatory,
  Inhib = inhibitory,
  Mic = microglia,
  Oligo = oligodendrocytes,
  OPC= str_to_upper(trimws(strsplit(OPC,',')[[1]])),
  Ependymal= str_to_upper(trimws(strsplit(Ependymal,',')[[1]])) ,
  peri = pericyte
)
gene_list = lapply(gene_list,function(x){
  unique( c(str_to_title(x),str_to_upper(x)))
} )
gene_list
p3=DotPlot(sce.all.int, assay = "RNA", features = gene_list, 
           group.by = 'RNA_snn_res.0.2'
) + theme(axis.text.x = element_text(angle = 45, 
                                     vjust = 0.5, hjust=0.5))
p3
ggsave('makers_brain_RNA_snn_res.0.2.pdf',width = 9)

table(sce.all.int$RNA_snn_res.0.2)

# 付费环节 800 元人民币
if(T){
  sce.all.int
  celltype=data.frame(ClusterID=0:9,
                      celltype= 0:9) 
  #定义细胞亚群        
  celltype[celltype$ClusterID %in% c( 3 ),2]='CD4+ Tcells'
  celltype[celltype$ClusterID %in% c( 2),2]='CD8+ Tcells'
  celltype[celltype$ClusterID %in% c(  4 ),2]='Myeloids'  
  
  celltype[celltype$ClusterID %in% c(  0,1,5,7 ),2]='microglia'  
  celltype[celltype$ClusterID %in% c(  8 ),2]='BCells'  
  celltype[celltype$ClusterID %in% c(  9 ),2]='pericyte'  
  
  celltype[celltype$ClusterID %in% c( 6 ),2]='oligodendrocytes'  
   
  head(celltype)
  celltype
  table(celltype$celltype)
  sce.all.int@meta.data$celltype = "NA"
  
  for(i in 1:nrow(celltype)){
    sce.all.int@meta.data[which(sce.all.int@meta.data$RNA_snn_res.0.2 == celltype$ClusterID[i]),'celltype'] <- celltype$celltype[i]}
  Idents(sce.all.int)=sce.all.int$celltype
  
  sel.clust = "celltype"
  sce.all.int <- SetIdent(sce.all.int, value = sel.clust)
  table(sce.all.int@active.ident) 
  
  dir.create('check-by-celltype')
  setwd('check-by-celltype')
  source('../scRNA_scripts/check-all-markers.R')
  setwd('../') 
  getwd()
}

phe=sce.all.int@meta.data
save(phe,file = 'phe.Rdata')