单细胞实战之inferCNV——入门到进阶(中级篇3）

rm(list = ls())
library(infercnv)
library(Seurat)
library(gplots)
library(ggplot2)
library(qs)
sce <- qread("./4-Corrective-data/sce.all.qs")

#检查一下自己导入进来的数据
DimPlot(sce,reduction = 'umap',
        label = TRUE,pt.size = 0.5) +NoLegend()

dir.create("11-infercnv")
setwd("./11-infercnv")

table(Idents(sce))
table(sce@meta.data$seurat_clusters)
table(sce@meta.data$orig.ident)

# 提取部分细胞
sce@meta.data
sce <- sce[,sce$celltype %in% c("plasma","fibroblasts","epithelial/cancer cells")]

# 文件制作1：表达量矩阵
# 除了用GetAssayData函数其实也可以直接sce@assays$RNA@count即可
dat <- GetAssayData(sce,layer = 'counts',assay = 'RNA')
dat[1:4,1:4]

# 文件制作2：样本的描述
groupinfo <- data.frame(v1 = colnames(dat),v2 = sce@meta.data$celltype)
head(groupinfo)
#                                  v1                      v2
# 1 GSM5688706_WGC_AACACACGTGCTTCAA-1 epithelial/cancer cells
# 2 GSM5688706_WGC_AACACACTCCCGTTCA-1 epithelial/cancer cells
# 3 GSM5688706_WGC_AACACACTCTGCTAGA-1                  plasma
# 4 GSM5688706_WGC_AACCATGAGGCACAAC-1                  plasma
# 5 GSM5688706_WGC_AACCATGCAAAGCTCT-1 epithelial/cancer cells
# 6 GSM5688706_WGC_AACCCAAGTGCCGTAC-1                  plasma

# 文件制作3：基因在染色体中的坐标
library(AnnoProbe)
#annoGene 函数返回一个数据框，包含输入基因的详细注释信息。
#注释信息可能包括基因名称、染色体位置、基因描述等。
geneInfor <- annoGene(rownames(dat),"SYMBOL","human") #物种需要小心哦
# 使用逻辑条件来移除包含 chrM, chrX, 和 chrY 的行
geneInfor <- geneInfor[!geneInfor$chr %in% c("chrM", "chrX", "chrY"), ]
#提取chr后后面的数字并转化为num，从而按这个num排序
#sub函数用于把chr替换为空
geneInfor$chr_num <- as.numeric(sub("chr", "", geneInfor$chr))
colnames(geneInfor)

#with函数的作用是简化写法，可问gpt
geneInfor <- geneInfor[with(geneInfor,order(chr_num,start)),c(1,4:6)]
geneInfor <- geneInfor[!duplicated(geneInfor[,1]),]

length(unique(geneInfor[,1]))
head(geneInfor)
#        SYMBOL  chr  start    end
# 23 AP006222.2 chr1 266855 268655
# 51     FAM87B chr1 817371 819837
# 53  LINC00115 chr1 826206 827522
# 55     FAM41C chr1 868071 876903
# 62     SAMD11 chr1 923928 944581
# 63      NOC2L chr1 944203 959309

#保留行名在geneInfor第一列中存在的行。
dat <- dat[rownames(dat) %in% geneInfor[,1],] 
#match(x, table)：match函数返回x中每个元素在table中的位置索引。
#获得位置后对dat进行重新排序使其跟geneInfor中的顺序一致
dat <- dat[match(geneInfor[,1],rownames(dat)),]
dim(dat)
head(groupinfo)
dat[1:4,1:4]
table(groupinfo$v2)
dim(groupinfo)

#为了节约计算机资源，直接抽样
#如果真实项目请不要用抽样
#kp <- sample(1:nrow(groupinfo),500)
#groupinfo <- groupinfo[kp,]
#dat <- dat[,kp]

# 保存矩阵文件
#expFile:是一个变量，存储写入的文件的文件名或路径。在这里文件名是expFile.txt。
expFile <- 'expFile.txt' #定义输出文件名
colnames(dat) <- gsub("-", "_", colnames(dat))
write.table(dat, file = expFile, sep = '\t', quote = F)
# 保存细胞注释文件
groupFiles <- 'groupFiles.txt'
groupinfo$v1 <- gsub("-", "_", groupinfo$v1)
write.table(groupinfo,file = groupFiles, sep = '\t',
            quote = F, col.names = F, row.names = F)
# 保存基因排序文件
head(geneInfor)
geneFile <- 'geneFile.txt'
write.table(geneInfor, file = geneFile, sep = '\t',
            quote = F, col.names = F, row.names = F)

#创建对象，请注意文件需要一一对应哦！
#ref_group_names 这里的细胞是正常对照，然后跟其他的细胞比较
infercnv_obj <- CreateInfercnvObject(raw_counts_matrix = expFile,
                                     annotations_file = groupFiles,
                                     delim = "\t",
                                     gene_order_file = geneFile,
                                     ref_group_names = c("plasma","fibroblasts") #这里用plasma和fibroblasts做参照
)

infercnv_obj <- infercnv::run(infercnv_obj,
                               cutoff =  0.1, #smart-seq选择1,10X选择0.1
                               out_dir = "infercnv_output", # dir is auto
                               # 是否根据细胞注释文件的分组
                               # 对肿瘤细胞进行分组
                               # 影响read.dendrogram, 如果有多个细胞类型，且设置为TRUE，
                               # 后续的read.dendrogram无法执行
                               cluster_by_groups =  TRUE, #是否根据细胞类型(由细胞注释文件中定义)
                               hclust_method = "ward.D2",# ward.D2 方法进行层次聚类
                               analysis_mode = "subclusters", # 默认是samples，推荐是subclusters
                               denoise = TRUE, # 去噪音
                               HMM = F,  ##特别耗时间,是否要去背景噪音
                               plot_steps = F, #不在每个步骤后生成图形。
                               leiden_resolution = "auto", #可以手动调参数
                               num_threads = 16 # 16线程工作，加快速度，一般电脑不建议哈
                               )

infercnv_obj_medianfiltered = infercnv::apply_median_filtering(infercnv_obj)

infercnv::plot_cnv(infercnv_obj_medianfiltered, 
                   out_dir='./example_output/',
                   output_filename='infercnv.median_filtered', 
                   x.range="auto",
                   x.center=1,
                   title = "infercnv", 
                   color_safe_pal = FALSE)
                   
setwd("..")                 

rm(list = ls())
options(stringsAsFactiors = F)
library(phylogram)
library(gridExtra)
library(grid)
require(dendextend)
require(ggthemes)
library(tidyverse)
library(Seurat)
library(infercnv)
library(miscTools)

sce <- qread("./4-Corrective-data/sce.all.qs")

dir.create("11-infercnv")
setwd("./11-infercnv")

# 提取想要的细胞
sce <- sce[,sce$celltype %in% c("plasma","fibroblasts","epithelial/cancer cells")]

# 提取infercnv_obj
infer_CNV_obj<-readRDS('./infercnv_output/run.final.infercnv_obj')
expr <- infer_CNV_obj@expr.data
expr[1:4,1:4]
data_cnv <- as.data.frame(expr)
dim(expr)
colnames(data_cnv)
rownames(data_cnv)

table(sce$orig.ident)
table(sce$seurat_clusters)
#记得要改一下sce中的列名
colnames(sce) <- gsub("-","_",colnames(sce))
meta <- sce@meta.data

#要根据参数修改哦，这里的对照组有plasma和fibroblast
if(T){
  tmp1 = expr[,infer_CNV_obj@reference_grouped_cell_indices$plasma]
  tmp = tmp1
  tmp2 = expr[,infer_CNV_obj@reference_grouped_cell_indices$fibroblast]
  tmp= cbind(tmp1,tmp2)
  down=mean(rowMeans(tmp)) - 2 * mean( apply(tmp, 1, sd))
  up=mean(rowMeans(tmp)) + 2 * mean( apply(tmp, 1, sd))
  oneCopy=up-down
  oneCopy
  a1= down- 2*oneCopy
  a2= down- 1*oneCopy
  down;up
  a3= up +  1*oneCopy
  a4= up + 2*oneCopy 
  
  cnv_score_table<-infer_CNV_obj@expr.data
  cnv_score_table[1:4,1:4]
  cnv_score_mat <- as.matrix(cnv_score_table)
  
  # Scoring
  cnv_score_table[cnv_score_mat > 0 & cnv_score_mat < a2] <- "A" #complete loss. 2pts
  cnv_score_table[cnv_score_mat >= a2 & cnv_score_mat < down] <- "B" #loss of one copy. 1pts
  cnv_score_table[cnv_score_mat >= down & cnv_score_mat <  up ] <- "C" #Neutral. 0pts
  cnv_score_table[cnv_score_mat >= up  & cnv_score_mat <= a3] <- "D" #addition of one copy. 1pts
  cnv_score_table[cnv_score_mat > a3  & cnv_score_mat <= a4 ] <- "E" #addition of two copies. 2pts
  cnv_score_table[cnv_score_mat > a4] <- "F" #addition of more than two copies. 2pts
  
  # Check
  table(cnv_score_table[,1])
  # Replace with score 
  cnv_score_table_pts <- cnv_score_mat
  rm(cnv_score_mat)
  # 
  cnv_score_table_pts[cnv_score_table == "A"] <- 2
  cnv_score_table_pts[cnv_score_table == "B"] <- 1
  cnv_score_table_pts[cnv_score_table == "C"] <- 0
  cnv_score_table_pts[cnv_score_table == "D"] <- 1
  cnv_score_table_pts[cnv_score_table == "E"] <- 2
  cnv_score_table_pts[cnv_score_table == "F"] <- 2
   
  cnv_score_table_pts[1:4,1:4]
  str(as.data.frame(cnv_score_table_pts[1:4,1:4])) 
  cell_scores_CNV <- as.data.frame(colSums(cnv_score_table_pts))
  
  colnames(cell_scores_CNV) <- "cnv_score" 
}

#可视化
head(cell_scores_CNV) 
score=cell_scores_CNV
head(score)
meta$totalCNV = score[match(colnames(sce),
                            rownames(score)),1] 
p <- ggplot(meta, aes(x= celltype, 
                 y=totalCNV, 
                 fill=celltype)) +
            geom_boxplot();print(p) 
ggsave("totalCNV.pdf",plot = p, width = 9,height = 7,dpi = 300)