為了完成學(xué)徒任務(wù)：HCC珍促，CHC铃辖，CC這3組，跟healthy的分開(kāi)比較猪叙，然后3個(gè)火山圖
具體題目參考http://www.reibang.com/p/d521459ae1d0
操作方法參見(jiàn)http://www.reibang.com/p/f635f60e1c8a

1. 初始配置

rm(list = ls())
options(stringsAsFactors = F)
library(limma)
library(Glimma)
library(edgeR)
BiocManager::install("Homo.sapiens")
library(Homo.sapiens)
library(biomaRt)#biomaRt主要使用Ensembl基因ID進(jìn)行檢索

2. 數(shù)據(jù)整合

2.1 讀入計(jì)數(shù)數(shù)據(jù)

#url <- "https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE84073&format=file"
#utils::download.file(url, destfile="GSE84073_RAW.tar", mode="wb") 
#如果下載全部數(shù)據(jù)是134.9 MB娇斩，下載時(shí)間有點(diǎn)長(zhǎng)，估計(jì)可能因?yàn)榫W(wǎng)絡(luò)原因出錯(cuò)穴翩，因此可以手動(dòng)下載需要分析的樣本并放在當(dāng)前文件夾
#手動(dòng)下載的是一個(gè)2.8兆的文件
utils::untar("GSE84073_RAW.tar", exdir = ".")#解壓縮文件
directory <- "."#設(shè)置需要操作的目錄
sampleFiles <- grep("Counts",list.files(directory),value=TRUE)#提取文件名稱
sampleNames <- sub("(.*)(.gz)", "\\1", sampleFiles)
files <- sampleNames
for(i in paste(files, ".gz", sep=""))
  R.utils::gunzip(i, overwrite=TRUE)
files[1]#查看第一個(gè)文件名
read.delim(files[1], nrow=5)#查看第一個(gè)文件內(nèi)容
x <- readDGE(files, columns=c(1,3))#將所有文件讀入對(duì)象DEGLlist
class(x)
dim(x)

#當(dāng)樣本量少的時(shí)候可以通過(guò)復(fù)制來(lái)賦值犬第，多的情況下就還是用老師推薦的函數(shù)了
#files <- c("GSM2653819_Counts_notmergedTR_Healthy1_Tissue_1.txt", "GSM2653820_Counts_notmergedTR_Healthy1_Tissue_2.txt", ...)

2.2 組織樣品信息

• 我們的分析僅包含了此實(shí)驗(yàn)中的Healthy、HCC芒帕、CHC和CC樣品歉嗓，從sampleFiles可以看到，這些樣品來(lái)自人的肝臟組織背蟆，可以提取不同的疾病表型鉴分，由于是不同時(shí)期取樣，可以提取批次信息淆储。

samplenames <- substring(colnames(x),31, nchar(colnames(x)))
samplenames
colnames(x) <- samplenames
group <- as.factor(rep(c('Healthy','HCC','CHC','CC'),c(4,3,3,4)))
group
x$samples$group <- group
tlist=rep(c("Tissue_1","Tissue_2","Tissue","Tissue"),4)
length(tlist)
tlist=tlist[-c(8,12)]
tlist
lane <- as.factor(tlist)
x$samples$lane <- lane
x$samples

2.3 組織基因注釋

• DGEList對(duì)象中的第二個(gè)數(shù)據(jù)框名為genes冠场，用于存儲(chǔ)與計(jì)數(shù)矩陣的行相關(guān)聯(lián)的基因水平的信息家浇。 為檢索這些信息本砰，我們可以使用包含特定物種信息的包。
• 人類的Homo.sapiens (Bioconductor Core Team 2016a)钢悲；或者也可以使用biomaRt 包 (Durinck et al. 2005, 2009)点额，它通過(guò)接入Ensembl genome數(shù)據(jù)庫(kù)來(lái)進(jìn)行基因注釋舔株。
• 可以檢索的信息類型包括基因符號(hào)（gene symbols）、基因名稱（gene names）还棱、染色體名稱和位置（chromosome names and locations）载慈、Entrez基因ID（Entrez gene IDs）、Refseq基因ID（Refseq gene IDs）和Ensembl基因ID（Ensembl gene IDs）等珍手。biomaRt主要使用Ensembl基因ID進(jìn)行檢索办铡。
• 我可以比較使用Homo.sapiens包和biomaRt 包，利用此次數(shù)據(jù)集中的Entrez基因ID來(lái)檢索相關(guān)的基因符號(hào)和染色體信息琳要。
  • 先了解一下Homo.sapiens包

  Based on genome: hg19
  The OrgDb gene id ENTREZID is mapped to the TxDb gene id GENEID .

#Examples
Homo.sapiens
cls <- columns(Homo.sapiens)
cls
cls <- cls[c(1,19,45)]
kts <- keytypes(Homo.sapiens)
kt <- kts[2]
kts
ks <- head(keys(Homo.sapiens, keytype=kts[2]))
ks
res <- select(Homo.sapiens, keys=ks, columns=cls, keytype=kt)
head(res)

經(jīng)過(guò)嘗試后發(fā)現(xiàn)基本跟小鼠的那個(gè)代碼差不多

geneid <- rownames(x)
genes <- select(Homo.sapiens, keys=geneid, columns=c("SYMBOL", "TXCHROM"), 
                keytype="ENSEMBL" )
head(genes)
# #試一下biomart
# ensembl = useMart( "ensembl", dataset = "hsapiens_gene_ensembl" )
# #str(ensembl)
# genemap <- getBM(listAttributes = c("ensembl_gene_id", "entrezgene", "hgnc_symbol"),
#                   filters = "ensembl_gene_id",
#                   values = rownames(x),
#                   mart = ensembl )
# genemap <- genemap[match(rownames(x),genemap$ensembl_gene_id),]
# 
# rownames(x) <-genemap$hgnc_symbol
# x <- x[!is.na(rownames(x)),]
# #39597 genes
# x <- x[rownames(x) != "", ]
# #27884 genes
# x <- x[!duplicated(rownames(x)),]

也嘗試了 OvCaBiobank/rnaseq/allPatients_deseq.R 沒(méi)有成功
估計(jì)是之前的代碼需要更新了寡具。由于我得到的矩陣的基因名稱是ENSG開(kāi)頭的，所以只需要把keytype="ENSEMBL"做一下更改即可稚补，雖然只是這一個(gè)小地方童叠，花了我整整1個(gè)小時(shí)才搞明白，就是因?yàn)椴欢@個(gè)包的內(nèi)容课幕。

• 與任何基因ID一樣厦坛，Entrez基因ID可能不能一對(duì)一地匹配我們想獲得的基因信息。在處理之前乍惊，檢查重復(fù)的基因ID和弄清楚重復(fù)的來(lái)源非常重要杜秸。
• • 為了處理重復(fù)的基因ID，我們可以合并來(lái)自多重匹配基因的所有染色體信息润绎，比如將基因Gm1987分配到chr4 and chr4_JH584294_random亩歹，或選取其中一條染色體來(lái)代表具有重復(fù)注釋的基因。為了簡(jiǎn)單起見(jiàn)凡橱，我們選擇后者小作，保留每個(gè)基因ID第一次出現(xiàn)的信息。

genes <- genes[!duplicated(genes$ENSEMBL),]

• 在此例子中稼钩，注釋與數(shù)據(jù)對(duì)象中的基因順序是相同的顾稀。
• 如果由于缺失和／或重新排列基因ID導(dǎo)致其順序不一致，可以用match來(lái)正確排序基因坝撑。
• 然后將基因注釋的數(shù)據(jù)框加入數(shù)據(jù)對(duì)象静秆，數(shù)據(jù)即被整潔地整理入一個(gè)DGEList對(duì)象中，它包含原始計(jì)數(shù)數(shù)據(jù)和相關(guān)的樣品信息和基因注釋巡李。

x$genes <- genes
x

3. 數(shù)據(jù)預(yù)處理

3.1 原始數(shù)據(jù)尺度轉(zhuǎn)換

• 在此處抚笔，使用edgeR中的cpm函數(shù)將原始計(jì)數(shù)轉(zhuǎn)換為CPM和log-CPM值。

cpm <- cpm(x)
lcpm <- cpm(x, log=TRUE, prior.count=2)

3.2 刪除低表達(dá)基因

• 所有數(shù)據(jù)集中都混有表達(dá)的基因與不表達(dá)的基因侨拦。盡管我們想要檢測(cè)在一種條件中表達(dá)但再另一種條件中不表達(dá)的基因殊橙，也有一些基因在所有樣品中都不表達(dá)。

table(rowSums(x$counts==0)==14)

• edgeR包中的filterByExpr函數(shù)提供了自動(dòng)過(guò)濾基因的方法，可保留盡可能多的有足夠表達(dá)計(jì)數(shù)的基因膨蛮。

keep.exprs <- filterByExpr(x, group=group)
x <- x[keep.exprs,, keep.lib.sizes=FALSE]
dim(x)

• 使用這個(gè)標(biāo)準(zhǔn)叠纹，基因的數(shù)量減少到了16088個(gè)，約為開(kāi)始時(shí)數(shù)量的69%敞葛。過(guò)濾后的log-CPM值顯示出每個(gè)樣本的分布基本相同（下圖B部分）誉察。需要注意的是，從整個(gè)DGEList對(duì)象中取子集時(shí)同時(shí)刪除了被過(guò)濾的基因的計(jì)數(shù)和其相關(guān)的基因信息惹谐。過(guò)濾后的DGEList對(duì)象為留下的基因保留了相對(duì)應(yīng)的基因信息和計(jì)數(shù)持偏。
• 下方給出的是繪圖所用代碼。參考之前的學(xué)習(xí)記錄

lcpm.cutoff <- log2(10/M + 2/L)
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.26), las=2, main="", xlab="")
title(main="A. Raw data", xlab="Log-cpm")
abline(v=lcpm.cutoff, 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.26), las=2, main="", xlab="")
title(main="B. Filtered data", xlab="Log-cpm")
abline(v=lcpm.cutoff, 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")

figure1

3.3 歸一化基因表達(dá)分布

• 在樣品制備或測(cè)序過(guò)程中氨肌，不具備生物學(xué)意義的外部因素會(huì)影響單個(gè)樣品的表達(dá)综液。比如說(shuō)，在實(shí)驗(yàn)中第一批制備的樣品會(huì)總體上表達(dá)高于第二批制備的樣品儒飒。假設(shè)所有樣品表達(dá)值的范圍和分布都應(yīng)當(dāng)相似谬莹，需要進(jìn)行歸一化來(lái)確保整個(gè)實(shí)驗(yàn)中每個(gè)樣本的表達(dá)分布都相似。
• 盡管如此桩了，我們依然需要使用edgeR中的calcNormFactors函數(shù)附帽，用TMM(Robinson and Oshlack 2010)方法進(jìn)行歸一化。此處計(jì)算得到的歸一化系數(shù)被用作文庫(kù)大小的縮放系數(shù)井誉。當(dāng)我們使用DGEList對(duì)象時(shí)蕉扮，這些歸一化系數(shù)被自動(dòng)存儲(chǔ)在x$samples$norm.factors。對(duì)此數(shù)據(jù)集而言颗圣，TMM歸一化的作用比較溫和喳钟，這體現(xiàn)在所有的縮放因子都相對(duì)接近1。

x <- calcNormFactors(x, method = "TMM")
x$samples$norm.factors

• 為了更好地可視化表現(xiàn)出歸一化的影響在岂，我們復(fù)制了數(shù)據(jù)并進(jìn)行了調(diào)整奔则，使得第一個(gè)樣品的計(jì)數(shù)減少到了其原始值的5%，而第二個(gè)樣品增大到了5倍蔽午，第11易茬，12增大了2倍。
• 歸一化之后就比較整齊了

x2 <- x
x2$samples$norm.factors <- 1
x2$counts[,1] <- ceiling(x2$counts[,1]*0.05)
x2$counts[,2] <- x2$counts[,2]*5
x2$counts[,11] <- x2$counts[,2]*2
x2$counts[,12] <- x2$counts[,2]*2
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
lcpm <- cpm(x2, log=TRUE)
boxplot(lcpm, las=2, col=col, main="")
title(main="B. Example: Normalised data",ylab="Log-cpm")

figure 2

3.4 對(duì)樣本的無(wú)監(jiān)督聚類

• 在我們看來(lái)及老，用于檢查基因表達(dá)分析的最重要的探索性圖表之一便是MDS圖或其余類似的圖抽莱。這種圖表使用無(wú)監(jiān)督聚類方法展示出了樣品間的相似性和不相似性，能讓我們?cè)谶M(jìn)行正式的檢驗(yàn)之前對(duì)于能檢測(cè)到多少差異表達(dá)基因有個(gè)大致概念骄恶。理想情況下食铐，樣本會(huì)在不同的實(shí)驗(yàn)組內(nèi)很好的聚類，且可以鑒別出遠(yuǎn)離所屬組的樣本僧鲁，并追蹤誤差或額外方差的來(lái)源虐呻。如果存在技術(shù)重復(fù)象泵，它們應(yīng)當(dāng)互相非常接近。

• 這樣的圖可以用limma中的plotMDS函數(shù)繪制铃慷。第一個(gè)維度表示能夠最好地分離樣品且解釋最大比例的方差的引導(dǎo)性的倍數(shù)變化（leading-fold-change）单芜，而后續(xù)的維度的影響更小蜕该，并與之前的維度正交犁柜。當(dāng)實(shí)驗(yàn)設(shè)計(jì)涉及到多個(gè)因子時(shí)，建議在多個(gè)維度上檢查每個(gè)因子堂淡。如果在其中一些維度上樣本可按照某因子聚類馋缅，這說(shuō)明該因子對(duì)于表達(dá)差異有影響，在線性模型中應(yīng)當(dāng)將其包括進(jìn)去绢淀。反之萤悴，沒(méi)有或者僅有微小影響的因子在下游分析時(shí)應(yīng)當(dāng)被剔除。

• 在這個(gè)數(shù)據(jù)集中皆的，可以看出樣本在維度1和2能很好地按照實(shí)驗(yàn)分組聚類覆履，隨后在維度3按照測(cè)序道（樣品批次）分離（如下圖所示）。請(qǐng)記住费薄，第一維度解釋了數(shù)據(jù)中最大比例的方差硝全，需要注意到，當(dāng)我們向高維度移動(dòng)楞抡，維度上的取值范圍會(huì)變小伟众。

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)
plotMDS(lcpm, labels=group, col=col.group)
title(main="A. subtype groups")
plotMDS(lcpm, labels=lane, col=col.lane, dim=c(3,4))
title(main="B.tissue times")

figure 3

• 從上圖可見(jiàn)，樣本的聚類不是太好召廷，先看看結(jié)果吧

4. 差異表達(dá)分析

4.1 創(chuàng)建設(shè)計(jì)矩陣和對(duì)比

• 在此研究中凳厢，我想知道哪些基因在不同亞型之間以不同水平表達(dá)。在我們的分析中竞慢，假設(shè)基礎(chǔ)數(shù)據(jù)是正態(tài)分布的先紫，為其擬合一個(gè)線性模型。在此之前筹煮，需要?jiǎng)?chuàng)建一個(gè)包含疾病分類亞型以及樣本提扰莺ⅰ（批次）信息的設(shè)計(jì)矩陣。

design <- model.matrix(~0+group+lane)
colnames(design) <- gsub("group", "", colnames(design))
design

先模仿看看：

contr.matrix <- makeContrasts(
  HealthyvsCC = Healthy - CC, 
  HealthyvsCHC = Healthy - CHC, 
  HealthyvsHCC = Healthy - HCC, 
  levels = colnames(design))
contr.matrix

par(mfrow=c(1,2))
v <- voom(x, design, plot=TRUE)
v
vfit <- lmFit(v, design)
vfit <- contrasts.fit(vfit, contrasts=contr.matrix)
efit <- eBayes(vfit)
plotSA(efit, main="Final model: Mean-variance trend")

summary(decideTests(efit))

tfit <- treat(vfit, lfc=1)
dt <- decideTests(tfit)
summary(dt)

de.common <- which(dt[,1]!=0 & dt[,2]!=0)
length(de.common)
head(tfit$genes$SYMBOL[de.common], n=20)
vennDiagram(dt[,1:2], circle.col=c("turquoise", "salmon"))
write.fit(tfit, dt, file="results.txt")

Healthy.vs.CC <- topTreat(tfit, coef=1, n=Inf)
Healthy.vs.CHC <- topTreat(tfit, coef=2, n=Inf)
head(Healthy.vs.CC)
head(Healthy.vs.CHC)

plotMD(tfit, column=1, status=dt[,1], main=colnames(tfit)[1], 
       xlim=c(-8,13))

figure 4

這個(gè)探索的過(guò)程中得到如下結(jié)果

> summary(decideTests(efit))
       HealthyvsCC HealthyvsCHC HealthyvsHCC
Down           849          408            4
NotSig       13948        15409        16084
Up            1291          271            0
> tfit <- treat(vfit, lfc=1)
> dt <- decideTests(tfit)
> summary(dt)
       HealthyvsCC HealthyvsCHC HealthyvsHCC
Down           176           64            0
NotSig       15432        15985        16088
Up             480           39            0

關(guān)于火山圖寺谤，明天再學(xué)習(xí)進(jìn)一步美化吧仑鸥，好困……