区区20万个单细胞居然超内存了

a=fread('inputs/GSE131907_Lung_Cancer_raw_UMI_matrix.txt.gz',
        data.table = F)
a[1:4,1:4]
rownames(a)=a[,1]
a=a[,-1]
# 208,506 cells derived from 58 lung adenocarcinomas from 44 patients, 

# 11 tumour, 11 distant normal lung, 10 normal lymph node, and 10 metastatic brain tissue samples 
# were obtained from LUAD patients conserving surgery without prior treatment. 
#  7 metastatic lymph node and 4 lung tumour tissue samples from advanced stage LUAD patient
# endobronchial ultrasound (EBUS) and bronchoscopy (BRONCHO).

kp1 = grepl('LN',colnames(a))
table(kp1)
kp2 = grepl('T',colnames(a))
table(kp2 )

a[1:4,1:4]
dim(a)

if(F){
  format(object.size(a) ,units='Gb')
  A <- as( a , "sparseMatrix") 
  B <- Matrix(as.matrix(a), sparse = TRUE) 
  library(Matrix)
  library(mltools)
  A = mltools::sparsify(a)
  format(object.size(A) ,units='Gb')
}

sampl = paste0(phe[,2],phe[,3])
 
sceList = lapply(unique(sampl) ,function(pro){ 
  # pro=unique(sampl)[1] 
  print(pro)  
  sce =CreateSeuratObject(counts = a[,sampl %in% pro],
                          project =  gsub( '_matrix','',pro)  ,
                          min.cells = 5,
                          min.features = 500 )
  
  return(sce)
})
names(sceList) 

sce.all=merge(x=sceList[[1]],
              y=sceList[ -1 ] )
as.data.frame(sce.all@assays$RNA@counts[1:10, 1:2])
head(sce.all@meta.data, 10)
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='human'

###### 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('../')

###### 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.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.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()

  celltype=data.frame(ClusterID=0:13 ,
                      celltype= 0:13) 
  #定义细胞亚群      
  celltype[celltype$ClusterID %in% c( 0, 11,13 ),2]='Tcells'
  celltype[celltype$ClusterID %in% c( 6 ),2]='Bcells'
  celltype[celltype$ClusterID %in% c( 5 ),2]='mast'
  
  celltype[celltype$ClusterID %in% c(  3 ),2]='cycle'
  celltype[celltype$ClusterID %in% c(  4 ),2]='fibro'
  celltype[celltype$ClusterID %in% c(  7 ),2]='endo'
  celltype[celltype$ClusterID %in% c(  9 ),2]='pDC'
  
  celltype[celltype$ClusterID %in% c(  1 ),2]='myeloids'
  celltype[celltype$ClusterID %in% c( 2,12,8 ),2]='epi'