轉(zhuǎn)載來(lái)自：劉小澤 鏈接：http://www.reibang.com/p/e659a2f164f7

劉小澤寫于19.7.4-第二單元第五講：重復(fù)平均表達(dá)量和變異系數(shù)相關(guān)性散點(diǎn)圖

筆記目的：根據(jù)生信技能樹的單細(xì)胞轉(zhuǎn)錄組課程探索smart-seq2技術(shù)相關(guān)的分析技術(shù)
課程鏈接在：https://www.bilibili.com/video/BV1dt411Y7nn?p=8

前言

這一次的目的是重復(fù)文章附件中的一張圖：

附件地址在：https://static-content.springer.com/esm/art%3A10.1038%2Fs41467-018-07582-3/MediaObjects/41467_2018_7582_MOESM1_ESM.pdf

存在這樣一張圖：

image

圖片分析

首先看橫坐標(biāo)咪惠，不論是RPKM還是原始count都是表達(dá)值让网，然后做了均值的log10處理诞丽；然后縱坐標(biāo)是CV值磅轻，它表示變異系數(shù)(coefficient of variation)，也是先求CV的平方吵取，然后做log10處理

來(lái)看一下什么是CV值禽额？

變異系數(shù)又稱離散系數(shù)或相對(duì)偏差(我們肯定都聽過(guò)標(biāo)準(zhǔn)偏差，也就是sd值皮官，它描述了數(shù)據(jù)值偏離算術(shù)平均值的程度)脯倒，這個(gè)相對(duì)偏差描述的是標(biāo)準(zhǔn)偏差與平均值之比，即：cv=sd/mean*100% 捺氢。

為何不用sd而用cv值呢藻丢？

先說(shuō)說(shuō)sd值，它和均值mean摄乒、方差var一樣悠反，都是對(duì)一維數(shù)據(jù)進(jìn)行的分析，需要數(shù)據(jù)滿足兩個(gè)條件：中部馍佑、單峰斋否。也就是說(shuō)數(shù)據(jù)集只存在一個(gè)峰值，并且這個(gè)峰值大致位于數(shù)據(jù)集的中部拭荤。另外當(dāng)比較兩組數(shù)據(jù)集的離散程度大小時(shí)茵臭，即使它們各自滿足"中部單峰"的條件，如果出現(xiàn)兩組數(shù)據(jù)測(cè)量尺度差別太大或數(shù)據(jù)量綱存在差異的話舅世，直接用標(biāo)準(zhǔn)差就不合適了

變異系數(shù)就可以解決這個(gè)問(wèn)題旦委，它利用原始數(shù)據(jù)標(biāo)準(zhǔn)差和原始數(shù)據(jù)平均值的比值來(lái)各自消除尺度與量綱的差異奇徒。

CV的意義

https://www.biostars.org/p/272801/

Genes which are stably expressed across replicates/experiments as the CV is low (0.5).

代碼計(jì)算

第一步：清空變量，加載數(shù)據(jù)

rm(list = ls())  ## 魔幻操作缨硝，一鍵清空~
options(stringsAsFactors = F)
load(file = '../input_rpkm.Rdata')
# 就會(huì)得到RPKM的過(guò)濾表達(dá)矩陣dat和細(xì)胞信息metadata(包括clust分群摩钙、細(xì)胞板信息扶镀、檢測(cè)到的基因數(shù)量)
> head(metadata)
               g plate  n_g all
SS2_15_0048_A3 1  0048 3065 all
SS2_15_0048_A6 2  0048 3036 all
SS2_15_0048_A5 1  0048 3742 all
SS2_15_0048_A4 3  0048 5012 all
SS2_15_0048_A1 1  0048 5126 all
SS2_15_0048_A2 3  0048 5692 all

第二步：計(jì)算mean、sd值

mean_per_gene <- apply(dat, 1, mean, na.rm = TRUE) #對(duì)表達(dá)矩陣每行求均值
sd_per_gene <- apply(dat, 1, sd, na.rm = TRUE) #同理求標(biāo)準(zhǔn)差

第三步：構(gòu)建數(shù)據(jù)框羔挡，計(jì)算cv值

cv_per_gene <- data.frame(mean = mean_per_gene,
  sd = sd_per_gene,
  cv = sd_per_gene/mean_per_gene)
# 然后賦予行名
rownames(cv_per_gene) <- rownames(dat)
> head(cv_per_gene)
                   mean       sd       cv
0610007P14Rik 25.634318 48.60264 1.895999
0610009B22Rik 27.203266 64.41918 2.368068
0610009L18Rik  3.864324 19.61355 5.075544
0610009O20Rik 10.808251 26.01667 2.407112
0610010F05Rik  8.194137 21.06394 2.570612
0610010K14Rik 31.812982 52.29101 1.643701

作圖敲才，一點(diǎn)點(diǎn)優(yōu)化

先畫第一張圖—縱坐標(biāo)為cv

with(cv_per_gene,plot(log10(mean),log10(cv)))
# with實(shí)現(xiàn)的功能就是將逗號(hào)后面的操作對(duì)象限定在逗號(hào)前面的數(shù)據(jù)框中，就不用每次都引用數(shù)據(jù)框了

image

然后畫第二張圖—縱坐標(biāo)為cv平方

with(cv_per_gene,plot(log10(mean),log10(cv^2)))

image

