最近聽說定量免疫浸潤很火，于是我也報名參加了果子老師的課程傲茄，跑了幾個R包毅访，你比如說xCell啊，GSVA啊盘榨，MCPcounter啊喻粹，ConsensusTME，ImmuneSubtypeClassifier草巡，當(dāng)然還有大名鼎鼎的CIBERSORT守呜。安裝R包跑實例文件這件事是我經(jīng)常做的，但是很少做的這么系統(tǒng)山憨。相比于跑R包查乒，對我的更大難點在于理解定量免疫浸潤這件事。
也有國產(chǎn)在線的：http://bioinfo.life.hust.edu.cn/ImmuCellAI#!/
定量免疫浸潤在說一件什么事呢郁竟？
從平均人說起
平均人（average man） 以人體測量結(jié)果平均值為基礎(chǔ)建立的人體模型玛迄。其身體結(jié)構(gòu)的各個尺寸均與特定總體的平均值相對應(yīng)，可反映特定群體身體尺寸的平均狀況枪孩。
我們知道憔晒，人有性別/年齡/身高/胖瘦/嗓門大小。蔑舞。拒担。。等等屬性攻询，每個人不同从撼，一般來看東北人和廣東人還是有區(qū)別的，雖然我們不需要一個一個地比較钧栖。
平均人就好比傳統(tǒng)的bulk RNA測序低零，拿到一塊組織去測序，得到一個樣本各個基因的表達情況拯杠。所以我們拿到bulk RNA的結(jié)果掏婶，其準確度相當(dāng)于我們對歪果仁的模糊印象。
定量免疫浸潤說的是我們憑著這個模糊印象潭陪，可以大致推斷出歪果仁某幾種特殊類型的比例（或在某方面的分值）雄妥。這個比例或者分值最蕾，在一定程度上也可以反應(yīng)該地方人口的異質(zhì)性：你看，本來只是一個平均值老厌，但是里面卻蘊含著可愛的異質(zhì)性瘟则。
定量免疫浸潤（Immune infiltration）
首先，“浸潤”是啥意思枝秤？
嚴陣 《牡丹園記》：“這蒙蒙的綠意醋拧，這團團的紅霧，真像剛滴到宣紙上的水彩一樣淀弹，慢慢地
浸潤開來丹壕。”
這里面的“浸潤”很好的詮釋了其“逐漸滲透垦页，引申為積久而發(fā)生作用”的意義雀费。所以“浸潤”有逐漸滲透的意義。
我們提取臨床癌癥組織去測序(bulk的)痊焊，并不是純腫瘤的結(jié)果（腫瘤有純的嗎？來一打）忿峻，在測序結(jié)果中我們會發(fā)現(xiàn)薄啥，一些屬于血管的基因，屬于免疫系統(tǒng)的基因也有表達逛尚。這時候我們就知道垄惧，測的是個混合體，有腫瘤細胞绰寞，有免疫細胞到逊，有血管，有細胞外基質(zhì)滤钱。
定量免疫浸潤通過簡單的基因表達矩陣數(shù)據(jù)觉壶，將免疫細胞類型含量計算出來。隨著算法的發(fā)展件缸，到后來也不僅是免疫細胞了铜靶，可以是很多細胞類型。
基本的套路是計算每個樣本在不同細胞類型中的分值争剿，主要的算法分為兩類，包括:
- 基于Marker gene的算法, 如GSE(V)A
- 基于 deconvolution的算法痊末， 如CIBERSORT
這倆算法前者大家盡管也不懂蚕苇，但是沒有那么嚇人，第二個反卷積(deconvolution)真的是聽著就頭大凿叠。那么涩笤，到底什么是反卷積嚼吞？
在理解這個反卷積之前我們來看看在數(shù)學(xué)上卷積是什么：
卷積就像加減乘除一樣是一種運算，只是低等數(shù)學(xué)用不到辆它，就沒教誊薄，它是指：把二元函數(shù) U(x,y) = f(x)g(y) 卷成一元函數(shù) V(t) ，俗稱降維打擊锰茉。
在我們的bulk RNA的例子中呢蔫，就是把平均值看作單個細胞表達值與某個函數(shù)（和檢測方法有關(guān)）的卷積。獲得樣本bulk RNA 表達量的過程飒筑，就是一個卷積的過程片吊，而反卷積，是卷積運算的逆運算协屡。
鑒于我們bulk RNA 的數(shù)據(jù)已經(jīng)做的很全了俏脊，TCGA數(shù)據(jù)庫里面有大量的已經(jīng)卷積好的表達量，我們可以用反卷積的方法肤晓，看看樣本中到底免疫細胞有著怎樣的比重（分布）爷贫。而這，也這不也是在說明樣本的異質(zhì)性嗎补憾？
比如一個定量免疫浸潤的套路是這樣的：
- 公共數(shù)據(jù)庫下載數(shù)據(jù)
- 基于表達譜做定量免疫浸潤（如CIBERSORT）
- 得到樣本細胞的異質(zhì)性
- 統(tǒng)計（富集分析）
- 驗證
在眾多工具中CIBERSORT的應(yīng)用是最為廣泛的漫萄，所以這篇文獻是有必要打印出來慢慢評鑒的。在這篇文獻中盈匾，我們也看到了一種把bulk RNA和scRNA結(jié)合起來的橋梁：
CIBERSORT函數(shù)在果子老師的課上講的很詳細腾务，這里僅做簡單的復(fù)現(xiàn)，看看CIBERSORT究竟做了什么削饵。更多介紹參見官網(wǎng)：https://cibersort.stanford.edu/
source("CIBERSORT/CIBERSORT.R")
.libPaths("C:\\Users\\86158\\Documents\\R\\win-library\\3.6")
讀入自帶的 signature matrix（細胞類型打分矩陣）岩瘦，并看看他的格式如何：
sig_matrix <- "CIBERSORT/LM22.txt"
sig_matrix
sdf<- readr::read_tsv(sig_matrix)
sdf[1:4,1:4]
# A tibble: 4 x 4
`Gene symbol` `B cells naive` `B cells memory` `Plasma cells`
<chr> <dbl> <dbl> <dbl>
1 ABCB4 556. 10.7 7.23
2 ABCB9 15.6 22.1 653.
3 ACAP1 215. 322. 38.6
4 ACHE 15.1 16.6 22.1
colnames(sdf)
[1] "Gene symbol" "B cells naive" "B cells memory" "Plasma cells"
[5] "T cells CD8" "T cells CD4 naive" "T cells CD4 memory resting" "T cells CD4 memory activated"
[9] "T cells follicular helper" "T cells regulatory (Tregs)" "T cells gamma delta" "NK cells resting"
[13] "NK cells activated" "Monocytes" "Macrophages M0" "Macrophages M1"
[17] "Macrophages M2" "Dendritic cells resting" "Dendritic cells activated" "Mast cells resting"
[21] "Mast cells activated" "Eosinophils" "Neutrophils"
也就是有22種細胞類型。其實我們是可以根據(jù)自己的樣本類型來利用scrna的數(shù)據(jù)來制作這個矩陣的窿撬。制作這個有什么用启昧？下面再說。
讀入表達譜數(shù)據(jù)尤仍。在這里了停一下箫津，三秒禪，應(yīng)該是什么格式的數(shù)據(jù)宰啦，count苏遥？TPM？
CIBERSORT接受的數(shù)據(jù)是不需要log的赡模，如果你一不小心取了log田炭，也不要害怕，他會幫你再去掉漓柑。
mixture_file = 'exprMat.txt'
exp<- read.table(mixture_file,header = T)
exp[1:4,1:4]
rownames.exprMat. LAU125 LAU355 LAU1255
1 A1BG 0.82 0.58 0.81
2 A1CF 0.00 0.01 0.00
3 A2M 247.15 24.88 2307.94
4 A2M-AS1 1.38 0.20 2.60
我們來執(zhí)行核心函數(shù)CIBERSORT：
res_cibersort <- CIBERSORT(sig_matrix, mixture_file, perm=10, QN=TRUE)
B cells naive B cells memory Plasma cells T cells CD8 T cells CD4 naive T cells CD4 memory resting T cells CD4 memory activated
LAU125 0.01146448 0.33118622 0.000000000 0.03587147 0 0.1465500 0.00000000
LAU355 0.18558675 0.36063872 0.005305268 0.06538556 0 0.1969691 0.07829503
LAU1255 0.04230761 0.07319686 0.014003042 0.37666902 0 0.0000000 0.14757024
LAU1314 0.24241384 0.34369856 0.015805166 0.06588054 0 0.1580854 0.05774793
T cells follicular helper T cells regulatory (Tregs) T cells gamma delta NK cells resting NK cells activated Monocytes
LAU125 0.09940070 0 0.00000000 0.04715444 0.00000000 0.09273837
LAU355 0.02004924 0 0.02277827 0.00000000 0.00000000 0.01053398
LAU1255 0.11192989 0 0.04067900 0.00000000 0.04482746 0.02217635
LAU1314 0.04853605 0 0.03254491 0.00000000 0.01240244 0.00000000
Macrophages M0 Macrophages M1 Macrophages M2 Dendritic cells resting Dendritic cells activated Mast cells resting
LAU125 0.03315481 0.00000000 0.12737478 0.00000000 0.004796157 0.000000000
LAU355 0.00000000 0.00000000 0.04398640 0.00000000 0.010471658 0.000000000
LAU1255 0.02427582 0.02090543 0.06168147 0.01180571 0.000000000 0.002012619
LAU1314 0.00000000 0.00000000 0.01644681 0.00000000 0.006438391 0.000000000
Mast cells activated Eosinophils Neutrophils P-value Correlation RMSE
LAU125 0.06608013 0.004228426 0.000000000 0.5 0.01785367 1.0921730
LAU355 0.00000000 0.000000000 0.000000000 0.0 0.66799046 0.7449124
LAU1255 0.00000000 0.004148568 0.001810919 0.0 0.20761899 1.0207522
LAU1314 0.00000000 0.000000000 0.000000000 0.0 0.64988147 0.7595212
可以看出22種細胞類型在每個樣本中的分布及其顯著性指標：P-value 教硫，Correlation RMSE叨吮。
我們簡單探索一下這個返回的結(jié)果：
library(pheatmap)
pheatmap(res_cibersort[,1:22])
apply(res_cibersort[,1:22] ,1,sum)
LAU125 LAU355 LAU1255 LAU1314
1 1 1 1
可見它是按滿分是1 來對樣本打分的，也就是每個樣本只有這22種細胞類型瞬矩。既然如此茶鉴，我們就可以：
library(ggplot2)
library(tidyverse)
allcolour=c("#DC143C","#0000FF","#20B2AA","#FFA500","#9370DB","#98FB98","#F08080","#1E90FF","#7CFC00","#FFFF00","#808000","#FF00FF","#FA8072","#7B68EE","#9400D3","#800080","#A0522D","#D2B48C","#D2691E","#87CEEB","#40E0D0","#5F9EA0","#FF1493","#0000CD","#008B8B","#FFE4B5","#8A2BE2","#228B22","#E9967A","#4682B4","#32CD32","#F0E68C","#FFFFE0","#EE82EE","#FF6347","#6A5ACD","#9932CC","#8B008B","#8B4513","#DEB887")
res_cibersort[,1:22] %>% reshape2::melt() %>%
ggplot(aes(x=Var1,y=value,fill=Var2))+
geom_bar(position = 'stack',stat="identity")+
scale_fill_manual(values =allcolour ) + theme_bw()
其他工具的思路：
我們試一下這個模型在單細胞中的應(yīng)用情況。在這之前我們需要有一個已經(jīng)注釋好的數(shù)據(jù)景用，以便我們看看細胞類型能不能對應(yīng)上涵叮。于是，我們請出pbmc3k.final數(shù)據(jù)集伞插，這個是可以安裝使用的割粮。
library(Seurat)
library(SeuratData)
head(pbmc3k.final@meta.data)
orig.ident nCount_RNA nFeature_RNA seurat_annotations percent.mt RNA_snn_res.0.5 seurat_clusters
AAACATACAACCAC pbmc3k 2419 779 Memory CD4 T 3.0177759 1 1
AAACATTGAGCTAC pbmc3k 4903 1352 B 3.7935958 3 3
AAACATTGATCAGC pbmc3k 3147 1129 Memory CD4 T 0.8897363 1 1
AAACCGTGCTTCCG pbmc3k 2639 960 CD14+ Mono 1.7430845 2 2
AAACCGTGTATGCG pbmc3k 980 521 NK 1.2244898 6 6
AAACGCACTGGTAC pbmc3k 2163 781 Memory CD4 T 1.6643551 1 1
我們看到是有seurat_annotations 的細胞類型注釋的。
為了模擬bulk RNA的數(shù)據(jù)我們對亞群取平均值：
avexpr <- AverageExpression(pbmc3k.final)
因為CIBERSORT直接讀電腦文件媚污，而我們又不想先把數(shù)據(jù)輸出舀瓢，于是改一下函數(shù)的讀數(shù)據(jù)方式，再source耗美。
source("CIBERSORT/CIBERSORT.R")
cibe<-CIBERSORT(sig_matrix, cbind(rownames(avexpr$RNA), avexpr$RNA), perm=10, QN=TRUE)
cibe[1:4,1:4]
B cells naive B cells memory Plasma cells T cells CD8
Naive CD4 T 0.0004745028 0.006527557 0.00000000 0.40124293
Memory CD4 T 0.0000000000 0.017245116 0.00000000 0.14626803
CD14+ Mono 0.0057763350 0.007812051 0.00000000 0.01540697
B 0.3429285094 0.450309369 0.03076117 0.00895134
同上京髓，我們來用熱圖看看對應(yīng)關(guān)系怎么樣：
library(pheatmap)
pheatmap(cibe[,1:22])
雖然我們用的并不是腫瘤樣本，也雖然熱圖的大部分都是天藍色商架，但是橘黃色的色塊告訴我們：對應(yīng)關(guān)系還算靠譜朵锣。
cibe[,1:22] %>% reshape2::melt() %>%
ggplot(aes(x=Var1,y=value,fill=Var2))+
geom_bar(position = 'stack',stat="identity")+
scale_fill_manual(values =allcolour ) + theme_bw() +theme(axis.text.x = element_text(angle = 90, hjust=1))
這個柱形圖顯然會得到令人困惑的結(jié)果，比如NK細胞里面為什么還會有其他21種細胞呢甸私？這當(dāng)然和pbmc3k.final做細胞注釋的方法有關(guān)，也和CIBERSORT這個算法有關(guān)飞傀，他總會給一個表達譜再這22種細胞類型中找比例關(guān)系皇型，除非自己準備一個 signature matrix。但是這些結(jié)果并不是一無用處砸烦，特別是在我們對細胞類型一無所知的時候弃鸦，還是有一些參考價值的。
xCell幢痘，GSVA唬格，MCPcounter，ConsensusTME颜说，ImmuneSubtypeClassifier等等购岗，我就不再一一演示了，基本上是簡單的基因表達矩陣數(shù)據(jù)门粪，將免疫細胞類型含量計算出來喊积。不同的軟件默認給定的細胞類型不同，算法有不同玄妈，但是結(jié)果都是得到每個樣本在不同細胞類型中的分值乾吻。
在單細胞中的應(yīng)用
基于以上的討論髓梅，我們知道了所謂的定量免疫浸潤，實在無法獲得單細胞水平異質(zhì)性的歷史條件下绎签，一種通過bulk RNA數(shù)據(jù)來推斷(腫瘤)樣本免疫細胞異質(zhì)性的手段枯饿。
我們知道了做定量免疫浸潤的關(guān)鍵是有表達譜和一個參數(shù)數(shù)據(jù)集，表達譜在我們單細胞這里是不缺的诡必，關(guān)鍵在于參數(shù)數(shù)據(jù)集如何獲得奢方。定量免疫浸潤的這套方法，在單細胞中至少有以下應(yīng)用方向：
- 驗證
單細胞數(shù)據(jù)科學(xué)是一個超指數(shù)發(fā)展的行業(yè)擒权，數(shù)據(jù)是爆發(fā)式的增長的袱巨，但是數(shù)據(jù)之間并不是孤立的，之前有著豐富的數(shù)據(jù)積累碳抄，其中TCGA愉老，GEO是比較詳細的。越來越多的文章剖效，開始利用公用數(shù)據(jù)庫來挖掘有用的信息來豐富自己的研究嫉入，甚至專門有個數(shù)據(jù)挖掘的方向。在這樣的背景下璧尸，免疫浸潤也是揭示異質(zhì)性的咒林，那么自然地，我們會聯(lián)想到：拿我們的單細胞數(shù)據(jù)通過免疫浸潤與公用數(shù)據(jù)聯(lián)系起來爷光，相互驗證垫竞。
- 細胞類型推斷
既然定量免疫浸潤可以得到每個樣本在不同細胞類型中的分值，我們可不可以根據(jù)這個分值來推斷手里表達譜的細胞類型呢蛀序？在學(xué)理上應(yīng)該是可以的欢瞪，反卷積的算法也許有些不當(dāng)，但是如果把求平均值的過程看作卷積徐裸，也不是不可以吧遣鼓。基于GSVA 的方法本質(zhì)就是自定義基因集做富集分析重贺，所以完全是可以的骑祟，關(guān)鍵在與marker gene 的選擇。
dtangle:基于反卷積的表達譜分解確定細胞組分
使用xcell根據(jù)表達譜推斷樣本組成細胞的類型
3到11分文章解讀(腫瘤免疫浸潤挖掘方向)
必讀|TCGA數(shù)據(jù)挖掘-肺癌腫瘤免疫浸潤分析
一文獻一技術(shù)路線：再來“免疫浸潤”
腫瘤免疫浸潤分析工具匯總-數(shù)據(jù)挖掘新高度
TIMER:腫瘤浸潤免疫細胞分析的綜合網(wǎng)站
如何用轉(zhuǎn)錄組數(shù)據(jù)定量腫瘤浸潤免疫細胞
視頻課程：TCGA數(shù)據(jù)免疫浸潤的量化方法
癌細胞浸潤是什么意思气笙？
Targeting adenosine receptor 2B in triple negative breast cancer
Identification of key biomarkers and immune infiltration in the synovial tissue of osteoarthritis by bioinformatics analysis
Profiles of Immune Infiltration and Prognostic Immunoscore in Lung Adenocarcinoma
Elisa Martini次企， et al.Single-Cell Sequencing of Mouse Heart Immune Infiltrate in Pressure Overload–Driven Heart Failure Reveals Extent of Immune Activation.Circulation. 2019；140：2089–2107
Single-cell immune landscape of human atherosclerotic plaques
如何通俗易懂地解釋卷積健民？
Understanding-Convolutions/
量化免疫浸潤時CIBERSORT的注意事項抒巢。
evaluation-of-methods-to-assign-cell-type-labels-to-cell-clusters-from-single-cell-rna-sequencing-data
