劉小澤寫于18.11.12 半年之前我和花花成立了“生信星球”，轉眼過得好快，我們每天都能和愛生信的小伙伴們交流，還辦了9期培訓，灰常開心
終于可以回歸寫推送了，連續(xù)半個月奮戰(zhàn)兩篇文章，好累，真不如學知識寫推送輕松。時隔半月，有點忘記，但是重新學一次，又能學到新知識，還是不錯的
數(shù)據(jù)準備:
數(shù)據(jù)地址在 http://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE63310&format=file
files <- c("GSE63310_RAW/GSM1545535_10_6_5_11.txt",
"GSE63310_RAW/GSM1545536_9_6_5_11.txt",
"GSE63310_RAW/GSM1545538_purep53.txt",
"GSE63310_RAW/GSM1545539_JMS8-2.txt",
"GSE63310_RAW/GSM1545540_JMS8-3.txt",
"GSE63310_RAW/GSM1545541_JMS8-4.txt",
"GSE63310_RAW/GSM1545542_JMS8-5.txt",
"GSE63310_RAW/GSM1545544_JMS9-P7c.txt",
"GSE63310_RAW/GSM1545545_JMS9-P8c.txt")
# 讀取第一個文件，看下前5行
read.delim(files[1], nrow=5)
#> EntrezID GeneLength Count
#> 1 497097 3634 1
#> 2 100503874 3259 0
#> 3 100038431 1634 0
#> 4 19888 9747 0
#> 5 20671 3130 1
批量讀取+轉換:
- 方法一：一次性讀取全部，用readDGE函數(shù)
suppressMessages(library(limma))
suppressMessages(library(edgeR))
#輸入基因名和count的列
x <- readDGE(files, columns = c(1,3))
# x包含了兩部分信息：sample和count
class(x)
#> [1] "DGEList"
#> attr(,"package")
#> [1] "edgeR"
dim(x)
#> [1] 27179 9
- 方法二：分別讀取，再用DGEList函數(shù)轉換全部
轉換下樣本信息
存儲樣本與cell type(basal,LP,ML)/ genotype(wild-type, knockout)/ phenotype(status, sex, age)/
sample treatment(drug, control)/ batch effect (date, lane)
# 輸入start、stop 位置
samplenames <- substring(colnames(x), 25, nchar(colnames(x)))
samplenames
#> [1] "10_6_5_11" "9_6_5_11" "purep53" "JMS8-2" "JMS8-3" "JMS8-4"
#> [7] "JMS8-5" "JMS9-P7c" "JMS9-P8c"
#去掉列名的重復項，重新命名
colnames(x) <- samplenames
#添加分組信息，比如細胞類型
group <- as.factor(c("LP", "ML", "Basal", "Basal", "ML", "LP",
"Basal", "ML", "LP"))
x$samples$group <- group
# 添加lane信息
lane <- as.factor(rep(c("L004","L006","L008"), c(3,4,2)))
x$samples$lane <- lane
x$samples
#> files group lib.size
#> 10_6_5_11 GSE63310_RAW/GSM1545535_10_6_5_11.txt LP 32863052
#> 9_6_5_11 GSE63310_RAW/GSM1545536_9_6_5_11.txt ML 35335491
#> purep53 GSE63310_RAW/GSM1545538_purep53.txt Basal 57160817
#> JMS8-2 GSE63310_RAW/GSM1545539_JMS8-2.txt Basal 51368625
#> JMS8-3 GSE63310_RAW/GSM1545540_JMS8-3.txt ML 75795034
#> JMS8-4 GSE63310_RAW/GSM1545541_JMS8-4.txt LP 60517657
#> JMS8-5 GSE63310_RAW/GSM1545542_JMS8-5.txt Basal 55086324
#> JMS9-P7c GSE63310_RAW/GSM1545544_JMS9-P7c.txt ML 21311068
#> JMS9-P8c GSE63310_RAW/GSM1545545_JMS9-P8c.txt LP 19958838
#> norm.factors lane
#> 10_6_5_11 1 L004
#> 9_6_5_11 1 L004
#> purep53 1 L004
#> JMS8-2 1 L006
#> JMS8-3 1 L006
#> JMS8-4 1 L006
#> JMS8-5 1 L006
#> JMS9-P7c 1 L008
#> JMS9-P8c 1 L008
增加一個基因注釋信息
基因注釋涉及到ID轉換，方法一：clusterProfiler
# 數(shù)據(jù)集中使用的是EntrezID
suppressMessages(library(Mus.musculus, supr))
suppressMessages(library(clusterProfiler))
geneid <- rownames(x)
# 但是好像沒有染色體信息
genes <- bitr(geneid, OrgDb = Mus.musculus, fromType = "ENTREZID", toType = "SYMBOL" )
#> 'select()' returned 1:1 mapping between keys and columns
#> Warning in bitr(geneid, OrgDb = Mus.musculus, fromType = "ENTREZID", toType
#> = "SYMBOL"): 3.52% of input gene IDs are fail to map...
head(genes)
#> ENTREZID SYMBOL
#> 1 497097 Xkr4
#> 2 100503874 Gm19938
#> 3 100038431 Gm10568
#> 4 19888 Rp1
#> 5 20671 Sox17
#> 6 27395 Mrpl15
方法二：biomaRt（屬于Ensembl）
suppressMessages(library(Mus.musculus))
geneid <- rownames(x)
# 增加Symbol名稱和染色體信息
genes <- select(Mus.musculus, keys=geneid,
columns=c("SYMBOL", "TXCHROM"),keytype="ENTREZID")
#> 'select()' returned 1:many mapping between keys and columns
head(genes)
#> ENTREZID SYMBOL TXCHROM
#> 1 497097 Xkr4 chr1
#> 2 100503874 Gm19938 <NA>
#> 3 100038431 Gm10568 <NA>
#> 4 19888 Rp1 chr1
#> 5 20671 Sox17 chr1
#> 6 27395 Mrpl15 chr1
需要注意的是：ID轉換中可能出現(xiàn)多個同樣的Entrez ID，比如某個基因出現(xiàn)在不同染色體，
就會在不同位置出現(xiàn)同一個ID。這種情況需要先檢查下
dup <- genes$ENTREZID[duplicated(genes$ENTREZID)]
genes[genes$ENTREZID %in% dup,][1:10,]
#> ENTREZID SYMBOL TXCHROM
#> 5360 100316809 Mir1906-1 chr12
#> 5361 100316809 Mir1906-1 chrX
#> 9563 12228 Btg3 chr16
#> 9564 12228 Btg3 chr17
#> 11350 433182 Eno1b chr4
#> 11351 433182 Eno1b chr18
#> 11543 100217457 Snord58b chr14
#> 11544 100217457 Snord58b chr18
#> 13226 735274 Mir684-1 chr2
#> 13227 735274 Mir684-1 chr5
# 出現(xiàn)多次的基因只統(tǒng)計出現(xiàn)第一次的，得到的是一個位置向量，即位于多少行
mat <- match(geneid, genes$ENTREZID)
genes <- genes[mat,]
#再次驗證
genes[genes$ENTREZID %in% dup,][1:5,]
#> ENTREZID SYMBOL TXCHROM
#> 5360 100316809 Mir1906-1 chr12
#> 9563 12228 Btg3 chr16
#> 11350 433182 Eno1b chr4
#> 11543 100217457 Snord58b chr14
#> 13226 735274 Mir684-1 chr2
x$genes <- genes
數(shù)據(jù)預處理
raw-count標準化
不管是差異分析還是其他用到count值的相關分析，如WGCNA，都基本不直接用原始count值，原因是這個數(shù)值容易受到外界的干擾。一般來講，影響raw count的有兩個因素：
- 測序深度
這個比較好理解，深度越大，得到的reads數(shù)就越多，當然count差異就越大 - RNA的組織/時間特異性
不同時期基因表達量是不一樣的，很難保證兩組樣本中RNA是同一時期的
因此需要先對測序深度做一個轉換，也就是所謂的“標準化”，減少外在影響
比較流行的標準化方式有：counts per million (CPM), log2-counts per million (log-CPM), reads per kilobase of transcript per million (RPKM), fragments per kilobase of transcript per million (FPKM); 其中CPM與logCPM沒有考慮基因長度（要求略低），直接對矩陣進行轉換。
# edegR中的CPM與log-CPM
cpm <- cpm(x)
lcpm <- cpm(x, log=TRUE)
# 如果用edegR中的rpkm，需要提供gene length信息
# 用TMM方法
x <- calcNormFactors(x, method = "TMM")
去掉不表達基因
所有的表達量數(shù)據(jù)都存在表達與非表達基因，拿到表達量數(shù)據(jù)先看下非表達基因占多少。
#這里有9組樣本，因此可以統(tǒng)計9組中表達量都為0的基因
table(rowSums(x$counts==0)==9)
#>
#> FALSE TRUE
#> 22026 5153
基因必須至少在一組處理中或者任意3個樣本中有表達量才能被保留。當然，非表達的需要去除，以減少后續(xù)分析壓力
keep.exprs <- rowSums(cpm>1)>=3
# 行代表基因，列代表樣本；這里選擇全部樣本
x <- x[keep.exprs,, keep.lib.sizes=FALSE]
dim(x)
#> [1] 14165 9
基因表達分布的標準化
使用trimmed mean of M-values（TMM），也就是目前認為最準確的標準化方式
為了解決樣本間由于組織、時間特異性而引起的誤差，TMM首先從所有樣本中挑選一個樣本作為參照，當對其他樣本進行歸一化時，就只考慮參照樣本和待歸一化的樣本間共有的RNA
一般來講，mRNA數(shù)據(jù)使用TMM就可以【但是單細胞測序及全轉錄組一般是不行的，因為它們的差異表達比較多】
x <- calcNormFactors(x, method = "TMM")
非監(jiān)督式樣本聚類
在正式開始差異分析之前，可以使用非監(jiān)督式學習處理一下：
什么是非監(jiān)督？
將數(shù)據(jù)按相似度聚類（clustering）成不同的分組或者用降維的方式，保留基本的數(shù)據(jù)結構，同時壓縮數(shù)據(jù)。包含了K 均值聚類、層次聚類、主成分分析（PCA）和奇異值分解（SVD）
- K均值聚類
目標是將數(shù)據(jù)點分組，讓不同聚類中的數(shù)據(jù)點不同，同一聚類中的數(shù)據(jù)點類似；使用這種方法希望結果的數(shù)據(jù)點被聚類為 K 組。K 大分組就小，就有更多粒度；K 小時，分組就大，粒度更少（可以理解成零散程度）
比如一個領導者一般氣場比較大，周圍會有一些像被吸附了的幫手；如果這樣的領導只有一個，那么群眾們都圍著他轉；如果英雄輩出，那么各自會形成一個小中心圈，分別集結一小幫人
- 層次聚類
我們?nèi)粘Ｕf的“做個聚類樹”，就是指的這個。目標是構建一個樹狀層次，讓相似的離得更近。它是這么做的：
- 先從 第N 個聚類開始，每個數(shù)據(jù)點一個聚類；
- 將彼此靠得最近的兩個聚類融合為一個，那么現(xiàn)在有 N-1 個聚類
- 重新計算這些聚類之間的距離（方法很多，其中一種“average-linkage clustering”將兩個聚類之間的距離看作是聚類中元素所有距離的平均值）
- 重復2-3，再選擇聚類的數(shù)量（就是要將整體分成幾部分），畫條水平線
- 降維
兩種辦法：主成分分析和奇異值分解
- 主成分分析：假設我們現(xiàn)在有1w維空間的數(shù)據(jù)，然后怎么樣才能知道數(shù)據(jù)間的聯(lián)系呢，哪些比較強勢？我們可以將原始數(shù)據(jù)映射到另一個更緊湊的空間中（比如只有100維），減少了復雜度（也就是維度），還保證了原始結構不變（即方差不變）
- 奇異值分解（SVD）
如果說我們的數(shù)據(jù)是一個a X b的大型矩陣，通過SVD將大矩陣拆分成幾個小矩陣的乘積，從中找一個對角矩陣叫做“奇異值”，之所以奇異，是因為它可以用于壓縮原來的矩陣。如果我們丟棄奇異值中最小的 20% 以及其他矩陣中相關的列，我們就可以節(jié)省大量空間，并能很好地表示原來的矩陣
這里怎么做？
這里采用limma的multidimensional scaling (MDS) plot即“多維尺度變換”，這個會用非監(jiān)督式顯示樣本之間的異同，在正式差異分析前幫助我們判斷潛在的差異樣本。它的結果會將所有樣本劃分成幾個維度，第一維度的樣本代表了最大的差異，維度越往后，相似度就越大，差異分析的結果越不理想
當實驗設計中存在多個實驗因子時，每個實驗因子都要進行這幾個維度的檢測。這樣做，還能看出哪個實驗因子對差異表達的貢獻最大，基本沒什么影響的因子，就無法參與后續(xù)分析。
suppressMessages(library(RColorBrewer))
# 這里就用最簡單的CPM值
lcpm <- cpm(x, log=TRUE)
par(mfrow=c(1,2))
col.group <- group
levels(col.group) <- brewer.pal(nlevels(col.group), "Set1")
col.group <- as.character(col.group)
col.lane <- lane
levels(col.lane) <- brewer.pal(nlevels(col.lane), "Set2")
col.lane <- as.character(col.lane)
# 1、2維度
plotMDS(lcpm, labels=group, col=col.group)
title(main="A. Sample groups")
# 3、4維度
plotMDS(lcpm, labels=lane, col=col.lane, dim=c(3,4))
title(main="B. Sequencing lanes")
歡迎關注我們的公眾號～_～
我們是兩個農(nóng)轉生信的小碩，打造生信星球，想讓它成為一個不拽術語、通俗易懂的生信知識平臺。需要幫助或提出意見請后臺留言或發(fā)送郵件到Bioplanet520@outlook.com
