劉小澤寫于19.6.18唆涝、23-第二單元第四講：得到表達矩陣后看內(nèi)部異質(zhì)性
筆記目的：根據(jù)生信技能樹的單細胞轉(zhuǎn)錄組課程探索smart-seq2技術(shù)相關(guān)的分析技術(shù)
課程鏈接在：http://jm.grazy.cn/index/mulitcourse/detail.html?cid=53
前言
依然是主要代碼跟著流程走宛篇，這里只寫一寫我認為比較重要的內(nèi)容
上次我們得到了原始表達矩陣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大聚類結(jié)果憾儒；plate為細胞板批次势誊；n_g是每個細胞樣本中有表達的基因數(shù)量苗傅；all暫時用不到
另外，注意最好每次運行代碼之前缸浦，都要清空一下變量夕冲，然后設(shè)置不要將字符型變成因子型向量
rm(list = ls())
options(stringsAsFactors = F)
想要做個熱圖
先構(gòu)建分組信息，也就是提取出層次聚類信息
需要注意一點裂逐，count的表達矩陣和rpkm表達矩陣的聚類分組結(jié)果是不一樣的
# 我們這里是count表達矩陣分組
> grp=meta$g
> table(grp)
grp
1 2 3 4
312 300 121 35
然后想想歹鱼，熱圖需要什么信息？
主要就是行卜高、列弥姻，行是基因，列是樣本掺涛。那么先對基因（行）進行設(shè)置：
因為dat矩陣相對于a雖然過濾掉了一萬多基因庭敦，但是依然還剩一萬多，然后我們有700多樣本薪缆，那么可以算一下秧廉，這樣的結(jié)果是10000*700的圖，相當(dāng)大矮燎，并且看不出什么含義定血。我們可以將基因數(shù)設(shè)置小一點，可以設(shè)置成1000诞外，先看一下
好了澜沟，基因數(shù)有了（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')
熱圖基礎(chǔ)上增加歸一化
上面??那個圖主经，可以看出基因絕對表達量 ，顏色越偏紅色表示絕對表達量越高庭惜，比如頂部那些基因的表達量就是要比底部那些基因的高
但是罩驻，有個問題，這樣會受到某些特高表達基因的影響护赊，導(dǎo)致其他基因的差異就不明顯惠遏；另外，我們真正關(guān)心的是一個基因在不同樣本中的差異骏啰，是一種相對的表達量节吮。
可以這么理解：有的基因本身就是表達量小，但不能因為小就認為它在每個樣本都是一樣的器一。雖然小课锌，也是有差異的小。但往往這種差異會由于"強者"的存在而被忽視祈秕。
因此，這里我們要做的雏胃，就是將這種"弱小"的差異給拎出來
主要利用scale() 请毛，先理解一下：
# 構(gòu)建一個測試數(shù)據(jù)
x=1:10;plot(x)
scale(x);plot(scale(x))# 結(jié)果就是變成從-1.4到+1.4的范圍
可以看到，scale后并不改變數(shù)據(jù)分布瞭亮，只是修改了坐標方仿，讓結(jié)果取值更加集中
注意：scale是對列進行操作，而我們是想對基因(也就是按行操作)统翩，這個函數(shù)有兩個主要的選項：center和scale 仙蚜，其中center是將每列的元素減去這一列的均值(這個選項是默認TRUE的)；scale 是在center操作后厂汗，再將處理過的元素除以標準差(同樣是默認TRUE的)委粉。另外，處理完別忘了再轉(zhuǎn)換回來
n=t(scale(t(dat[top_g,])))
畫個新的熱圖
可以看到之前熱圖顯示的坐標范圍是0-15娶桦，當(dāng)然這里我們可以設(shè)置上限贾节、下限汁汗，比如可以將上限設(shè)為2，下限設(shè)為-2
n[n>2]=2
n[n< -2]= -2
范圍設(shè)置好以后栗涂，可以再將分組信息grp加上去
# 先構(gòu)建一個分組的數(shù)據(jù)框知牌，列是原來的分組信息
anno_col=data.frame(g=grp)
# 行名是樣本名，也就是歸一化后的n矩陣的列名
rownames(anno_col)=colnames(n)
# 看一下結(jié)果
> head(anno_col)
g
SS2_15_0048_A3 1
SS2_15_0048_A6 1
SS2_15_0048_A5 2
SS2_15_0048_A4 3
SS2_15_0048_A1 2
SS2_15_0048_A2 3
最后設(shè)置pheatmap的選項：
pheatmap(n,
show_colnames =F, #不顯示列名
show_rownames = F, #不顯示行名
annotation_col=anno_col, # 在列上加注釋(也就是分組信息)
filename = 'all_cells_top_1000_new.png')
可以看到這張圖和畫的第一張圖的區(qū)別是：出現(xiàn)了一些紅色斤程，說明新出現(xiàn)了一些基因在某些樣本中高表達角寸，并且很明顯。但是仍然很有可能它們的實際表達量并不高忿墅，僅僅是玩了一個"樣本排位賽“(即使數(shù)值再小扁藕，也有甲乙丙丁)
關(guān)于分組有一點奇怪
可以看到這里的分組信息有點散亂，想到：這里使用的anno_col 是利用grp得到的球匕，而grp是基于一萬多基因的dat矩陣得到的(回憶下第5篇內(nèi)容)
因此這里的分組信息可以更新一下纹磺，基于我們這里的top1000基因，只需要將原來的dat換成現(xiàn)在的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
看一下當(dāng)前基于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個基因?qū)儆谛路纸M的1號組役耕，但其中有295個屬于原來分組的2號組(這個數(shù)量超過了原來分組的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矩陣瞬痘，不能因為下面分析的失誤被"污染"故慈，因此再進行下一個分析之前先做一個數(shù)據(jù)備份是個好習(xí)慣
dat_bk=dat
# 然后我們就能放心對dat進行操作了
dat=t(dat)
dat=as.data.frame(dat)
dat=cbind(dat,grp)
PCA分析需要行是樣本，列是基因表達量的數(shù)據(jù)框(和聚類一樣框全，是對行/樣本進行操作察绷，最后做的圖中一個點就表示一個樣本/細胞)
最后用PCA進行計算分析，用fviz_pca_ind函數(shù)進行可視化
這里用到的分組還是之前基于全部基因進行聚類的cutree結(jié)果
library("FactoMineR")
library("factoextra")
dat.pca <- PCA(dat[,-ncol(dat)], graph = FALSE) #'-'表示“非”
fviz_pca_ind(dat.pca,repel =T,
geom.ind = "point",
col.ind = dat$grp, # 按組分顏色
# palette = c("#00AFBB", "#E7B800"),
addEllipses = TRUE,
legend.title = "Groups"
)
