infercnv输入文件的制作 / 开普饭

有粉丝反映跟着我们的教程：使用inferCNV分析单细胞转录组中拷贝数变异，但是第一步3个输入文件就制作失败，值得单独写教程强调一下这个解决方案。当然了，如果你还卡在第一步安装R包，请看我昨天在生信菜鸟团的教程：有些R包是你的电脑操作系统缺东西，但也有一些不是。然后就可以查看https://github.com/broadinstitute/inferCNV/wiki 的示例代码：

首先查看3个示例文件

307K Dec 30 15:21 gencode_downsampled.EXAMPLE_ONLY_DONT_REUSE.txt 5.3K Dec 30 15:21 oligodendroglioma_annotations_downsampled.txt 2.8M Dec 30 15:21 oligodendroglioma_expression_downsampled.counts.matrix.gz

制作基因坐标文件文件

示例文件 gencode_downsampled.EXAMPLE_ONLY_DONT_REUSE.txt 里面是：

WASH7P chr1 14363 29806 LINC00115 chr1 761586 762902 NOC2L chr1 879584 894689 MIR200A chr1 1103243 1103332 SDF4 chr1 1152288 1167411 UBE2J2 chr1 1189289 1209265 SCNN1D chr1 1215816 1227409 ACAP3 chr1 1227756 1244989 PUSL1 chr1 1243947 1247057

很简单，就是基因ID，染色体及起始终止坐标，当然，是需要排序的。

制作样本分组矩阵文件

示例文件 oligodendroglioma_annotations_downsampled.txt 里面：

MGH36_P3_H06 Microglia/Macrophage MGH54_P2_F03 Microglia/Macrophage MGH54_P16_F12 Oligodendrocytes (non-malignant) MGH54_P12_C10 Oligodendrocytes (non-malignant) MGH54_P11_C11 Oligodendrocytes (non-malignant) MGH54_P15_D06 Oligodendrocytes (non-malignant) MGH54_P16_A03 Oligodendrocytes (non-malignant) MGH53_P7_B09 Oligodendrocytes (non-malignant)

也是很简单， 第一列是细胞的名字，第二列是其所属的分组。

制作表达矩阵文件

我比较好奇，这个时候，broad出品的软件居然抛弃了他们自己长久以来推广的gct格式表达矩阵。就是GSEA要求的基因表达芯片数据，文本文件格式（* .gct）

这里面的表达矩阵，就是 R里面的 read.table 可以读取的即可。当然了，需要跟上面两个文件保持一致，比如细胞在表达矩阵的列，所以顺序跟样本分组矩阵示例文件 oligodendroglioma_annotations_downsampled.txt 里面的顺序一致，然后基因在表达矩阵的行，所以顺序跟基因坐标文件的示例文件 gencode_downsampled.EXAMPLE_ONLY_DONT_REUSE.txt 里面保持一致。

文件写出

基本上来说，大家自己制作好表达矩阵，分组信息这两个R语言里面的数据框是没有问题的，比较麻烦的可能是基因信息文件。我这里给出一个我自己的解决方案：

library(infercnv) expFile=system.file("extdata", "oligodendroglioma_expression_downsampled.counts.matrix.gz", package = "infercnv") geneFile=system.file("extdata", "oligodendroglioma_annotations_downsampled.txt", package = "infercnv") groupFiles=system.file("extdata", "gencode_downsampled.EXAMPLE_ONLY_DONT_REUSE.txt", package = "infercnv")

library(AnnoProbe) geneInfor=annoGene(rownames(dat),"SYMBOL",'mouse') 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[rownames(dat) %in% geneInfor[,1],] dat=dat[match( geneInfor[,1], rownames(dat) ),] dim(dat) expFile='expFile.txt' write.table(dat,file = expFile,sep = '\t',quote = F) groupFiles='groupFiles.txt' head(groupinfo) 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)

这个代码还是蛮考验大家的技术实力的。

最后运行inferCNV程序

infercnv_obj = CreateInfercnvObject(raw_counts_matrix=expFile, annotations_file=groupFiles, delim="\t", gene_order_file= geneFile, ref_group_names=c("WT")) ## 这个取决于自己的分组信息里面的

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=tempfile(), cluster_by_groups=TRUE, denoise=TRUE, HMM=TRUE)

然后就成功啦!

当然了，值得注意的是我里面调用了一个函数annoGene，来自于AnnoProbe包，获取基因的坐标哦。具体介绍见：基因类型注释根据基因ID就好了

表达芯片探针ID转换大全

在2019年的尾巴，我推出3个R包，

第一个是整合全部的bioconductor里面的芯片探针注释包。
第二个是整合全部GPL的soft文件里面的芯片探针注释包。
第三个是下载全部的GPL的soft文件里面的探针碱基序列比对后注释包。

配合着详细的介绍：

因为这些包暂时托管在GitHub平台，但是非常多的朋友访问GitHub困难，尤其是我打包了好几百个GPL平台的注释信息后，我的GitHub包变得非常臃肿，大家下载安装困难，所以我重新写一个精简包。也在：芯片探针ID的基因注释以前很麻烦和：芯片探针序列的基因注释已经无需你自己亲自做了, 里面详细介绍了。最重要的是idmap函数，安装方法说到过：芯片探针序列的基因注释已经无需你自己亲自做了, 使用起来也非常简单：

library(AnnoProbe) ids=idmap('GPL570',type = 'soft') head(ids)

仅仅是一句话，就拿到了这个平台的探针的注释信息。需要注意的是，这个函数的type参数，其实是有3个选择，这里我演示的是选择soft这个来源的基因注释信息。

并不是所有的平台都是有soft注释，也不是所有的平台都被我的这个工具囊括哦。

infercnv输入文件的制作