由表達矩陣看內部異質性
- 2020 年 3 月 31 日
- 筆記
依然是主要程式碼跟著流程走，這裡只寫一寫我認為比較重要的內容
上次我們得到了原始表達矩陣a，過濾後的表達矩陣dat ，樣本資訊meata
簡單看下：
> dim(a) [1] 24582 768 > dim(dat) [1] 12198 768 > head(meta) g plate n_g all SS2_15_0048_A3 1 0048 2624 all SS2_15_0048_A6 1 0048 2664 all SS2_15_0048_A5 2 0048 3319 all SS2_15_0048_A4 3 0048 4447 all SS2_15_0048_A1 2 0048 4725 all SS2_15_0048_A2 3 0048 5263 all # 樣本資訊中的g表示cutree分的4大聚類結果；plate為細胞板批次；n_g是每個細胞樣本中有表達的基因數量；all暫時用不到
另外，注意最好每次運行程式碼之前，都要清空一下變數，然後設置不要將字元型變成因子型向量
rm(list = ls()) options(stringsAsFactors = F)
熱 · 圖
先構建分組資訊，也就是提取出層次聚類資訊
需要注意一點，count的表達矩陣和rpkm表達矩陣的聚類分組結果是不一樣的
# 我們這裡是count表達矩陣分組 > grp=meta$g > table(grp) grp 1 2 3 4 312 300 121 35
然後想想，熱圖需要什麼資訊？
主要就是行、列，行是基因，列是樣本。那麼先對基因（行）進行設置：
因為dat矩陣相對於a雖然過濾掉了一萬多基因，但是依然還剩一萬多，然後我們有700多樣本，那麼可以算一下，這樣的結果是10000*700的圖，相當大，並且看不出什麼含義。我們可以將基因數設置小一點，可以設置成1000，先看一下
好了，基因數有了（1000個），但是取哪1000個基因呢？很顯然，利用head或tail直接取前/後1000個基因是不能使人信服的，這裡可以用sd 進行篩選，也就是取表達量標準差最大的1000個基因(也即是說，這1000個基因在所有的樣本中表達差異最大，這樣更像差異表達基因)
tail(sort(apply(dat,1,sd)),1000) # 解釋下程式碼：從裡向外看=》apply對dat矩陣的每一行求sd值，然後用sort排序，默認從小到大，然後用tail從後到前，也即是從大到小取1000個 # 最後取出基因名 top_g=names(tail(sort(apply(dat,1,sd)),100)) > head(top_g) [1] "Comt" "Mrgprf" "Stard13" "Cdipt" "Ifnar1" "Pdcd6ip"
畫第一張熱圖
1000基因 X 768個細胞 = 70多萬的點
這個熱圖反映了log2(cpm+1)的表達量，範圍是0-15
library(pheatmap) pheatmap(dat[top_g,],show_colnames =F,show_rownames = F, filename = 'all_cells_top_1000_sd.png')
熱圖基礎上增加歸一化
上面?那個圖，可以看出基因絕對表達量 ，顏色越偏紅色表示絕對表達量越高，比如頂部那些基因的表達量就是要比底部那些基因的高
但是，有個問題，這樣會受到某些特高表達基因的影響，導致其他基因的差異就不明顯；另外，我們真正關心的是一個基因在不同樣本中的差異，是一種相對的表達量。
可以這麼理解：有的基因本身就是表達量小，但不能因為小就認為它在每個樣本都是一樣的。雖然小，也是有差異的小。但往往這種差異會由於"強者"的存在而被忽視。因此，這裡我們要做的，就是將這種"弱小"的差異給拎出來
主要利用scale() ，先理解一下：
# 構建一個測試數據 x=1:10;plot(x) scale(x);plot(scale(x))# 結果就是變成從-1.4到+1.4的範圍
可以看到，scale後並不改變數據分布，只是修改了坐標，讓結果取值更加集中
注意：scale是對列進行操作，而我們是想對基因(也就是按行操作)，這個函數有兩個主要的選項：center和scale ，其中center是將每列的元素減去這一列的均值(這個選項是默認TRUE的)；scale 是在center操作後，再將處理過的元素除以標準差(同樣是默認TRUE的)。另外，處理完別忘了再轉換回來
n=t(scale(t(dat[top_g,])))
畫個新的熱圖
可以看到之前熱圖顯示的坐標範圍是0- 15，當然這裡我們可以設置上限、下限，比如可以將上限設為2，下限設為-2
n[n>2]=2 n[n< -2]= -2
範圍設置好以後，可以再將分組資訊grp加上去
最後設置pheatmap的選項：
pheatmap(n, show_colnames =F, #不顯示列名 show_rownames = F, #不顯示行名 annotation_col=anno_col, # 在列上加註釋(也就是分組資訊) filename = 'all_cells_top_1000_new.png')
可以看到這張圖和畫的第一張圖的區別是：出現了一些紅色，說明新出現了一些基因在某些樣本中高表達，並且很明顯。但是仍然很有可能它們的實際表達量並不高，僅僅是玩了一個"樣本排位賽「(即使數值再小，也有甲乙丙丁)
關於分組有一點奇怪
可以看到這裡的分組資訊有點散亂，想到：這裡使用的anno_col 是利用grp得到的，而grp是基於一萬多基因的dat矩陣得到的(回憶下第5篇內容)
因此這裡的分組資訊可以更新一下，基於我們這裡的top1000基因，只需要將原來的dat換成現在的n矩陣就好，依然選取前4個聚類分群
# 將原來dat換為n hc=hclust(dist(t(n))) clus = cutree(hc, 4) top1000_grp=as.factor(clus) > table(top1000_grp) top1000_grp 1 2 3 4 462 166 42 98
看一下當前基於1000個基因的前4組和原來基於所有基因的前4組之間重合度，雖然他們總和一樣，都是1000，而且也都是按照1-4的順序排列，但實際上背後的意義千差萬別
> table(top1000_grp,grp) grp top1000_grp 1 2 3 4 1 167 295 0 0 2 7 3 121 35 3 42 0 0 0 4 96 2 0 0
舉個例子，有462個基因屬於新分組的1號組，但其中有295個屬於原來分組的2號組(這個數量超過了原來分組的1號組)，可以看出新分組和原分組的重合度並不高，因此更加說明我們重新分組的重要性
new_anno_col=data.frame(g=top1000_grp) rownames(new_anno_col)=colnames(n) pheatmap(n, show_colnames =F, #不顯示列名 show_rownames = F, #不顯示行名 annotation_col=new_anno_col, # 在列上加註釋(也就是分組資訊) filename = 'all_cells_top_1000_new_2.png')
PCA · 圖
之前好不容易過濾得到的dat矩陣，不能因為下面分析的失誤被"污染"，因此再進行下一個分析之前先做一個數據備份是個好習慣
dat_bk=dat # 然後我們就能放心對dat進行操作了 dat=t(dat) dat=as.data.frame(dat) dat=cbind(dat,grp)
PCA分析需要行是樣本，列是基因表達量的數據框(和聚類一樣，是對行/樣本進行操作，最後做的圖中一個點就表示一個樣本/細胞)
最後用PCA進行計算分析，用fviz_pca_ind函數進行可視化
這裡用到的分組還是之前基於全部基因進行聚類的cutree結果
