查看原文
其他

CNS图表复现14—检查文献的inferCNV流程

生信技能树 单细胞天地 2022-06-06

分享是一种态度

前言

CNS图表复现之旅前面我们已经进行了13讲,你可以点击图表复现话题回顾。如果你感兴趣也想加入交流群,自己去:你要的rmarkdown文献图表复现全套代码来了(单细胞)找到我们的拉群小助手哈。

前面我们的教程讲到了,自己取全部的上皮细胞,以及部分Fibroblasts和Endothelial_cells细胞来一起运行inferCNV流程,但是得到的结果很诡异,明明是作为二倍体正常细胞参考集的Fibroblasts和Endothelial_cells细胞居然也是在某些染色体上面有明显的CNV情况。为了解决这个问题,让我们一起看看文献自己的inferCNV流程是如何使用的,以及对应的数据集。

首先运行作者自己的代码和数据

那,我们就看看作者自己的代码和数据吧,运行他们的inferCNV流程,看看我们的差异究竟是在哪了?

我注意到文章的脚本里面有这样的一句话:

Save all inferCNV files and run inferCNV in previous version of R

看了看作者准备的3个文件,如下:

183K Aug 30 11:33 NI03_CNV_cell_metadata_shuffle_largefile.txt
544M Aug 30 12:02 NI03_CNV_data_out_all_cells_raw_counts_largefile.txt
671K Aug 30 11:31 NI03_CNV_hg19_genes_ordered_correct_noXY.txt

上面的3个文件作者的制作方式,跟我的大同小异,就不过多介绍啦,然后运行作者的inferCNV代码,如下;

library(infercnv)
infercnv_obj = CreateInfercnvObject(raw_counts_matrix = "NI03_CNV_data_out_all_cells_raw_counts_largefile.txt"
                                    annotations_file = "NI03_CNV_cell_metadata_shuffle_largefile.txt"
                                    gene_order_file = "NI03_CNV_hg19_genes_ordered_correct_noXY.txt"
                                    ref_group_names = c("endothelial_normal""fibroblast_normal"), delim = "\t")
# Make sure that chrmosomes are ordered correctly 
slot(infercnv_obj, "gene_order")[,"chr"] <- factor(slot(infercnv_obj, "gene_order")[,"chr"], 
                                                   levels = c("chr1""chr2","chr3","chr4""chr5""chr6","chr7""chr8""chr9","chr10""chr11""chr12","chr13""chr14""chr15","chr16""chr17""chr18","chr19""chr20""chr21","chr22"))
# Run infer CNV
infercnv_all = infercnv::run(infercnv_obj,
                             cutoff=1,  # use 1 for smart-seq, 0.1 for 10x-genomics
                             out_dir= "myresults",  # dir is auto-created for storing outputs
                             cluster_by_groups=F,   # cluster
                             hclust_method="ward.D2", plot_steps=F)

我仔细看了看作者运行inferCNV的代码,差异真的很小,其中cluster_by_groups这个参数仅仅是可视化的选项,不会影响重要的结论。而hclust_method通常呢,影响细胞之间的距离,按照道理并不影响CNV,那么应该是我前面的那些其它参数导致的。

让我们看看这个函数的默认参数:

run(infercnv_obj, cutoff = 1, min_cells_per_gene = 3, out_dir = NULL,
  window_length = 101, smooth_method = c("pyramidinal""runmeans",
  "coordinates"), num_ref_groups = NULL,
  ref_subtract_use_mean_bounds = TRUE, cluster_by_groups = FALSE,
  cluster_references = TRUE, k_obs_groups = 1,
  hclust_method = "ward.D2", max_centered_threshold = 3,
  scale_data = FALSE, HMM = FALSE, HMM_transition_prob = 1e-06,
  HMM_report_by = c("subcluster""consensus""cell"),
  HMM_type = c("i6""i3"), HMM_i3_pval = 0.05, HMM_i3_use_KS = TRUE,
  BayesMaxPNormal = 0.5, sim_method = "meanvar",
  sim_foreground = FALSE, reassignCNVs = TRUE,
  analysis_mode = c("samples""subclusters""cells"),
  tumor_subcluster_partition_method = c("random_trees""qnorm",
  "pheight""qgamma""shc"), tumor_subcluster_pval = 0.1,
  denoise = FALSE, noise_filter = NA, sd_amplifier = 1.5,
  noise_logistic = FALSE, outlier_method_bound = "average_bound",
  outlier_lower_bound = NA, outlier_upper_bound = NA,
  final_scale_limits = NULL, final_center_val = NULL, debug = FALSE,
  num_threads = 4, plot_steps = FALSE, resume_mode = TRUE,
  png_res = 300, plot_probabilities = TRUE, save_rds = TRUE,
  save_final_rds = TRUE, diagnostics = FALSE,
  remove_genes_at_chr_ends = FALSE, prune_outliers = FALSE,
  mask_nonDE_genes = FALSE, mask_nonDE_pval = 0.05,
  test.use = "wilcoxon", require_DE_all_normals = "any",
  hspike_aggregate_normals = FALSE, no_plot = FALSE,
  no_prelim_plot = FALSE, output_format = "png", useRaster = TRUE,
  up_to_step = 100)

多到让人头皮发麻!

其中文献运行infercnv::run的时候,下面两个参数,都是默认值:

HMM参数 when set to True, runs HMM to predict CNV level (default: FALSE)
denoise  If True, turns on denoising according to options below (default: FALSE) 

而我运行的时候,把这两个参数都设置为了T,运行该文献他自己的数据集和文献代码后,运行的日志文件如下所示:

INFO [2020-10-19 11:17:44] ::process_data:Start
INFO [2020-10-19 11:17:44] Creating output path myresults
INFO [2020-10-19 11:17:44] Checking for saved results.
INFO [2020-10-19 11:17:44] 

 STEP 1: incoming data

INFO [2020-10-19 11:18:19] 

 STEP 02: Removing lowly expressed genes

INFO [2020-10-19 11:18:19] ::above_min_mean_expr_cutoff:Start
INFO [2020-10-19 11:18:19] Removing 5929 genes from matrix as below mean expr threshold: 1
INFO [2020-10-19 11:18:20] validating infercnv_obj
INFO [2020-10-19 11:18:20] There are 14467 genes and 7181 cells remaining in the expr matrix.
INFO [2020-10-19 11:18:24] no genes removed due to min cells/gene filter

INFO [2020-10-19 12:19:51] plot_cnv_references:Start
INFO [2020-10-19 12:19:51] Reference data size: Cells= 1000 Genes= 14467
INFO [2020-10-19 12:20:07] plot_cnv_references:Number reference groups= 2
INFO [2020-10-19 12:20:07] plot_cnv_references:Plotting heatmap.
INFO [2020-10-19 12:20:10] Colors for breaks:  #00008B,#24249B,#4848AB,#6D6DBC,#9191CC,#B6B6DD,#DADAEE,#FFFFFF,#EEDADA,#DDB6B6,#CC9191,#BC6D6D,#AB4848,#9B2424,#8B0000
INFO [2020-10-19 12:20:10] Quantiles of plotted data range: 0.544327010335531,0.938892738431902,1,1.06175861606737,1.53099452887365
INFO [2020-10-19 12:20:10] plot_cnv_references:Writing reference data to myresults/infercnv.references.txt

耗时约1个小时(主要的时间花在了第15步),关键问题是,他得到的CNV非常漂亮。也就是说如果不考虑数据集的差异,这个时候得到的结论是HMM参数和 denoise参数都应该是默认值才行啊。

infercnv

然后运行我的代码在作者的数据

跟上一讲我们的代码大同小异,如下:

rm(list = ls())

dat=read.table('NI03_CNV_data_out_all_cells_raw_counts_largefile.txt',
               header = T,sep = '\t')
dim(dat)
library(AnnoProbe)
geneInfor=annoGene(rownames(dat),"SYMBOL",'human')
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)