發(fā)現(xiàn)縱坐標(biāo)的區(qū)間發(fā)生改變橘霎，點(diǎn)的位置沒(méi)有改變，好了，初見雛形谬返，和原圖最大的差別是趨勢(shì)線

然后畫第三張圖—添加趨勢(shì)線

為了更方便地模擬原始數(shù)據(jù)，先在CV的數(shù)據(jù)框中添加兩列：log10cv2和log10mean日杈；另外看到我們上面作的圖中x軸范圍是-1~5遣铝，而文章的圖中是0~4

cv_per_gene$log10cv2=log10(cv_per_gene$cv^2)
cv_per_gene$log10mean=log10(cv_per_gene$mean)
# 修改x軸坐標(biāo)范圍為0-4
cv_per_gene=cv_per_gene[cv_per_gene$log10mean < 4, ]
cv_per_gene=cv_per_gene[cv_per_gene$log10mean > 0, ]

library(ggpubr)
ggscatter(cv_per_gene, x = 'log10mean', y = 'log10cv2',
          color = "black", shape = 16, size = 1, # Points color, shape and size
          xlab = 'log10(mean RPKM)', ylab = "log10(cv^2)",
          add = "loess", #添加擬合曲線
          add.params = list(color = "red",fill = "lightgray"),
          cor.coeff.args = list(method = "spearman"), 
          label.x = 3,label.sep = "\n",
          dot.size = 2,
          ylim=c(-0.5, 2.5),
          xlim=c(0,4) 
)

image

發(fā)現(xiàn)和原文的差別還是蠻大的，為什么呢莉擒？趨勢(shì)線添加loess不應(yīng)該有這么明顯的上升趨勢(shì)才對(duì)酿炸，想到原文的注釋可能是最好的答案，又返回看了原文涨冀，發(fā)現(xiàn)這么一句話

image

原來(lái)趨勢(shì)線不是指全部的基因填硕，而是ERCC spike-in

> grep("ERCC",rownames(cv_per_gene))
 [1] 11475 11476 11477 11478 11479 11480 11481 11482 11483 11484 11485
[12] 11486 11487 11488 11489 11490 11491 11492 11493 11494 11495 11496
[23] 11497 11498 11499 11500 11501 11502 11503 11504 11505 11506 11507
[34] 11508 11509 11510 11511 11512 11513 11514 11515 11516 11517 11518
[45] 11519 11520 11521 11522 11523 11524 11525 11526 11527 11528 11529
[56] 11530
> tail(rownames(cv_per_gene))
[1] "ERCC-00160" "ERCC-00162" "ERCC-00163" "ERCC-00165" "ERCC-00170"
[6] "ERCC-00171"

于是需要做兩個(gè)圖層，其中一個(gè)需要挑出來(lái)ERCC鹿鳖，然后單獨(dú)做扁眯，然后搜索了ggscatter的幫助信息，發(fā)現(xiàn)只需要增加一個(gè)ggp參數(shù)就可以做到

image

ggscatter(cv_per_gene[-grep("ERCC",rownames(cv_per_gene)),], x = 'log10mean', y = 'log10cv2',
          color = "black", shape = 16, size = 1, # Points color, shape and size
          xlab = 'log10(mean RPKM)', ylab = "log10(CV^2)",
          mean.point=T,
          cor.coeff.args = list(method = "spearman"), 
          label.x = 3,label.sep = "\n",
          dot.size = 2,
          ylim=c(-0.5, 2.5),
          xlim=c(0,4),
          # 這里ggp參數(shù)的意思就是：增加一個(gè)ggplot圖層
          ggp = ggscatter(cv_per_gene[grep("ERCC",rownames(cv_per_gene)),], x = 'log10mean', y = 'log10cv2',
                          color = "red", shape = 16, size = 1, # Points color, shape and size
                          xlab = 'log10(mean RPKM)', ylab = "log10(CV^2)",
                          add = "loess", #添加擬合曲線
                          mean.point=T,
                          add.params = list(color = "red",fill = "lightgray"),
                          cor.coeff.args = list(method = "pearson"), 
                          label.x = 3,label.sep = "\n",
                          dot.size = 2,
          )
)

image

寫寫自己的理解

再看原文翅帜，發(fā)現(xiàn)這個(gè)附圖2.a主要是用來(lái)說(shuō)明ERCC的姻檀，也就是做了一個(gè)技術(shù)誤差檢測(cè)，變異系數(shù)對(duì)兩個(gè)或多個(gè)數(shù)據(jù)集進(jìn)行比較時(shí)涝滴，如果度量單位保持平均數(shù)相同绣版，那么可以直接比較標(biāo)準(zhǔn)差，也就是說(shuō)歼疮，設(shè)定同一個(gè)橫坐標(biāo)(比如1)杂抽，然后紅線是ERCC，衡量技術(shù)誤差腋妙，如果我們測(cè)得基因在ERCC以下默怨，說(shuō)明我們測(cè)得基因sd值小于ERCC標(biāo)準(zhǔn)的，說(shuō)明基因的技術(shù)誤差也是在可接受范圍之內(nèi)

作者：劉小澤
鏈接：http://www.reibang.com/p/3525e624946a
來(lái)源：簡(jiǎn)書
著作權(quán)歸作者所有骤素。商業(yè)轉(zhuǎn)載請(qǐng)聯(lián)系作者獲得授權(quán)匙睹，非商業(yè)轉(zhuǎn)載請(qǐng)注明出處愚屁。