“清晰而又吸引人——這無疑是學習R的有趣方式郁副!”
???????????????????????????????????????????????????????????????——Amos A. Folarin,倫敦大學學院
1.對readscount進行l(wèi)n()轉(zhuǎn)換
表達矩陣的行代表feature(如基因绳锅、外顯子等)耕肩,列代表sample;對于原始的每一個基因每一個樣本對應(yīng)的reads數(shù)進行以e為底的轉(zhuǎn)換:ln(counts)ln()函數(shù)在R語言中是log()
log(counts)
2.對每一行的值進行均值計算
對該基因?qū)?yīng)的所有樣本的轉(zhuǎn)換后的readscount計算均值
rowMeans(log(counts))
3.ln(counts)減去每一行的均值
log(counts)-rowMeans(log(counts))
4.篩選出step3得到的矩陣中不包含-Inf董朝,NaN的基因行
(log(counts)-rowMeans(log(counts)))[is.finite(rowMeans(log(counts))),]
5.計算step4得到的矩陣中每列的中位數(shù)
median((log(counts)-rowMeans(log(counts)))[is.finite(rowMeans(log(counts))),])
6.計算每列的e中位數(shù),該數(shù)即為sizefactor
exp(median((log(counts)-rowMeans(log(counts)))[is.finite(rowMeans(log(counts))),]))
7.原始每一個readscount除以對應(yīng)列的sizefactor即為歸一化后的基因表達量
original readscount/size factors
即為最終的每個樣本每個基因?qū)?yīng)的歸一化的表達量干跛。每列中的任意original readscount除以該列的sizefactor子姜。
以上代碼只具有表面含義,因為計算時各種數(shù)據(jù)類型的存在,實際在R中運算較為復雜哥捕。當然牧抽,以上步驟只是用DESeq2去estimating size factors的內(nèi)部算法,實際可以直接利用sizeFactors(dds)
去計算size factors遥赚。