写在前面，其实现在网上很多对于一个包(一个/种分析）更多是讲分析代码或流程，但是对于分析结果可视化的教程讲的比较少，这导致我这种画图小白无法对结果进行plot，或者可视化不如意。所以我会在这次分析的后面加上一些可视化的教程，希望大家喜欢。
首先还是讲一下inferCNV这个包能干嘛，个人的理解就是它能根据参比，推断(infer)出某一个细胞的拷贝数变异(CNV)的程度，而拷贝数变异与肿瘤细胞的发生发展有着密切的关系，一般来说，肿瘤细胞的拷贝数会比正常的细胞的高，所以inferCNV主要应用于肿瘤的研究中。

inferCNV分析

分析流程
1.第一步：加载包+准备数据集 (这里的scRNA.data加载后的scRNA是我自己的seurat对象)

#加载包，清空
library(Seurat)
library(tidyverse)
library(ggplot2)
library(infercnv)
library(ComplexHeatmap)
library(ggpubr)
rm(list=ls())
#加载数据
load("scRNA.Rdata")
# 随机提取2000个细胞演示    这里可以不比学，我是想跑一下流程而已
tmp <- sample(colnames(scRNA),2000) 
scRNA_harmony <- scRNA[,tmp]

2.第二步：整理分析所需要的数据

#########################################################################################################################################################
#思考一下 如何推断cnv的  gene的表达量（count表达矩阵），gene的位置（基因染色体信息），比较信息（分组信息）
#inferCNV需要三个文件1.count表达矩阵，2.分组信息，3.基因染色体信息

#制作基因染色体位置信息 和提取表达矩阵
dat <- GetAssayData(scRNA_harmony,assay = "RNA",slot = "counts")
dat <- as.data.frame(dat)
expFile=system.file("extdata", "oligodendroglioma_expression_downsampled.counts.matrix.gz", package = "infercnv")  #这个是inferCNV自带参比数据集  用来做infer
geneFile=system.file("extdata", "oligodendroglioma_annotations_downsampled.txt", package = "infercnv")   #这个是inferCNV自带参比数据集  用来做infer
groupFiles=system.file("extdata", "gencode_downsampled.EXAMPLE_ONLY_DONT_REUSE.txt", package = "infercnv")    #这个是inferCNV自带参比数据集  用来做infer

library(AnnoProbe)    #用jimmy老师的包  这个包很强大，如同数据库，这一次我们用这个包做一个基因与其染色体位置
geneInfor=annoGene(rownames(dat),"SYMBOL",'human')    #可能有一些gene 无法anno
colnames(geneInfor)
geneInfor=geneInfor[with(geneInfor, order(chr, start)),c(1,4:6)]      
geneInfor=geneInfor[!duplicated(geneInfor[,1]),]
length(unique(geneInfor[,1]))
head(geneInfor)

dat=dat[match(geneInfor[,1], rownames(dat)),]    #将表达矩阵的基因排列顺序与geneInfor的基因排列顺序弄成一致
rownames(geneInfor) <- geneInfor$SYMBOL   
geneInfor <- geneInfor[,-1]     #这样我们就制作好了染色体位置信息和排列顺序好的count表达矩阵

#制作mate信息
meta <- subset(scRNA_harmony@meta.data,select = c("celltype"))   #假如你后面是想分析每一个群的CNV的话  这里要改成seruat_cluster

#验证1   表达矩阵的列名要与meta的行名一致
identical(colnames(dat),rownames(meta))  
#验证2   表达矩阵的行名要与geneInfor的行名一致
identical(rownames(dat),rownames(geneInfor))

#因此三个输入数据准备好了   dat-表达矩阵  meta-分组信息  geneInfor-基因染色体信息

我大概有这么多基因没有注释到染色体的位置信息

53b6f635ac1c5bdbeb05d1ff02541b1.png
基因的位置信息
83359298cf84d6ebdca3da4082d9bc5.png
meta数据
ea492ca8913e001154d02b05104aedd.png
表达矩阵
1799b93e8dba35564afc811225c442a.png
3.第三步 inferCNV分析

#二步法构建对象
infercnv_obj = CreateInfercnvObject(raw_counts_matrix=dat,
                                    annotations_file=meta,
                                    delim="\t",
                                    gene_order_file=geneInfor,
                                    ref_group_names=c("T cells"))   #选基础的细胞  或者样本 看meta的输入类型   也可以不选择 算法根据平均值来自己算

infercnv_obj = infercnv::run(infercnv_obj,
                             cutoff=1, # cutoff=1 works well for Smart-seq2, and cutoff=0.1 works well for 10x Genomics
                             out_dir="cnv1/", 
                             cluster_by_groups=TRUE,  # 选择TRUE是按样本分组 改为FALSE会进行按另一个参数k_obs_groups给出的分组数（默认为1）进行分组
                             denoise=F,     #是否去噪
                             HMM=T)   # 是否基于HMM预测CNV,True的话时间很久
#结果输出在当前工作路径下的out_dir文件夹下（最终会输出很多文件在out_dir的目录下）   可以直接用里面的热图

#也就可以对热图进行改造  换颜色（用inferCNV的官方的画图函数）
infercnv::plot_cnv(infercnv_obj, #上两步得到的infercnv对象
                   plot_chr_scale = T, #画染色体全长，默认只画出（分析用到的）基因
                   output_filename = "better_plot",output_format = "pdf", #保存为pdf文件
                   custom_color_pal =  color.palette(c("#8DD3C7","white","#BC80BD"), c(2, 2))) #改颜色

按照很多网上的教程，inferCNV已经结束了，但是身为小白的我很不情愿呀，难道我就只能够用他给出来的结果了吗，不！我要自己主导自己想要展示的内容

结果可视化第一种

提取inferCNV的结果 用相线图或者小提琴图画不同组或者细胞中的cnv的结果，就可以比较不同的细胞群或者不同的样本的CNV的差异，对！最重要的就是差异

b5a9fa6b4cf6adc9c8e62d024d6e4ac.png
这样就可以看出哪种细胞的CNV变异最大了，就最有可能是肿瘤细胞
那么这种图应该怎么做呢？图不就是数据的体现吗，我们可以从图中看出X轴为细胞类型，Y轴是该细胞类型的CNV值，那么我们只需要这两个数据就可以了。要做成的表格如下，每一个细胞的CNV数(obs),这个细胞属于什么cluster，什么类型的细胞，来源哪里(orig.ident)
e900ae60061929d6d791b509409b19c.png

#第一步：提取inferCNV的结果   
obs <- read.table("cnv1/infercnv.observations.txt",header = T,check.names = F)     #首先这个obs是每一个基因在每一个细胞的拷贝信息,相当于该基因的拷贝量
if(T){              #可以通过定义obs的元素数值来让差异变大   在后面画图就能够更大的差异   也可以不运行，更真实体现拷贝数（但是可能没有啥差异）
  max(obs)          #根据最大最小值来定义
  min(obs)     
  obs[obs > 0.6 & obs < 0.7] <- -2   #把0.6-0.7定义为数值2  后面依此类推
  obs[obs >= 0.7 & obs < 0.8] <- -1
  obs[obs >= 0.8 & obs < 1.0] <- 0
  obs[obs >= 1.0 & obs <= 1.1] <- 1
  obs[obs > 1.1 & obs <= 1.3] <- 2
  obs[obs > 1.3] <- 2
}

score=as.data.frame(colSums(obs))   #把obs的每一个基因拷贝量加起来，就是这个细胞的总拷贝数obs

#提取meta信息  celltype,seurat_clusters,orig.ident
meta <- subset(scRNA_harmony@meta.data,select = c("celltype","seurat_clusters","orig.ident"))  #提取该细胞的其他的meta信息

#将meta信息添加给score
meta <- rownames_to_column(meta)
score <- rownames_to_column(score)
score <- merge(score,meta,by.x = "rowname",by.y = "rowname")    #这里会可能损失一些细胞   为什么呢，因为在前面infer时，有一些细胞无法推断，会被删掉，但是总体上问题不大

#我们现在可以用ggplot画图了     但是直接这样出图很丑   因为他会根据你的Y轴的大小  所以建议定义y轴范围
ggboxplot(score,"orig.ident","colSums(obs)",fill="orig.ident")   #也可以将orig.ident换成celltype
#+scale_y_continuous(limits = c(0, 1000))

这样就可以知道哪一个样本来源的细胞拷贝数变化大了（我这里只是一个示例数据集，不符合趋势勿喷）

结果可视化第二种

有一天我看了一篇文献，觉得做的还蛮有意思的，也能讲一些东西，我将图片show出来，然后我们就朝着那个目标前进吧

10726e95d30bcb42c8e6743b912eaf4.png
3cfb5a6a30ec811e506d50f3aef5c57.png
文献出处《Single-cell RNA sequencing highlights the role of inflammatory cancer-associated fibroblasts in bladder urothelial carcinoma》

就是相当于给每一个细胞定义了一个CNV的程度，并plot
所以大家无奖竞猜一下，我们画图需要用到什么样的数据呢，有两个哦，一个是细胞每一个基因的表达量，一个是细胞的注释(相当于热图上面的anno)

#第一步
#我们先解决细胞注释的表格
#由于我们要定义CNV的程度，所以首先要将数值变成注释
#加change列,标记inferCNV的变化
max(score$`colSums(obs)`)          #根据最大最小值来定义
min(score$`colSums(obs)`)          #一般变化倍数大于2倍的可以定义为拷贝数变异
k1 = score$`colSums(obs)` < 250*2;table(k1)      #这里的250是我min(score$`colSums(obs)`)的结果 
k2 = score$`colSums(obs)` >= 250*2.5;table(k2)
score <- mutate(score,change = ifelse(k1,"normal",ifelse(k2,"High","Low")))  #要看懂这段代码  ifelse的条件句  假如符合K1标准为normal，不符合就按照ifelse(k2,"High","Low")这个条件再分，假如符合K2标准就High，不符合就Low
table(score$change)   #这样就做好了一个anno
rownames(score) <- score$rowname
score <- score[order(score$`colSums(obs)`),]   #对数据框进行升序排列  如果加decreasing = T就实现降序排列,这样就会对后面热图的细胞顺序进行排序，为什么想对细胞排序呢，因为我想让CNV属于normal，High，Low的细胞排一起，
score <- score[,-c(1,2)]  #删除多余的

#第二步，解决细胞每一个基因的表达量
#由于前面score跟meta merge了，所以这里要对dat进行匹配
dat=dat[,match(rownames(score), colnames(dat))]
identical(rownames(score), colnames(dat))#这样每一个细胞的anno和它的基因表达量就对应了，记住一定要对应，要想明白哈，有些东西虽然不报错，但是你对应错了，结果也是错的。
dat <- log2(dat+1)#这一步的目的是给表达矩阵取log2，增加基因之间的差异

#第三步 画图啦
#想自定义annotation的颜色       
ann_colors = list(
  change = c(normal="black", Low="white",High="red"))   #这里也可以改变多个anno的颜色  

#画热图
library(pheatmap)
pheatmap(dat,
         show_colnames =F,
         show_rownames = F,
         scale = "row",
         cluster_cols = F,                                        #这些参数建议去Google
         cluster_rows = F, 
         annotation_col=score,
         annotation_colors = ann_colors,
         colorRampPalette(c("#00FF00", "white", "#DC143C"))(75),   #改热图的颜色
         breaks = seq(-3,3,length.out = 100))