轉(zhuǎn)載來(lái)自：劉小澤 鏈接：http://www.reibang.com/p/e659a2f164f7
劉小澤寫(xiě)于19.6.17-第二單元第三講：熟悉文獻(xiàn)作者提供的兩個(gè)表達(dá)矩陣
筆記目的：根據(jù)生信技能樹(shù)的單細(xì)胞轉(zhuǎn)錄組課程探索smart-seq2技術(shù)相關(guān)的分析技術(shù)
鏈接在：https://www.bilibili.com/video/BV1dt411Y7nn?p=6
過(guò)濾后的操作
上次得到的dat表達(dá)矩陣過(guò)濾掉低表達(dá)基因后，剩下12198個(gè)基因
看看其中的spike-in情況
grep('ERCC',rownames(dat))
[1] 12139 12140 12141 12142 12143 12144 12145 12146 12147 12148 12149 12150
[13] 12151 12152 12153 12154 12155 12156 12157 12158 12159 12160 12161 12162
[25] 12163 12164 12165 12166 12167 12168 12169 12170 12171 12172 12173 12174
[37] 12175 12176 12177 12178 12179 12180 12181 12182 12183 12184 12185 12186
[49] 12187 12188 12189 12190 12191 12192 12193 12194 12195 12196 12197 12198
一、去除細(xì)胞文庫(kù)大小差異 (歸一化 cpm法)
每個(gè)細(xì)胞測(cè)得數(shù)據(jù)大小不同轨奄，這樣是沒(méi)辦法看高表達(dá)還是低表達(dá)的芜飘，必須先保證基數(shù)一樣才能比較，cpm（counts per million）這個(gè)算法就是做這個(gè)事情的埠对。
cpm是歸一化的一種方法下愈，代表每百萬(wàn)堿基中每個(gè)轉(zhuǎn)錄本的count值
注意：這個(gè)算法只是校正文庫(kù)差異，而沒(méi)有校正基因長(zhǎng)度差異要尔。要注意我們分析的目的就是：比較一個(gè)基因在不同細(xì)胞的表達(dá)量差異，而不是考慮一個(gè)樣本中不同兩個(gè)基因的差異新娜，因?yàn)?沒(méi)有兩片相同的樹(shù)葉”這個(gè)差異是正常的赵辕。但是同一個(gè)基因由于某種條件發(fā)生了改變，背后的生物學(xué)意義是更值得探索的概龄。
用起來(lái)很簡(jiǎn)單还惠，有現(xiàn)成的函數(shù)cpm() ，然后我們?cè)儆胠og將數(shù)據(jù)降個(gè)維度私杜，但保持原有數(shù)據(jù)形狀不變：log2(edgeR::cpm(dat)+1)
意思就是：cpm需要除以測(cè)序總reads數(shù)蚕键，而這個(gè)值作為分母會(huì)導(dǎo)致結(jié)果千差萬(wàn)別，有的特別大有的很小衰粹。為了后面可視化不受極值的影響锣光，用log轉(zhuǎn)換一下可以將數(shù)值變小，并且原來(lái)大的數(shù)值最后還是大铝耻，并不改變這個(gè)現(xiàn)實(shí)
那么具體這個(gè)函數(shù)做了什么事誊爹，才是真正需要了解的：
# 先看看前4行4列的數(shù)據(jù)
> dat[1:4,1:4]
SS2_15_0048_A3 SS2_15_0048_A6 SS2_15_0048_A5 SS2_15_0048_A4
0610007P14Rik 0 0 18 11
0610009B22Rik 0 0 0 0
0610009L18Rik 0 0 0 0
0610009O20Rik 0 0 1 1
# 比如先計(jì)算一下第三個(gè)樣本的總測(cè)序量
> sum(dat[,3])
[1] 206831 #結(jié)果是0.2M
# 那么對(duì)于第三個(gè)樣本SS2_15_0048_A5的第一個(gè)基因0610007P14Rik(結(jié)果是18)
# 計(jì)算它的cpm值：count值*1000000/總測(cè)序reads
> 18*1000000/206831
[1] 87.02757
# 和標(biāo)準(zhǔn)公式比較看看，結(jié)果完全相同
> edgeR::cpm(dat[,3])[1]
[1] 87.02757
# 因此最后就是
dat=log2(edgeR::cpm(dat)+1)
二、歸一化后聚類(lèi)
主要步驟為：
1频丘、用dist函數(shù)計(jì)算細(xì)胞直接的距離办成，用hclus進(jìn)行層次聚類(lèi)，用cutree提取分群的數(shù)目（比如分成4個(gè)cluster）搂漠。注釋到每個(gè)細(xì)胞所在的cluster數(shù)迂卢。
2、用strsplit提取每個(gè)細(xì)胞的板號(hào)状答，用apply計(jì)算每個(gè)細(xì)胞的基因表達(dá)總數(shù)冷守。
第一步：理解dist函數(shù)
首先理解，它是計(jì)算距離用的惊科，正如函數(shù)名稱所描述的一樣：distance
# 先構(gòu)建一個(gè)測(cè)試矩陣
x=1:5
y=2*x
z=52:56
tmp=data.frame(x,y,z)
> tmp
x y z
1 1 2 52
2 2 4 53
3 3 6 54
4 4 8 55
5 5 10 56
# 可以看到拍摇，x和y是有一定相關(guān)性的，而z和它們很難扯上關(guān)系
# 然后嘗試計(jì)算x馆截、y充活、z之間的距離，來(lái)驗(yàn)證我們的猜想
> dist(tmp)
1 2 3 4
2 2.449490
3 4.898979 2.449490
4 7.348469 4.898979 2.449490
5 9.797959 7.348469 4.898979 2.449490
# 好像得到的不是我們想要的蜡娶。我們想要的是x混卵、y、z距離結(jié)果窖张，而計(jì)算給出的是以"行"為單位的結(jié)果
# 因此幕随，猜測(cè)dist應(yīng)該是以行為輸入。因此修改一下tmp宿接，讓x赘淮、y、z為行睦霎，其實(shí)也就是轉(zhuǎn)置一下梢卸，轉(zhuǎn)置函數(shù)用t()
> dist(t(tmp))
x y
y 7.416198
z 114.039467 107.377838
> cor(tmp)
x y z
x 1.00000000 1.00000000 -0.09844392
y 1.00000000 1.00000000 -0.09844392
z -0.09844392 -0.09844392 1.00000000
> dist(t(scale(tmp)))
x y
y 0.000000
z 4.446571 4.446571
# 可以看到dist函數(shù)計(jì)算樣本直接距離和cor函數(shù)計(jì)算樣本直接相關(guān)性，是完全不同的概念副女。雖然我都沒(méi)有調(diào)它們兩個(gè)函數(shù)的默認(rèn)的參數(shù)蛤高。
#
# 總結(jié)：
#
# - dist函數(shù)計(jì)算行與行（樣本）之間的距離
# - cor函數(shù)計(jì)算列與列（樣本）之間的相關(guān)性
# - scale函數(shù)默認(rèn)對(duì)每一列（樣本）內(nèi)部歸一化
# - 計(jì)算dist之前，最好是對(duì)每一個(gè)樣本（列）進(jìn)行scale一下
#
同樣的碑幅，我們這里的dat數(shù)據(jù)戴陡，是要計(jì)算細(xì)胞間的距離，也就是列與列之間的距離沟涨，使用dist(t(dat)) 計(jì)算猜欺。數(shù)據(jù)中有768個(gè)細(xì)胞，也就是要計(jì)算768個(gè)細(xì)胞核768個(gè)細(xì)胞之間的距離拷窜，計(jì)算量還是很大的。
關(guān)于dist計(jì)算距離的方法：主要有6種：”歐式euclidean”, “切比雪夫距離maximum”, “絕對(duì)值距離manhattan”, “Lance距離canberra”, “定型變量距離binary” or “明可夫斯基距離minkowski（使用時(shí)要指定p值）”。
默認(rèn)使用第一種歐氏距離篮昧，它計(jì)算的是：幾何空間中兩點(diǎn)之間的距離赋荆。思想類(lèi)似于勾股定理求第三條斜邊的長(zhǎng)度=》平方和再開(kāi)方。
第二步：理解hclust函數(shù)
它是進(jìn)行層次聚類(lèi)（系譜聚類(lèi)）的方法
關(guān)于hclust聚類(lèi)的方法：”離差平方和法ward”, “最短距離法single”, “最長(zhǎng)距離法complete”,”類(lèi)平均法average”, “相似法mcquitty”, “中間距離法median” or “重心法centroid”懊昨。默認(rèn)使用complete算法
hc=hclust(dist(t(dat)))
# 如果要進(jìn)行可視化
plot(hc,labels = FALSE) #labels這個(gè)選項(xiàng)的意思是不顯示各個(gè)樣本名稱窄潭，因?yàn)闃颖咎啵瑫?huì)讓圖看起來(lái)很亂
另外hclust函數(shù)還有一個(gè)親戚：cutree酵颁，顧名思義嫉你，就是對(duì)聚類(lèi)樹(shù)進(jìn)行修剪。我們知道聚類(lèi)結(jié)果是分群的躏惋，cutree就是指定輸出哪些群(結(jié)果是從大群到小群排列)
# 例如要看看分的4大群
clus = cutree(hc, 4)
group_list= as.factor(clus) #得到的這個(gè)因子型變量group_list中樣本順序和輸入的順序一致幽污，并且屬于第幾類(lèi)都有記錄
> table(group_list)
group_list
1 2 3 4
312 300 121 35
提取批次信息
在上一步操作結(jié)果中，可以看到簿姨，樣本名都是有規(guī)律的距误，例如：
> head(colnames(dat))
[1] "SS2_15_0048_A3" "SS2_15_0048_A6" "SS2_15_0048_A5" "SS2_15_0048_A4"
[5] "SS2_15_0048_A1" "SS2_15_0048_A2"
其中SS2_15都是一樣的，Pxx也不需要管扁位，重要的是中間的0048准潭、0049，表示兩個(gè)384孔板編號(hào)
那么如何提扔虺稹刑然？
使用strsplit函數(shù)，strsplit(x, split, fixed = FALSE) 暇务，需要注意兩點(diǎn)：
字符串切分后泼掠，返回的是一個(gè)列表，如果要再還原成字符串般卑，需要用
unlist()-
默認(rèn)情況下它是使用正則表達(dá)式的武鲁，如果不想用，可以指定
fixed = TRUE> unlist(strsplit("a.b.c", ".")) [1] "" "" "" "" "" > unlist(strsplit("a.b.c", ".", fixed = TRUE)) [1] "a" "b" "c"
# 方法一：純base包(思路就是：將拆分得到的list變成數(shù)據(jù)框)
options(stringsAsFactors = F)
plate=do.call(rbind.data.frame,strsplit(colnames(dat),"_"))[,3]
# 方法二：stringr包
library(stringr)
plate=str_split(colnames(dat),'_',simplify = T)[,3]
統(tǒng)計(jì)每個(gè)樣本的基因表達(dá)數(shù)目
主要使用cutree剪下來(lái)的層次聚類(lèi)信息蝠检、細(xì)胞板批次信息沐鼠、每個(gè)樣本的基因表達(dá)信息
前兩個(gè)已經(jīng)具備，下面進(jìn)行第三個(gè)：每個(gè)樣本的基因表達(dá)信息
# 還記得之前對(duì)基因進(jìn)行過(guò)濾時(shí)叹谁，我們是對(duì)行進(jìn)行操作
apply(a,1, function(x) sum(x>1) > floor(ncol(a)/50))
# 這里檢測(cè)每個(gè)樣本中有多少基因是表達(dá)的饲梭，count值以1為標(biāo)準(zhǔn)，rpkm值可以用0為標(biāo)準(zhǔn)
n_g = apply(a,2,function(x) sum(x>1))
# 對(duì)于單細(xì)胞轉(zhuǎn)錄組焰檩，一般會(huì)有超過(guò)半數(shù)的基因不會(huì)表達(dá)(這個(gè)在下面構(gòu)建完數(shù)據(jù)框還可以再看一下)
最后新建細(xì)胞的屬性信息
可以構(gòu)建數(shù)據(jù)框了：
meta=data.frame(g=group_list,plate=plate,n_g=n_g)
# 然后再添加一列憔涉，目前用不到，后續(xù)會(huì)介紹
meta$all='all'
可以看到細(xì)胞中檢測(cè)到表達(dá)的基因最多有7372個(gè)析苫，最少才幾十個(gè)兜叨，而我們總共有12000多個(gè)基因
作者：劉小澤
鏈接：http://www.reibang.com/p/33a7eb26bd31
來(lái)源：簡(jiǎn)書(shū)
著作權(quán)歸作者所有穿扳。商業(yè)轉(zhuǎn)載請(qǐng)聯(lián)系作者獲得授權(quán)，非商業(yè)轉(zhuǎn)載請(qǐng)注明出處国旷。
