辛丑年初始，風風專欄打算開一個新的系列——從零開始單細胞系列一。看到後面帶了個一，大家就可以猜到這是一個長期系列，如果後面有新的內容或者你們有新的單細胞相關的需求，我們可以再開個二，與時俱進嘛！今天大年初一嘛，所以主要是跟大家聊一聊這個系列會有什麼內容。
這次的系列推文總共分為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))pvalues_mat <- matrix(NA, nrow=length(pathway_names), ncol=length(cell_types), dimnames = (list(pathway_names, cell_types)))
for(p in pathway_names){ genes <- pathways[[p]] genes_comm <- intersect(genes, rownames(norm_tpm)) if(length(genes_comm) < 5) next pathway_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)) keep <- colnames(mean_exp_eachCellType)[colAlls(mean_exp_eachCellType > 0.001)] if(length(keep) < 3) next pathway_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 * 3 col_lower <- col_q1 / 3 outliers <- apply(ratio_exp_eachCellType, 1, function(x) {any( (x>col_upper)|(x<col_lower) )} ) if(sum(!outliers) < 3) next keep <- 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:5000 weight_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) <- times mean_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) <- times colnames(shuffle_results) <- cell_types for(c in cell_types){ if(is.na(mean_expression_shuffle[p,c])) next if(mean_expression_shuffle[p,c]>1){ pval <- sum(shuffle_results[,c] > mean_expression_shuffle[p,c]) / 5000 }else if(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(i in 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)] <- 1
mybreaks <- 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」獲得本次數據和代碼，
牛年快樂，希望大家牛年文章發得越來越牛，我們下回見吧！
撥開數據迷雲霧，得見缺失值真面目
師兄用臨床數據一年發五篇文章，原來靠這個訣竅！（建議收藏！）
單個組學玩膩了，來試試20分文章喜歡的多組學聯合聚類吧！
頂刊JCO的ROC曲線繪製秘訣，讓你的臨床研究再上一個level！
手把手教你復現一篇CNS級別美圖！附代碼，建議收藏！
CNS級別的美圖復現起來難不難？小白也能快速上手（附代碼）（二）
小白課堂 | 學個R包，復現一張CNS級別美圖（附代碼）
學一個R包給文章加一張CNS級別的美圖（四）
學一個R包給文章加一張CNS級別的美圖（五）
學一個R包給文章加一張CNS級別的美圖（六）
拆和疊？原來你和高階論文的圖片就差這點
可繁可簡，可鹽可甜，這才是upset的真面目吧！
高分文章不講xx，竟然用這種圖片
盪氣迴腸一柄劍，看T細胞捨生取義
Cancer Cell的一張美圖送給你，新年也發個CC
嘀嘀嘀，美圖系列末班車已到站，來收拾最後的圖片下車
歡迎大家關註解螺旋生信頻道-挑圈聯靠公號~
2020
感謝所有小夥伴的一路陪伴
開心這一路和大家共同成長
2021
我們仍要一起並肩前行
朝更新的目標一起努力
為了感謝大家一路的支持
在春節大年初五迎財神時
酸談將進行一場福利抽獎直播
純抽獎part、全新福利周邊
大家一定記得來觀看直播奧
直播信息
直播時間：大年初五
直播地點：B站解螺旋直播間
直播內容：福利直播抽獎party
直播地址：https://live.bilibili.com/8116225
掃碼直達直播間
長按識別二維碼免費包郵領取！