在inferCNV完成之后，会得到很多的分析结果，其中最重要的肯定是inferCNV.png，但是这张热图中仅有概括性的基因表达情况和分布情况，如何利用inferCNV的分析结果，找到这部分缺失的基因呢？
1. 准备工作
-
infercnv.observations.txt，相关内容描述可在inferCNV——单细胞RNA中的拷贝数变异中查看 -
gencode_v19_gene_pos.txt，基因排序文件，基因排序文件提供每个基因的染色体位置。该格式以制表符分隔且没有列标题，仅提供基因名称、染色体和基因跨度
> head(genecodev1)
V1 V2 V3 V4
1 DDX11L1 chr1 11869 14412
2 WASH7P chr1 14363 29806
3 MIR1302-11 chr1 29554 31109
4 FAM138A chr1 34554 36081
5 OR4G4P chr1 52473 54936
6 OR4G11P chr1 62948 63887
-
inferCNV.png，把这张图片打开，然后数数，可以看到，是第1列和第19列缺失，记住数字即可
image
- 需要的包
library(dplyr)
library(ggplot2)
library(reshape2)
2. 开始操作
library(dplyr)
library(ggplot2)
library(reshape2)
# 定义缺失列或感兴趣的列
chr_num<-c("chr1","chr19") #
# 读取准备文件
exp <- read.table("./file/infercnv.observations.txt",sep = " ",header = T)
genecodev1 <- read.table("./file/gencode_v19_gene_pos.txt",sep = "\t",header = F)
# 定位缺失基因
chr_genes<-genecodev1[genecodev1$V2%in%chr_num,1] #
exp <- exp[chr_genes,]%>%na.omit()
selected_exp <- exp[rowMeans(exp)<1,]%>%as.matrix() # 样本的Expression均值筛选，可以修改为0.9或1左右，根据热图结果自己控制，数字越小越严格
heat_data<-melt(selected_exp) #长宽矩阵转化，为作图准备
# 查看缺失列的expression热图
ggplot(heat_data, aes(Var1, Var2, fill= value)) +
geom_tile()+
xlab("")+
ylab("")+
theme(axis.text.x = element_blank(),
axis.text.y = element_blank())+
scale_fill_gradient2(low ="#00008b",mid = "#ffffff",high = "#8b0000",midpoint = 1)
# 输出genes
write.table(x = row.names(selected_exp),file = "interest_genes.txt",sep = "\t",row.names = F,col.names = F,quote = F)
3. 最终效果
阈值为1时
可以看到基本expression都是小于1的，也有接近与1的情况
image
阈值为0.9时
得到的基因基本都满足条件了
image
当然还有输出基因的结果
这个文件中的基因即可进行后续分析
> interest_genes.txt
[1] "WASH7P" "LINC00115" "NOC2L" "SDF4" "UBE2J2"
[6] "SLC35E2B" "SLC35E2" "GNB1" "PRKCZ" "C1orf86"
[11] "RER1" "PEX10" "PANK4" "WRAP73" "LRRC47"
[16] "DFFB" "C1orf174" "RPL22" "ICMT" "ACOT7"
......
感谢观看，如果有用还请点赞，关注！
网友评论