美文网首页单细胞测序Single Cell RNA Sequencing Analysis
【​单细胞转录组】inferCNV第二次总结

【​单细胞转录组】inferCNV第二次总结

作者: Geekero | 来源:发表于2020-03-23 20:30 被阅读0次

InferCNV的目的:

InferCNV用于探索肿瘤单细胞RNA-Seq数据,以鉴定大规模染色体拷贝数变异的证据,例如整个染色体或染色体大片段的获得或缺失。

理论上是可以协助判断哪些亚群是肿瘤细胞,但实际上效果可能没想象中那么好

一、R包安装

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("infercnv")

下面可选安装:

install.packages("tibble")
install.packages("devtools")
devtools::install_github("bmbroom/tsvio")
devtools::install_github("bmbroom/NGCHMR", ref="stable")
devtools::install_github("broadinstitute/inferCNV_NGCHM")

在常规shell中键入以下内容来下载NGCHM Java应用程序:

wget http://tcga.ngchm.net/NGCHM/ShaidyMapGen.jar

二、输入文件准备

需要准备3个输入数据:

  • 单细胞RNA-seq表达量的原始矩阵
  • 注释文件,记录肿瘤和正常细胞
  • 基因或染色体位置文件
1. 原始表达矩阵(matrix)

行名是基因,列名是细胞编号。



矩阵可以制表符分隔的文件形式提供。(请注意,还支持稀疏矩阵-请参阅Running-InferCNV)

2. 样本注释文件(分组信息)

样本注释文件用于定义不同的细胞类型,并可选地指示应如何根据样本(即患者)对细胞进行分组。格式仅是两列,以制表符分隔,并且没有列标题。

命名为"cellAnnotations.txt"。一共两列,第一列是对应第一个文件的列名,第二列是细胞的分组

MGH54_P2_C12    Microglia/Macrophage
MGH36_P6_F03    Microglia/Macrophage
MGH54_P16_F12   Oligodendrocytes (non-malignant)
MGH54_P12_C10   Oligodendrocytes (non-malignant)
MGH36_P1_B02    malignant_MGH36
MGH36_P1_H10    malignant_MGH36
3.基因位置信息文件,

命名为"geneOrderingFile.txt"。这个说白了就是基因GTF注释文件,然后继续排序和去重后得到的

一共四列,第一列对应第一个文件的行名,其余三列则是基因的位置。注:基因名不能有重复

基因排序文件提供了每个基因的染色体位置。格式是制表符分隔的,没有列标题,只提供基因名称,染色体和基因跨度:

WASH7P  chr1    14363   29806
LINC00115       chr1    761586  762902
NOC2L   chr1    879584  894689
MIR200A chr1    1103243 1103332
SDF4    chr1    1152288 1167411
UBE2J2  chr1    1189289 1209265

要分析的计数矩阵中的每个基因都应该在这个基因排序文件中提供相应的基因名称和位置信息。

三、运行InferCNV

用10X单细胞转录组数据运行inferCNV实操中踩过的坑:
  1. 必须要输出到一个文件,再构造inferCNV对象,否则会报错
  2. 用于构造inferCNV对象的三个文件,必须查看一下格式有没有问题,必须是否与官网一致
2个月前写的比较偏向官方说明的流程:

#Robin 2020/01/18
#For inferCNV plot:Epithelial cell(reference B cell)
library(infercnv)
library(Seurat)
library(dplyr)

################ File prepare ########
# load data
integrated_1<-readRDS('/share/nas1/Data/Users/luohb/Personalization/result3/seurat_obj/integrated_1_new2.rds')

##file1: cellAnnotations.txt
cellann1<-read.table(file='./integrated_1/03.Cell_Cluster_Recognition/cellAnn.txt', header=F, sep="\t", col.names=c('cell', 'cell_type'))
cellann1$cluster<-integrated_1@meta.data$seurat_clusters
rownames(cellann1)<-cellann1[,1]


clusterAnn1<-read.table("./integrated_1/03.Cell_Cluster_Recognition/clusterAnn.txt", header=F, sep='\t', col.names=c("cluster", "cell_type"))

###mapping
map1=as.character(clusterAnn1[,2])
names(map1)<-as.character(clusterAnn1[,1])
cellann1$cluster_type<-plyr::revalue(as.character(cellann1$cluster), map1)
cellann1$cluster<-paste0('cluster', cellann1$cluster)
# ann1<-ann1[ann1$cluster_type=="Epithelial_cells", c('cell', "cluster")]

###reshape
ann1<-cellann1[cellann1$cluster_type=="Epithelial_cells"|cellann1$cluster_type=="B_cell", c('cell', 'cluster_type', 'cluster')]
ann1$cluster<-ifelse(ann1$cluster_type=="B_cell", "B_cell", ann1$cluster)
ann1<-ann1[, c('cell', 'cluster')]

##file2 geneOrderingFile.txt
gtf<-read.table(file='/share/nas1/Data/Seurate/geneOrderingFile_sort.txt', header=F, sep='\t', col.names=c('gene','chr', 'start', 'end'))
tail(gtf)
gtf1<-gtf[gtf$gene %in% rownames(integrated_1@assays$RNA@counts), ]

##file3 counts.matrix
count1<-integrated_1@assays$RNA@counts
# count1_sort <- count1[rownames(count1) %in% gtf1$gene, colnames(count1) %in% ann1$cell]



