大年初一，单细胞分析从nature communication开始

辛丑年，万事胜意

大家好，我是风风。今天是大年初一，先给大家拜个年，祝大家辛丑年心想事成，万事胜意！年年岁岁花相似，岁岁年年人不同，不知道今年的你，又是在跟谁团圆？是在家中？在宿舍？在单位？还是，在路上独自一人呢？不管相伴还是远离，不管匆忙还是从容，希望在这一年的第一天，你都是快乐的心情，把过去一年的忧愁和烦恼都忘掉，因为，今年花落颜色改，明年花开复谁在？好好珍惜现在，过去不必挂怀。我总是觉得元旦不算“过年”，只有到了春节才是真正意义上的“过年”，大概是因为只有到了春节，小孩才会有压岁钱，门外才会响起鞭炮，当然，最重要的是，和父母的每年一聚，今年，你回家了吗？《目送》一书中写道：“我慢慢地、慢慢地了解到，所谓父女母子一场，只不过意味着，你和他的缘分就是今生今世不断地在目送他的背影渐行渐远。你站立在小路的这一端，看着他逐渐消失在小路转弯的地方，而且，他用背影默默告诉你：不必追。”我想，千言万语都概括在了这段话里。所以，今天欢庆之余，即使是打电话也好，别忘了跟家中长辈说一句“新年快乐，我很好！”

辛丑年初始，风风专栏打算开一个新的系列——从零开始单细胞系列一。看到后面带了个一，大家就可以猜到这是一个长期系列，如果后面有新的内容或者你们有新的单细胞相关的需求，我们可以再开个二，与时俱进嘛！今天大年初一嘛，所以主要是跟大家聊一聊这个系列会有什么内容。

这次的系列推文总共分为3个部分：

第一部分：单细胞文章解读。

计划为大家准备3篇论著和两篇综述，3篇论著已经找好了，分为三个阶段，分别来自：Cell Report（没错，就是前段时间被某个事件误伤的那本Cell 子刊）、Nature Communication以及Nature Medicine。综述的文章我还没找好，所以没办法放出来截图。

我们先进行单细胞论著的文章解读，对单细胞系列的文章有一个大致的认识之后，接着开始通过两篇综述学习单细胞分析的一些常识和前沿知识。这也是酸菜校长在《三十六策》教给大家的一种方法：通过阅读综述学习一个领域的常识。两篇综述解读完如果觉得不够的话，可以再加两篇，这个我们后续再说。

第二部分：模块分析。

按照生信全书中酸菜校长的四字心法：挑、圈、联、靠，我们分模块学习在单细胞分析中对应的内容，内容包括但不局限于GEO单细胞数据下载、数据质控和过滤、批次效应和数据整合、聚类和分类分析、细胞注释、鉴定标志物和拟时序分析等等内容。

第三部分：实战演练与文章复现。

学习完第二部分的模块分析内容之后，我们需要把模块组合起来应用到自己的分析中从而产出自己的文章，最好的办法就是复现别人文章，然后对比和作者结果的差异。目前也是计划复现两篇文章的部分图片，也是分别来自：Cell Report和Nature Communication，先说好，不是完整复现啊哈哈哈，不然我会罢工的！如果大家有觉得比较容易上手又跟第二部分所学的知识比较接近的文章想要复现的话，也欢迎来提供文献，看看我能不能驾驭得住（驾驭不住那我也是爱莫能助/(ㄒoㄒ)/~~）。

好啦，就是以上三部分的内容，考虑到大家可能R也还不太熟悉，甚至可能跟我一样“菜”，更别说其他的编程语言了，所以这个系列全部使用R语言来完成，不涉及到python和Linux（Flag我先立在这，就看会不会倒了），如果R语言是零基础，欢迎参加我们的R语言基础训练营或者参加B站直播教小白学习R语言的直播课，直播课目前用的讲义是风风打赏营第一期的一部分讲义，所以参加过打赏营的各位就不要去浪费时间啦！如果还是对自己没信心的话，在留言区告诉我，我们可以用2-3期推文来一起从0开始R语言（安装软件不算，以前有推文了我记得）。对了对了对了，最后，单细胞数据都很大，对电脑要求也高，希望你们有个心里准备，比如我的电脑，写教程时候竟然崩溃了两次？？？完蛋！

