inferCNV | 如何找到染色体缺失的那部分基因

在inferCNV完成之后，会得到很多的分析结果，其中最重要的肯定是inferCNV.png，但是这张热图中仅有概括性的基因表达情况和分布情况，如何利用inferCNV的分析结果，找到这部分缺失的基因呢？

1. 准备工作

1. infercnv.observations.txt，相关内容描述可在inferCNV——单细胞RNA中的拷贝数变异中查看
2. 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

3. inferCNV.png，把这张图片打开，然后数数，可以看到，是第1列和第19列缺失，记住数字即可
   image
4. 需要的包

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"
......

感谢观看，如果有用还请点赞，关注！