groupFiles='groupFiles.txt'
groupinfo=read.table('NI03_CNV_cell_metadata_shuffle_largefile.txt',header = F,sep = '\t')
table(groupinfo$V2)
dim(groupinfo)
head(groupinfo)
table(groupinfo$V1 %in% colnames(dat))
write.table(groupinfo,file = groupFiles,sep = '\t',quote = F,col.names = F,row.names = F)
dat=dat[, colnames(dat) %in% groupinfo$V1]

expFile='expFile.txt'
write.table(dat,file = expFile,sep = '\t',quote = F)

head(geneInfor)
geneFile='geneFile.txt'
write.table(geneInfor,file = geneFile,sep = '\t',quote = F,col.names = F,row.names = F)


infercnv_obj = CreateInfercnvObject(raw_counts_matrix=expFile,
                                    annotations_file=groupFiles,
                                    delim="\t",
                                    gene_order_file= geneFile,
                                    ref_group_names = c("endothelial_normal""fibroblast_normal") )  
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='jimmy_results'
                             cluster_by_groups=TRUE
                             denoise=TRUE,
                             HMM=TRUE)

这个时候,我的时间主要是花费在了第STEP 18: Run Bayesian Network Model on HMM predicted CNV's

INFO [2020-10-20 09:25:35] Creating the following Directory:  jimmy_results/BayesNetOutput.HMMi6.hmm_mode-samples
INFO [2020-10-20 09:25:35] Initializing new MCM InferCNV Object.
INFO [2020-10-20 09:25:35] validating infercnv_obj
INFO [2020-10-20 09:25:36] Total CNV's:  1230
INFO [2020-10-20 09:25:36] Loading BUGS Model.
INFO [2020-10-20 09:25:38] Running Sampling Using Parallel with  4 Cores

中间调用了我MAC电脑的4个核心去计算,值得一提的是,因为等待时间过长,经常出现错误!!!,如下所示:

INFO [2020-10-19 15:46:13] Initializing new MCM InferCNV Object.
INFO [2020-10-19 15:46:13] validating infercnv_obj
INFO [2020-10-19 15:46:14] Total CNV's:  1239
INFO [2020-10-19 15:46:14] Loading BUGS Model.
INFO [2020-10-19 15:46:16] Running Sampling Using Parallel with  4 Cores
INFO [2020-10-19 18:21:03] Obtaining probabilities post-sampling
Error in do.call(rbind, mcmc[[j]]) : second argument must be a list
In addition: Warning message:
In parallel::mclapply(seq_along(obj@cell_gene), FUN = par_func,  :
  scheduled cores 1, 2, 3, 4 did not deliver results, all values of the jobs will be affected

运行了6次,都失败,让我很恼火,差不多的数据和代码,为什么我自己运行十多分钟即可,文章的这个需要十几个小时。

最后

后来我仔细比较了,发现自己的数据里面,是因为  366 genes and 7044 cells  , 得到是CNV数量太少了(第18步写的是:Total CNV's:  31 )计算量比较小,所以十几分钟就结束了。

但是文章的这个数据集呢, Total CNV's:  1229 太多了,耗费计算时间和资源有点过分了。这个数据量:14869 genes and 7181 cells 其实不能选择  denoise=TRUE以及HMM=TRUE,都应该是用默认的FALSE即可。

所以我真正需要比较的是,为什么我自己运行inferCNV的时候的输入数据跟作者的差异这么大!!!

咱们明明都是取全部的上皮细胞,以及部分Fibroblasts和Endothelial_cells细胞来一起运行inferCNV流程啊!!!



往期回顾

CNS图表复现13—使用inferCNV来区分肿瘤细胞的恶性与否

细胞身份何以在分裂中得以保持?

年薪40万起诚聘高级生物信息学经理(迪安诊断-国内率先上市的第三方医学诊断机构)






如果你对单细胞转录组研究感兴趣,但又不知道如何入门,也许你可以关注一下下面的课程



看完记得顺手点个“在看”哦!


生物 | 单细胞 | 转录组丨资料每天都精彩

长按扫码可关注



您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存