如果你能看到这里，那我相信你是“真爱”了，那总要有点惊喜嘛，我们提前来复现下面这张Nature Communication的图片吧，体验一下单细胞的图片和分析代码吧，年初一，我们整个简单点的，大家先看一看，图个开心，详细的代码讲解我们后面复现的时候会再具体讲§(*￣▽￣*)§

数据读取

我们先把数据读取进来，并且对数据进行相应处理：

rm(list=ls())source("utils.R")library(stringr)library(reshape2)library(scales)library(scater)library(pheatmap)pathway_file<- "KEGG_metabolism.gmt"# 加载数据selected_impute_sce<- readRDS("selected_impute_sce.rds")pathways<- gmtPathways(pathway_file)pathway_names<- names(pathways)all_cell_types<- as.vector(selected_impute_sce$cellType)cell_types<- unique(all_cell_types)

我们知道，每个通路都有很多基因，并且这些通路之间的基因可能有交叉，因此，我们计算一下每个基因参与的通路数目，同时计算通路活性：

# 查看每个基因参与的代谢通路数量gene_pathway_number<- num_of_pathways(pathway_file,rownames(selected_impute_sce)[rowData(selected_impute_sce)$metabolic])set.seed(0413)norm_rds_file<- file.path("Deconvolution_tpm.rds")norm_tpm<- readRDS(norm_rds_file)## 计算通路活性mean_expression_shuffle<- matrix(NA,nrow=length(pathway_names),ncol=length(cell_types),dimnames = list(pathway_names,cell_types))mean_expression_noshuffle<- matrix(NA,nrow=length(pathway_names),ncol=length(cell_types),dimnames = list(pathway_names,cell_types))# 计算p值pvalues_mat<- matrix(NA,nrow=length(pathway_names),ncol=length(cell_types),dimnames = (list(pathway_names, cell_types)))for(pin pathway_names){genes<- pathways[[p]]genes_comm<- intersect(genes, rownames(norm_tpm))if(length(genes_comm)< 5) nextpathway_metabolic_tpm<- norm_tpm[genes_comm, ]pathway_metabolic_tpm<- pathway_metabolic_tpm[rowSums(pathway_metabolic_tpm)>0,]mean_exp_eachCellType<- apply(pathway_metabolic_tpm, 1,function(x)by(x,all_cell_types,mean)) # 去除表达量在所有细胞中都为0的基因keep<- colnames(mean_exp_eachCellType)[colAlls(mean_exp_eachCellType > 0.001)]if(length(keep)< 3) next # 使用最小值代替0pathway_metabolic_tpm<- pathway_metabolic_tpm[keep,]pathway_metabolic_tpm<- t( apply(pathway_metabolic_tpm,1,function(x){x[x<=0] <- min(x[x>0]);x} ))pathway_number_weight = 1 / gene_pathway_number[keep,]mean_exp_eachCellType<- apply(pathway_metabolic_tpm, 1,function(x)by(x,all_cell_types,mean))ratio_exp_eachCellType<- t(mean_exp_eachCellType) / colMeans(mean_exp_eachCellType) # 去除异常比值col_quantile<- apply(ratio_exp_eachCellType,2,function(x)quantile(x,na.rm=T))col_q1<- col_quantile["25%",]col_q3<- col_quantile["75%",]col_upper<- col_q3 * 3col_lower<- col_q1 / 3outliers<- apply(ratio_exp_eachCellType,1,function(x){any( (x>col_upper)|(x<col_lower) )} )if(sum(!outliers)< 3) nextkeep<- names(outliers)[!outliers]pathway_metabolic_tpm<- pathway_metabolic_tpm[keep,]pathway_number_weight = 1 / gene_pathway_number[keep,]mean_exp_eachCellType<- apply(pathway_metabolic_tpm,1,function(x)by(x,all_cell_types,mean))ratio_exp_eachCellType<- t(mean_exp_eachCellType) / colMeans(mean_exp_eachCellType)mean_exp_pathway<- apply(ratio_exp_eachCellType,2,function(x)weighted.mean(x, pathway_number_weight/sum(pathway_number_weight)))mean_expression_shuffle[p,] <- mean_exp_pathway[cell_types]mean_expression_noshuffle[p,] <- mean_exp_pathway[cell_types]group_mean<- function(x){sapply(cell_types,function(y)rowMeans(pathway_metabolic_tpm[,shuffle_cell_types_list[[x]]==y,drop=F]))}column_weigth_mean<- function(x){apply(ratio_exp_eachCellType_list[[x]],2,function(y) weighted.mean(y, weight_values))}times<- 1:5000weight_values<- pathway_number_weight/sum(pathway_number_weight)shuffle_cell_types_list<- lapply(times,function(x) sample(all_cell_types)) names(shuffle_cell_types_list)<- timesmean_exp_eachCellType_list<- lapply(times,function(x) group_mean(x))ratio_exp_eachCellType_list<- lapply(times,function(x)mean_exp_eachCellType_list[[x]] / rowMeans(mean_exp_eachCellType_list[[x]]))mean_exp_pathway_list<- lapply(times,function(x) column_weigth_mean(x))shuffle_results<- matrix(unlist(mean_exp_pathway_list),ncol=length(cell_types),byrow = T) rownames(shuffle_results)<- timescolnames(shuffle_results)<- cell_typesfor(cin cell_types){if(is.na(mean_expression_shuffle[p,c]))nextif(mean_expression_shuffle[p,c]>1){pval<- sum(shuffle_results[,c] > mean_expression_shuffle[p,c]) / 5000 }elseif(mean_expression_shuffle[p,c]<1){pval<- sum(shuffle_results[,c] < mean_expression_shuffle[p,c]) / 5000}if(pval>0.01)mean_expression_shuffle[p, c] <- NA pvalues_mat[p,c]<- pval}}all_NA<- rowAlls(is.na(mean_expression_shuffle))mean_expression_shuffle<- mean_expression_shuffle[!all_NA,]

