单细胞转录组之拷贝数变异分析 |
您所在的位置：网站首页 › 基因snv和cnv › 单细胞转录组之拷贝数变异分析 |
|
1.什么是拷贝数变异 拷贝数变异(Copy number variation, CNV)：基因组发生重排而导致的,一般指长度1 kb 以上的基因组片段的拷贝数增加或者减少, 主要表现为亚显微水平的重复或者缺失。因此称为“微”缺失或重复变异。 人类基因组含有约31.6亿个DNA碱基对，组成了23对染色体，每对染色体中一条遗传自父亲，一条遗传自母亲，2条染色体亦称为同源染色体，大小和基因组成基本一致，22对为常染色体，还有一对性染色体，女性为XX，男性为XY。正常⼈类基因组成分通常是以2个拷贝存在，分别来⾃⽗母。常染色体和女性X染色体正常拷贝数值为2，男性X和Y染色体正常拷贝数值为1，当拷贝数检测数值大于正常值时即为重复，小于正常值时为缺失。 CNV在基因组中的存在形式主要有以下⼏种： 2条同源染⾊体拷贝数同时出现缺失；1条同源染⾊体发⽣缺失，1条正常；1条同源染⾊体出现拷贝数重复，另1条正常；1条同源染⾊体出现缺失，另1条出现拷贝数重复；2条同源染⾊体同时出现拷贝数重复。异常的DNA拷贝数变异（CNV）是许多⼈类疾病（如癌症、遗传性疾病、⼼⾎管疾病）的⼀种重要分⼦机制。作为疾病的⼀项⽣物标志，染⾊体⽔平的缺失、扩增等变化已成为许多疾病研究的热点，然⽽传统的⽅法（⽐如G显带，FISH，CGH等）存在操作繁琐，分辨率低等问题，难以提供变异区段的具体信息，单细胞测序为我们提供了一种新的工具和视野去分析CNV。 2.使用R进行CNV分析2.1 数据的准备#加载需要的包和数据 library(Seurat) # devtools::install_github('satijalab/seurat-data') library(SeuratData) library(ggplot2) library(patchwork) library(dplyr) #以之前pbmc的seurat标准流程为基础，进行分析 DimPlot(pbmc) sce=pbmc table( Idents(sce )) table([email protected]$seurat_clusters) table([email protected]$orig.ident) #由于CNV分析要用到infercnv包，其CreateInfercnvObject函数需要三个文件创建对象：The raw_counts_matrix；the gene_order_file, contains chromosome, start, and stop position for each gene；The annotations_file, containing the cell name and the cell type classification #获取counts矩阵 dat=GetAssayData(sce,slot='counts',assay='RNA') dat[1:4,1:4] #获取细胞亚群信息 groupinfo=data.frame(v1=colnames(dat),v2= Idents(sce ) ) head(groupinfo) #这里使用了AnnoProbe包 https://github.com/jmzeng1314/AnnoProbe/ #AnnoProbe是生信技能树健明老师开发的包，用于表达芯片数据分析，但也可以下载GEO数据，进行进行基因注释等功能，可以注释基因并标记其在染色体上的位置 library(AnnoProbe) geneInfor=annoGene(rownames(dat), "SYMBOL",'human') colnames(geneInfor) head(geneInfor) #安装染色体和基因位置排序，并将SYMBOL，chr，start，end这几列留下 geneInfor=geneInfor[with(geneInfor,order(chr, start)),c(1,4:6)] #去重（这里也可以操作去除性染色体，也可以把染色体排序方式改变） geneInfor=geneInfor[!duplicated(geneInfor[,1]),] length(unique(geneInfor[,1])) head(geneInfor) #整理数据 dat=dat[rownames(dat) %in% geneInfor[,1],] dat=dat[match( geneInfor[,1], rownames(dat) ),] dim(dat) head(groupinfo) dat[1:4,1:4] table(groupinfo$v2) head(groupinfo) # 为了节约计算机资源，直接抽样 kp=sample(1:nrow(groupinfo),500) groupinfo=groupinfo[kp,] dat=dat[,kp] # 如果是真实项目，且计算资源足够,忽略这个抽样的操作 #保存文件 expFile='expFile.txt' write.table(dat,file =expFile,sep = '\t',quote = F) groupFiles='groupFiles.txt' 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)2.2 创建CNV对象options(stringsAsFactors = F) expFile='expFile.txt' groupFiles='groupFiles.txt' geneFile='geneFile.txt' infercnv_obj = CreateInfercnvObject(raw_counts_matrix=expFile, annotations_file=groupFiles,delim="\t",gene_order_file= geneFile,ref_group_names=c("NK",'DC','Platelet')) #这里出现了报错，经查看，是expFile和groupFiles保存为txt时，列名出现了不一致，一个保存成了TCGAGCCTTGTGAC-1，一个保存成了TCGAGCCTTGTGAC.1,即expFile所有的-都变成了.，经过查找，并不清楚write.table函数那个参数导致的 #所以，干脆将groupFiles中的-先变成.，再保存 library(do) head(groupinfo$v1) head(Replace(data=groupinfo$v1,from='-1',to='.1')) groupinfo$v1 |
今日新闻 |
推荐新闻 |
| CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |