劉小澤寫于18.10.26
今晚10.15開始寫，先寫這么多吧轻抱，明天繼續(xù)
第一部分 背景介紹
選用GEO數(shù)據(jù)：GSE63310
從雌性小鼠乳腺中提取了三種類型的細(xì)胞：basal旦部、luminal progenitor（LP-乳腺癌早期前體細(xì)胞）士八、mature luminal （ML），各三個(gè)重復(fù)蘸秘。RNA利用Hiseq 2000+100bp單端測序蝗茁；參考基因組選擇mm10，Rsubread比對+featureCounts定量
轉(zhuǎn)錄組表達(dá)數(shù)據(jù)地址：
http://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE63310&format=file
第二部分 數(shù)據(jù)準(zhǔn)備
讀取表達(dá)量數(shù)據(jù) =>得到DEGList x
cd ~/Download
wget -c http://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE63310&format=file
tar xvf GSE63310_RAW.tar
unzip *
### 設(shè)置好rstudio 的project
files <- c("GSM1545535_10_6_5_11.txt", "GSM1545536_9_6_5_11.txt",
"GSM1545538_purep53.txt","GSM1545539_JMS8-2.txt",
"GSM1545540_JMS8-3.txt","GSM1545541_JMS8-4.txt",
"GSM1545542_JMS8-5.txt","GSM1545544_JMS9-P7c.txt",
"GSM1545545_JMS9-P8c.txt")
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
#接下來利用edgeR的readDEG將所有表達(dá)矩陣讀入并且組合成一個(gè)大的矩陣
library(limma)
library(edgeR)
x <- readDGE(files, columns=c(1,3)) # 只要它的Entrez ID 和count值
dim(x) # x中有表達(dá)量和樣本信息，樣本信息又包括了group粘舟、libsize、lane等信息
#[1] 27179 9
# 如果一開始所有的表達(dá)量都存在一個(gè)矩陣中霞揉，那么直接用DEGList()轉(zhuǎn)換就好
整理樣本分組信息 =》 對DEGList x的列進(jìn)行操作
下游分析之前晰骑，需要將整合的矩陣的列名與實(shí)驗(yàn)樣本設(shè)計(jì)關(guān)聯(lián)起來（包括不同的分組信息以及不同組的重復(fù)信息）
比如：細(xì)胞的basal、LP秽荞、ML分組抚官；基因型的wild-type凌节、knock-out分組；表型的disease朴上、status卒煞、sex、age分組衣撬、樣本處理的drug淮韭、control分組贴届、關(guān)于批次Batch的不同樣本采集蜡吧、測序日期昔善、測序lane等
這里的DEGList中包含的樣本信息有：細(xì)胞類型、批次信息（測序的lane）翩概，每種信息都是因子型factor變量，但是各自因子的level不同
# 先簡化GEO的ID名稱牍鞠，也就是把GSM去掉
samplenames <- substring(colnames(x), 12, nchar(colnames(x)))
samplenames
##[1] "10_6_5_11" "9_6_5_11" "purep53"
##[4] "JMS8-2" "JMS8-3" "JMS8-4"
##[7] "JMS8-5" "JMS9-P7c" "JMS9-P8c"
colnames(x) <- samplenames
#然后我們可以自己設(shè)計(jì)分組信息
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 norm.factors lane
##10_6_5_11 GSM1545535_10_6_5_11.txt LP 32863052 1 L004
##9_6_5_11 GSM1545536_9_6_5_11.txt ML 35335491 1 L004
##purep53 GSM1545538_purep53.txt Basal 57160817 1 L004
##JMS8-2 GSM1545539_JMS8-2.txt Basal 51368625 1 L006
##JMS8-3 GSM1545540_JMS8-3.txt ML 75795034 1 L006
##JMS8-4 GSM1545541_JMS8-4.txt LP 60517657 1 L006
##JMS8-5 GSM1545542_JMS8-5.txt Basal 55086324 1 L006
##JMS9-P7c GSM1545544_JMS9-P7c.txt ML 21311068 1 L008
##JMS9-P8c GSM1545545_JMS9-P8c.txt LP 19958838 1 L008
整理基因注釋信息=》 對DEGList x的行進(jìn)行操作
source("https://bioconductor.org/biocLite.R")
options(BioC_mirror="http://mirrors.ustc.edu.cn/bioc/")
biocLite("Mus.musculus")
library(Mus.musculus)
geneid <- rownames(x)
genes <- select(Mus.musculus, keys=geneid, columns=c("SYMBOL", "TXCHROM"),
keytype="ENTREZID") #轉(zhuǎn)換基因id难述，用clusterProfiler的bitr函數(shù)也可以;另外這個(gè)還增加了染色體信息
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
需要注意的是：Entrez ID可能并不是和基因信息一一匹配的胁后，可能同樣的ID會(huì)匹配到不同染色體嗦枢，因此需要檢查有沒有重復(fù)出現(xiàn)的Entrez ID，以保證注釋和我們的DEGList之間的基因順序是一致的
dup <- genes$ENTREZID[duplicated(genes$ENTREZID)]
genes[genes$ENTREZID %in% dup,][1:5,]
# 果然發(fā)現(xiàn)了不同染色體上重復(fù)的ID
## 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
# 然后把重復(fù)的基因挑出來【重復(fù)的只統(tǒng)計(jì)一次】
mat <- match(geneid, genes$ENTREZID)
genes <- genes[mat,]
genes[genes$ENTREZID %in% dup,][1:5,]
這時(shí)在看x這個(gè)DEGList敲才，就是一個(gè)標(biāo)準(zhǔn)的包含原始count數(shù)據(jù)紧武、樣本信息的表達(dá)量矩陣
第三部分 數(shù)據(jù)預(yù)處理
原始數(shù)據(jù)的轉(zhuǎn)換
進(jìn)行差異表達(dá)一般都不會(huì)用raw counts的敏储，因?yàn)榇嬖跍y序深度已添、文庫大小的差別，這樣的結(jié)果是不準(zhǔn)確的
一般的做法是：利用標(biāo)準(zhǔn)化算法畦幢，如CPM（counts per million), log-CPM (log2-counts per million), RPKM (reads per kilobase of transcript per million), FPKM(fragments per kilobase of
transcript per million)等去除文庫大小缆蝉、深度的影響。和RPKM黍瞧、FPKM不同的是原杂，CPM和log-CPM不需要考慮feature length的差異穿肄，也就是說基因長度在統(tǒng)計(jì)時(shí)被當(dāng)成常數(shù)际看，只考慮不同處理下的不同矢否，而不會(huì)受長度的影響
cpm使用cpm()函數(shù)兴喂；RPKM使用rpkm函數(shù)，都屬于edegR
cpm <- cpm(x)
lcpm <- cpm(x, log=TRUE)
去掉不感興趣的基因
所有的數(shù)據(jù)集中都會(huì)存在表達(dá)的和不表達(dá)的基因畏鼓，我們感興趣的是在一個(gè)條件下表達(dá)壶谒，另一個(gè)條件不表達(dá)的汗菜。
看下有多少基因在所有樣本中表達(dá)量都為0
table(rowSums(x$counts==0)==9)
##
## FALSE TRUE
## 22026 5153
表達(dá)量都為0的占比達(dá)到了19%
過濾基因=>標(biāo)準(zhǔn)就是cpm至少一組或整個(gè)實(shí)驗(yàn)中有三個(gè)樣本大于1
keep.exprs <- rowSums(cpm>1)>=3
x <- x[keep.exprs,, keep.lib.sizes=FALSE]
dim(x)
作圖檢查：虛線就是cpm為1（即log cpm為0的時(shí)候）作為判斷閾值，看超過這個(gè)線的有幾個(gè)
library(RColorBrewer)
nsamples <- ncol(x)
col <- brewer.pal(nsamples, "Paired")
par(mfrow=c(1,2))
plot(density(lcpm[,1]), col=col[1], lwd=2, ylim=c(0,0.21), las=2,
main="", xlab="")
title(main="A. Raw data", xlab="Log-cpm")
abline(v=0, lty=3)
for (i in 2:nsamples){
den <- density(lcpm[,i])
lines(den$x, den$y, col=col[i], lwd=2)
}
legend("topright", samplenames, text.col=col, bty="n")
lcpm <- cpm(x, log=TRUE)
plot(density(lcpm[,1]), col=col[1], lwd=2, ylim=c(0,0.21), las=2,
main="", xlab="")
title(main="B. Filtered data", xlab="Log-cpm")
abline(v=0, lty=3)
for (i in 2:nsamples){
den <- density(lcpm[,i])
lines(den$x, den$y, col=col[i], lwd=2)
}
legend("topright", samplenames, text.col=col, bty="n")
基因表達(dá)分布標(biāo)準(zhǔn)化
在準(zhǔn)備試驗(yàn)樣品或者測序的過程中陨界，外界因素非逞沧幔可能會(huì)引入誤差，影響樣本的基因表達(dá)水平菌瘪。比如：第一批測序的樣本可能比第二批測序的深度要深腮敌。標(biāo)準(zhǔn)化就是為了讓每個(gè)樣本的表達(dá)量分布在整個(gè)實(shí)驗(yàn)中是相似的
如何判斷是否標(biāo)準(zhǔn)了呢？
可以通過密度圖density或者箱線圖boxplot俏扩，比如上圖就是密度分布糜工，其中B圖的log-CPM分布就比較一致，并且都在閾值的右側(cè)
如何標(biāo)準(zhǔn)化录淡？
使用TMM（trimmed mean of M-values）算法捌木，利用edgeR中函數(shù)calNormFactors()
標(biāo)準(zhǔn)化用到的normalisation factors 就在DEGList中嫉戚，x$samples$norm.factors.調(diào)取
x <- calcNormFactors(x, method = "TMM")
x$samples$norm.factors
## [1] 0.896 1.035 1.044 1.041 1.032 0.922 0.984 1.083 0.979
畫一個(gè)箱線圖可以看到標(biāo)準(zhǔn)化前后差異【模擬數(shù)據(jù)】
#模擬一個(gè)數(shù)據(jù)
x2 <- x
x2$samples$norm.factors <- 1
x2$counts[,1] <- ceiling(x2$counts[,1]*0.05) #第一個(gè)樣本count縮小到原來5%
x2$counts[,2] <- x2$counts[,2]*5 # 第二個(gè)樣本擴(kuò)大到原來5倍
#畫圖
par(mfrow=c(1,2))
lcpm <- cpm(x2, log=TRUE)
boxplot(lcpm, las=2, col=col, main="")
title(main="A. Example: Unnormalised data",ylab="Log-cpm")
x2 <- calcNormFactors(x2)
x2$samples$norm.factors
## [1] 0.0547 6.1306 1.2293 1.1705 1.2149 1.0562 1.1459 1.2613 1.1170
lcpm <- cpm(x2, log=TRUE)
boxplot(lcpm, las=2, col=col, main="")
title(main="B. Example: Normalised data",ylab="Log-cpm")
歡迎關(guān)注我們的公眾號(hào)～_～
我們是兩個(gè)農(nóng)轉(zhuǎn)生信的小碩刨裆，打造生信星球，想讓它成為一個(gè)不拽術(shù)語彼水、通俗易懂的生信知識(shí)平臺(tái)崔拥。需要幫助或提出意見請后臺(tái)留言或發(fā)送郵件到Bioplanet520@outlook.com