好了，现在准备数据完成，我们来绘制一下原文中的热图，需要说明的是，由于有随机种子，所以我们复现的结果可能跟原文有出入，总体结果没有影响，我们来看看绘图：

# 热图绘制dat<- mean_expression_shufflesort_row<- c()sort_column<- c()for(iin colnames(dat)){select_row<- which(rowMaxs(dat,na.rm = T) == dat[,i])tmp<- rownames(dat)[select_row][order(dat[select_row,i],decreasing = T)]sort_row<- c(sort_row,tmp)}sort_column<- apply(dat[sort_row,],2,function(x) order(x)[nrow(dat)])sort_column<- names(sort_column)dat[is.na(dat)]<- 1mybreaks<- c(seq(0,0.5, length.out=33),seq(0.51,1.5, length.out=33),seq(1.51,max(dat),length.out=34))color<- colorRampPalette(c("skyblue","white","red")) (100)pheatmap(dat[sort_row,sort_column],cluster_cols = F,cluster_rows = F,color=color,breaks=mybreaks)

好了，我知道很多人都觉得看不懂，不过没关系，我也不懂！！！那就后面一起学习吧，希望经过第二部分的学习，在第三部分文章复现的时候，你能够看懂代码，并且用自己的数据画图吧！Flag立在这里，希望不会倒下！如果还是觉得有困难，并且大家也比较感兴趣的话，后面等空下来再给大家录个精品课，当然能文字解决的的事我们尽量别麻烦视频了是不O(∩_∩)O？

按照惯例，后台回复“feng37”获得本次数据和代码，

牛年快乐，希望大家牛年文章发得越来越牛，我们下回见吧！