轉(zhuǎn)載來(lái)自：劉小澤 鏈接：http://www.reibang.com/p/e659a2f164f7
劉小澤寫(xiě)于19.6.18吓妆、23-第二單元第四講：得到表達(dá)矩陣后看內(nèi)部異質(zhì)性
筆記目的：根據(jù)生信技能樹(shù)的單細(xì)胞轉(zhuǎn)錄組課程探索smart-seq2技術(shù)相關(guān)的分析技術(shù)
課程鏈接在：http://jm.grazy.cn/index/mulitcourse/detail.html?cid=53
前言
依然是主要代碼跟著流程走赊时，這里只寫(xiě)一寫(xiě)我認(rèn)為比較重要的內(nèi)容
上次我們得到了原始表達(dá)矩陣a，過(guò)濾后的表達(dá)矩陣dat 行拢，樣本信息meata
簡(jiǎn)單看下：
> dim(a)
[1] 24582 768
> dim(dat)
[1] 12198 768
> metadata=df
> 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為細(xì)胞板批次；n_g是每個(gè)細(xì)胞樣本中有表達(dá)的基因數(shù)量剂陡；all暫時(shí)用不到
另外狈涮，注意最好每次運(yùn)行代碼之前，都要清空一下變量鸭栖，然后設(shè)置不要將字符型變成因子型向量
rm(list = ls())
options(stringsAsFactors = F)
想要做個(gè)熱圖
先構(gòu)建分組信息歌馍，也就是提取出層次聚類信息
需要注意一點(diǎn)，count的表達(dá)矩陣和rpkm表達(dá)矩陣的聚類分組結(jié)果是不一樣的
# 我們這里是count表達(dá)矩陣分組
> grp=meta$g
> table(grp)
grp
1 2 3 4
312 300 121 35
然后想想晕鹊，熱圖需要什么信息松却？
主要就是行、列溅话，行是基因晓锻，列是樣本。那么先對(duì)基因（行）進(jìn)行設(shè)置：
因?yàn)?code>dat矩陣相對(duì)于a雖然過(guò)濾掉了一萬(wàn)多基因飞几，但是依然還剩一萬(wàn)多砚哆，然后我們有700多樣本，那么可以算一下屑墨，這樣的結(jié)果是10000*700的圖躁锁，相當(dāng)大纷铣，并且看不出什么含義。我們可以將基因數(shù)設(shè)置小一點(diǎn)战转，可以設(shè)置成1000搜立，先看一下
好了，基因數(shù)有了（1000個(gè)）槐秧，但是取哪1000個(gè)基因呢啄踊？很顯然，利用head或tail直接取前/后1000個(gè)基因是不能使人信服的刁标，這里可以用sd 進(jìn)行篩選颠通，也就是取表達(dá)量標(biāo)準(zhǔn)差最大的1000個(gè)基因(也即是說(shuō)，這1000個(gè)基因在所有的樣本中表達(dá)差異最大命雀，這樣更像差異表達(dá)基因)
tail(sort(apply(dat,1,sd)),1000)
# 解釋下代碼：從里向外看=》apply對(duì)dat矩陣的每一行求sd值蒜哀，然后用sort排序，默認(rèn)從小到大吏砂，然后用tail從后到前，也即是從大到小取1000個(gè)
# 最后取出基因名
top_g=names(tail(sort(apply(dat,1,sd)),100))
> head(top_g)
[1] "Comt" "Mrgprf" "Stard13" "Cdipt" "Ifnar1" "Pdcd6ip"
畫(huà)第一張熱圖
1000基因 X 768個(gè)細(xì)胞 = 70多萬(wàn)的點(diǎn)
這個(gè)熱圖反映了log2(cpm+1)的表達(dá)量乘客，范圍是0-15
library(pheatmap)
pheatmap(dat[top_g,],show_colnames =F,show_rownames = F,
filename = 'all_cells_top_1000_sd.png')
熱圖基礎(chǔ)上增加歸一化
上面??那個(gè)圖狐血，可以看出基因絕對(duì)表達(dá)量 ，顏色越偏紅色表示絕對(duì)表達(dá)量越高易核，比如頂部那些基因的表達(dá)量就是要比底部那些基因的高
但是匈织，有個(gè)問(wèn)題，這樣會(huì)受到某些特高表達(dá)基因的影響牡直，導(dǎo)致其他基因的差異就不明顯缀匕；另外，我們真正關(guān)心的是一個(gè)基因在不同樣本中的差異碰逸，是一種相對(duì)的表達(dá)量乡小。
可以這么理解：有的基因本身就是表達(dá)量小，但不能因?yàn)樾【驼J(rèn)為它在每個(gè)樣本都是一樣的饵史。雖然小满钟，也是有差異的小。但往往這種差異會(huì)由于"強(qiáng)者"的存在而被忽視胳喷。
因此湃番，這里我們要做的，就是將這種"弱小"的差異給拎出來(lái)
主要利用scale() 吭露，先理解一下：
# 構(gòu)建一個(gè)測(cè)試數(shù)據(jù)
x=1:10;plot(x)
scale(x);plot(scale(x))# 結(jié)果就是變成從-1.4到+1.4的范圍
可以看到吠撮，scale后并不改變數(shù)據(jù)分布，只是修改了坐標(biāo)讲竿，讓結(jié)果取值更加集中
注意：scale是對(duì)列進(jìn)行操作泥兰，而我們是想對(duì)基因(也就是按行操作)择浊，這個(gè)函數(shù)有兩個(gè)主要的選項(xiàng)：center和scale ，其中center是將每列的元素減去這一列的均值(這個(gè)選項(xiàng)是默認(rèn)TRUE的)逾条；scale 是在center操作后琢岩，再將處理過(guò)的元素除以標(biāo)準(zhǔn)差(同樣是默認(rèn)TRUE的)。另外师脂，處理完別忘了再轉(zhuǎn)換回來(lái)
n=t(scale(t(dat[top_g,])))
畫(huà)個(gè)新的熱圖
可以看到之前熱圖顯示的坐標(biāo)范圍是0-15担孔，當(dāng)然這里我們可以設(shè)置上限、下限吃警，比如可以將上限設(shè)為2糕篇，下限設(shè)為-2
n[n>2]=2
n[n< -2]= -2
范圍設(shè)置好以后，可以再將分組信息grp加上去
# 先構(gòu)建一個(gè)分組的數(shù)據(jù)框酌心，列是原來(lái)的分組信息
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的選項(xiàng)：
pheatmap(n,
show_colnames =F, #不顯示列名
show_rownames = F, #不顯示行名
annotation_col=anno_col, # 在列上加注釋(也就是分組信息)
filename = 'all_cells_top_1000_new.png')
可以看到這張圖和畫(huà)的第一張圖的區(qū)別是：出現(xiàn)了一些紅色，說(shuō)明新出現(xiàn)了一些基因在某些樣本中高表達(dá)安券，并且很明顯墩崩。但是仍然很有可能它們的實(shí)際表達(dá)量并不高，僅僅是玩了一個(gè)"樣本排位賽“(即使數(shù)值再小侯勉，也有甲乙丙丁)
關(guān)于分組有一點(diǎn)奇怪
可以看到這里的分組信息有點(diǎn)散亂鹦筹，想到：這里使用的anno_col 是利用grp得到的，而grp是基于一萬(wàn)多基因的dat矩陣得到的(回憶下第5篇內(nèi)容)
因此這里的分組信息可以更新一下址貌，基于我們這里的top1000基因铐拐，只需要將原來(lái)的dat換成現(xiàn)在的n矩陣就好，依然選取前4個(gè)聚類分群
# 將原來(lái)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個(gè)基因的前4組和原來(lái)基于所有基因的前4組之間重合度练对，雖然他們總和一樣遍蟋，都是1000，而且也都是按照1-4的順序排列螟凭，但實(shí)際上背后的意義千差萬(wàn)別
> 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
舉個(gè)例子虚青，有462個(gè)基因?qū)儆谛路纸M的1號(hào)組，但其中有295個(gè)屬于原來(lái)分組的2號(hào)組(這個(gè)數(shù)量超過(guò)了原來(lái)分組的1號(hào)組)赂摆，可以看出新分組和原分組的重合度并不高挟憔，因此更加說(shuō)明我們重新分組的重要性
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')
做個(gè)PCA
之前好不容易過(guò)濾得到的dat矩陣，不能因?yàn)橄旅娣治龅氖д`被"污染"烟号，因此再進(jìn)行下一個(gè)分析之前先做一個(gè)數(shù)據(jù)備份是個(gè)好習(xí)慣
dat_bk=dat
# 然后我們就能放心對(duì)dat進(jìn)行操作了
dat=t(dat)
dat=as.data.frame(dat)
dat=cbind(dat,grp)
#列是基因绊谭，行是樣本名，最后一列加上其分組
> dat[1:4,1:4]
0610007P14Rik 0610009B22Rik 0610009L18Rik 0610009O20Rik
SS2_15_0048_A3 0.000000 0 0 0.000000
SS2_15_0048_A6 0.000000 0 0 0.000000
SS2_15_0048_A5 6.459884 0 0 2.544699
SS2_15_0048_A4 6.313884 0 0 3.025273
> ## 表達(dá)矩陣可以隨心所欲的取行列汪拥，基礎(chǔ)知識(shí)需要打牢达传。
> dat[1:4,12197:12199]
ERCC-00170 ERCC-00171 group_list
SS2_15_0048_A3 0.000000 12.30686 1
SS2_15_0048_A6 7.262888 12.18923 1
SS2_15_0048_A5 0.000000 11.68331 1
SS2_15_0048_A4 5.672132 10.68257 2
PCA分析需要行是樣本，列是基因表達(dá)量的數(shù)據(jù)框(和聚類一樣，是對(duì)行/樣本進(jìn)行操作宪赶，最后做的圖中一個(gè)點(diǎn)就表示一個(gè)樣本/細(xì)胞)
最后用PCA進(jìn)行計(jì)算分析宗弯，用fviz_pca_ind函數(shù)進(jìn)行可視化
這里用到的分組還是之前基于全部基因進(jìn)行聚類的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"
)
作者：劉小澤
鏈接：http://www.reibang.com/p/2b16c40088db
來(lái)源：簡(jiǎn)書(shū)
著作權(quán)歸作者所有。商業(yè)轉(zhuǎn)載請(qǐng)聯(lián)系作者獲得授權(quán)搂妻，非商業(yè)轉(zhuǎn)載請(qǐng)注明出處蒙保。