############# write table #########
setwd('/share/nas1/Data/Users/luohb/Personalization/result6/')
write.table(count1, file = "count1.matrix", sep = '\t', quote = F, row.names = T, col.names = T)
# write.table(count1_sort, file = "count1.matrix", sep = '\t', quote = F, row.names = T, col.names = T)
write.table(ann1, file = "cellAnnotations1.txt", sep = '\t', quote = F, row.names = F, col.names = F)
write.table(gtf1, file = "geneOrderingFile1.txt", sep = '\t', quote = F, row.names = F, col.names = F)


######## Create inferCNV object ##############
path<-"./result6/"
dir.create(path)
setwd(path)

##Create inferCNV object
infercnv_obj = CreateInfercnvObject(
  raw_counts_matrix='./count1.matrix',
  annotations_file='./cellAnnotations1.txt',
  delim="\t",
  gene_order_file="./geneOrderingFile1.txt",
  ref_group_names=c("B_cell"))

# perform infercnv operations to reveal cnv signal
infercnv_obj = infercnv::run(infercnv_obj,
                             cutoff=0.1,  # use 1 for smart-seq, 0.1 for 10x-genomics
                             out_dir="./integrated_1",  # dir is auto-created for storing outputs
                             cluster_by_groups=F,   # 选择TRUE是按样本分组 改为FALSE会进行按另一个参数k_obs_groups给出的分组数(默认为1)进行分组
                             plot_steps=FALSE,
                             denoise=T,
                             HMM=T  # 是否基于HMM预测CNV,True的话时间很久
)
今天用函数做了几次小封装后的流程:

在后台运行脚本

#Robin 2020/03/23
#For inferCNV plot:all cluster cell(reference B cell)
library(infercnv)
library(Seurat)
library(dplyr)

######## Run inferCNV #############
# Create inferCNV object
two_steps_inferCNV<-function(raw_counts_matrix, annotations_file, gene_order_file,
                             ref_group_names, cutoff=0.1, out_dir, cluster_by_groups=F,
                             plot_steps=F, denoise=T, HMM=F, delim="\t"){
  
  
                             infercnv_obj = CreateInfercnvObject(
                                         raw_counts_matrix=raw_counts_matrix,
                                         annotations_file=annotations_file,
                                         delim=delim,
                                         gene_order_file=gene_order_file,
                                         ref_group_names=ref_group_names)  #change background cell type
  
                                   ## perform infercnv operations to reveal cnv signal
                            infercnv_obj = infercnv::run(infercnv_obj,
                                          cutoff=cutoff,  # use 1 for smart-seq, 0.1 for 10x-genomics
                                          out_dir=out_dir,  # dir is auto-created for storing outputs
                                          cluster_by_groups=cluster_by_groups,   # 选择TRUE是按样本分组 改为FALSE会进行按另一个参数k_obs_groups给出的分组数(默认为1)进行分组
                                          plot_steps=plot_steps,  
                                          denoise=denoise,
                                          HMM=HMM  # 是否基于HMM预测CNV,True的话时间很久
  )
}

##cluster16
two_steps_inferCNV(raw_counts_matrix = "./cluster16_count.matrix", 
                   annotations_file = "./cluster16_cellAnnotations.txt", 
                   gene_order_file = "./geneOrderingFile.txt",
                   ref_group_names = "B-cells",
                   out_dir = "./cluster16"
)
参数说明

Infercnv::run的参数非常之多,总体上分为如下几类:

  • 基本设置
  • 平滑参数
  • 基于肿瘤亚群的下游分析(HMM或non-DE-masking)
  • 根据 BayesNet P(Normal) 过滤低可信度HMM预测结果
  • 肿瘤亚群细分
  • 去噪参数
  • 离群值修剪
  • 其他选项
  • 实验性参数(不稳定)
  • 差异表达分析的实验性参数

按照具体的需求修改不同步骤的参数:

  • 例如聚类默认cluster_by_groups=FALSE会根据k_obs_groups聚类成指定的组数,而层次聚类方法用于计算组间相似度的参数则是hclust_method.
  • 设置HMM=TRUE的计算时间会长于HMM=FALSE,因此可以先设置HMM=FALSE快速查看结果。(其实加不加HMM个人感觉区别不是很大,时间上却长了很多)

四、输出结果

最终会输出很多文件在out_dir的目录下,而实际有用的是下面几个:

  • infercnv.preliminary.png : 初步的inferCNV展示结果(未经去噪或HMM预测)

  • infercnv.png : 最终inferCNV产生的去噪后的热图.

  • infercnv.references.txt : 正常细胞矩阵.

  • infercnv.observations.txt : 肿瘤细胞矩阵.

  • infercnv.observation_groupings.txt : 肿瘤细胞聚类后的分组关系.

  • infercnv.observations_dendrogram.txt : 与热图匹配的肿瘤细胞的newick格式的树状图

五、结果图的解释

得到的quickstart.pdf文件:
  • 正常细胞的表达值绘制在顶部热图中,肿瘤细胞绘制在底部热图中

  • 基因在整个染色体上从左到右排列。

  • 有效地从肿瘤细胞表达数据中减去正常细胞表达数据以产生差异

  • 其中染色体区域扩增显示为红色块,而染色体区域缺失显示为蓝色块

相关文章

网友评论

    本文标题:【​单细胞转录组】inferCNV第二次总结

    本文链接:https://www.haomeiwen.com/subject/xeyiyhtx.html