生信-記一次NCBI-R語(yǔ)言-淋巴癌突變與未突變基因的差異分析

關(guān)鍵詞：基因芯片、R、篩選、預(yù)處理、差異分析

NCBI-淋巴癌突變與未突變基因的差異分析
PS：好久沒分享生信了，這是一年前做的一次生信task（準(zhǔn)確來說是2018年11月了），這里分享一下給大家，有助于一些小伙伴們想通過常規(guī)的，使用NCBI科研數(shù)據(jù)庫(kù)+R編程語(yǔ)言方式，進(jìn)行對(duì)某種癌的差異分析。最近用心做了一些更棒的生信task，相信不久會(huì)分享出來~
PS2：如果這篇筆記有什么不足，或者疑惑不解的地方，可以提出來，看到會(huì)回答（留個(gè)TODO給大家）

如果感覺對(duì)你有幫助，可以關(guān)注：專欄-生物信息學(xué)-小白成長(zhǎng)記

簡(jiǎn)單介紹

通過使用計(jì)算機(jī)對(duì)NCBI下載的基因數(shù)據(jù)，進(jìn)行篩選，其中包括質(zhì)量控制、質(zhì)量分析，聚類分析，對(duì)數(shù)據(jù)進(jìn)行預(yù)處理，并使用差異分析，獲取突變與未突變基因的差異表達(dá)基因，進(jìn)行注釋，最后生成研究結(jié)果。
本次研究主要以R作為數(shù)據(jù)處理的語(yǔ)言，多次對(duì)38組細(xì)胞進(jìn)行篩選、其中編寫了5個(gè)版本的實(shí)驗(yàn)?zāi)_本，研究結(jié)果主要有：基因的差異分析，需要對(duì)海量的數(shù)據(jù)進(jìn)行十分嚴(yán)謹(jǐn)?shù)暮Y選與預(yù)處理，突變與未突變基因在組間差異較為明顯。

一、研究對(duì)象

淋巴癌突變與未突變基因的差異分析

二、研究工具

1. NCBI

NCBI, The National Center for Biotechnology Information biomedical genomic，美國(guó)國(guó)家生物技術(shù)旨在通過提供生物醫(yī)學(xué)和基因組信息提供給熱愛科學(xué)與健康事業(yè)的人們。NCBI功能強(qiáng)大，提供了海量豐富的基因數(shù)據(jù)庫(kù)、計(jì)算機(jī)分析工具、文獻(xiàn)等，研究者可以根據(jù)自己的研究方向，選擇自己所需要的數(shù)據(jù)與處理方式，推動(dòng)我們對(duì)基因遺傳的理解，從而推進(jìn)健康和疾病的理解。

2. R

R語(yǔ)言是用于統(tǒng)計(jì)分析，圖形表示和報(bào)告的編程語(yǔ)言和軟件環(huán)境。 R語(yǔ)言由Ross Ihaka和Robert Gentleman在新西蘭奧克蘭大學(xué)創(chuàng)建，目前由R語(yǔ)言開發(fā)核心團(tuán)隊(duì)開發(fā)。 R語(yǔ)言在GNU通用公共許可證下免費(fèi)提供，并為各種操作系統(tǒng)（如Linux，Windows和Mac）提供預(yù)編譯的二進(jìn)制版本。 這種編程語(yǔ)言被命名為R語(yǔ)言，基于兩個(gè)R語(yǔ)言作者的名字的第一個(gè)字母（Robert Gentleman和Ross Ihaka），并且部分是貝爾實(shí)驗(yàn)室語(yǔ)言S的名稱。

三、研究流程

1、 數(shù)據(jù)源載入

1） 針對(duì)參考的文獻(xiàn)，通過基因芯片研究突變/未突變的套細(xì)胞淋巴瘤患者的關(guān)系
2） 對(duì)應(yīng)的數(shù)據(jù)下載源：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE36000
3） 數(shù)據(jù)下載源：

數(shù)據(jù)源下載

數(shù)據(jù)源下載到本地

4） 參考數(shù)據(jù)集：


參考數(shù)據(jù)集

2、 使用R讀取下載的數(shù)據(jù)

在使用R進(jìn)行前期準(zhǔn)備的過程中，需要先將應(yīng)用到的一些常用的lib包加載到本機(jī)（如果看不懂，那可以先跳過，繼續(xù)往下看），不同的庫(kù)會(huì)依賴到其他庫(kù)，那么也需要提前下載好。
并為了增加可讀性，我們需要將基因芯片重命名，并根據(jù)已有原本數(shù)據(jù)的分類進(jìn)行簡(jiǎn)單的分類與命名（例如，本次研究的是未突變與突變淋巴細(xì)胞）

PS：下載package的時(shí)候，更新了Bioconductor version 3.8，更新到了R的Version 3.5.1（2018-07-02），所以需要使用到BiocManager::install的方式，bioCLite已經(jīng)提示廢棄。

具體流程：

  1. 設(shè)置工作目錄
  2. 先下載并加載必要的package，已安裝并加載的，可以不用重復(fù)這些步驟
  3. 可以通過GEOquery直接讀取基因芯片的編號(hào)（例如本次差異分析中使用到的是“GSE36000”）；也可以直接在官網(wǎng)下載raw包，并通過ReadAffy的方式，這里會(huì)有AffyBatch特定的格式（本差異分析采取的是這種方式）。
  4. 查看數(shù)據(jù)類型、基本信息，并重命名芯片名稱
    （溫馨提示：R代碼部分后面會(huì)提供）

代碼參考：

# 本文字符格式基于UTF-8編碼
# 【廢棄的前數(shù)據(jù)-不合適】：1.下載數(shù)據(jù)：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE47811
# 1.下載數(shù)據(jù)：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE36000
# 2.讀取下載的數(shù)據(jù)
# 2.1 先下載并加載必要的package，已安裝并加載的，可以不用重復(fù)這些步驟
# Bioconductor用于生物芯片數(shù)據(jù)分析和基因組數(shù)據(jù)分析的軟件包
source("https://bioconductor.org/biocLite.R")
# 2.1.1 Bioconductor version 3.7版本以前，使用bioCLite下載
biocLite("affy")
biocLite("simpleaffy")
biocLite("pheatmap")
biocLite("RColorBrewer")
biocLite("affyPLM")
biocLite("GOstats")

# 2.1.2 Bioconductor version 3.8 其中biocLite已在3.8廢棄，改用BiocManager::install
BiocManager::install("affy")
BiocManager::install("simpleaffy")
BiocManager::install("pheatmap")
BiocManager::install("RColorBrewer")
BiocManager::install("affyPLM")
BiocManager::install("GOstats")
BiocManager::install("survival")

# 2.2 設(shè)置工作目錄——放置源文件
setwd("C:/libtobrain/R/R_Package/GSE36000_RAW")

# 2.3 將GSE36000_RAW.tar解壓到GSE36000_RAW目錄下
# 2.4 讀取基因數(shù)據(jù)
dir.files <- list.files(path = "C:/libtobrain/R/R_Package/GSE36000_RAW/", pattern = "*.CEL.gz")
dir.files
library(parallel)
library(BiocGenerics)
library(Biobase)
library(affy)
affy.data <- ReadAffy(filenames = dir.files)

# 2.4.1 查看數(shù)據(jù)類型
data.class(affy.data)
class(affy.data)
# 2.4.2 查看芯片的基本信息
show(affy.data)
# 2.4.3 查看芯片的樣本名
sampleNames(affy.data)

# 2.4.4 重命名樣品名
# GSM879118 MCL_IGHV_UNMUT（R1394）       mcl.ighv.unmut.1
# GSM879119 MCL_IGHV_UNMUT（R1328）       mcl.ighv.unmut.2
# GSM879120 MCL_IGHV_UNMUT（R1329）       mcl.ighv.unmut.3
# GSM879121 MCL_IGHV_UNMUT（R1306）       mcl.ighv.unmut.4
# GSM879122 MCL_IGHV_UNMUT（R1332）       mcl.ighv.unmut.5
# GSM879123 MCL_IGHV_MUT（R1399_M）       mcl.ighv.mut.1
# GSM879124 MCL_IGHV_MUT（R1400）         mcl.ighv.mut.2
# GSM879125 MCL_IGHV_UNMUT（R1587）       mcl.ighv.unmut.6
# GSM879126 MCL_IGHV_UNMUT（R1680）       mcl.ighv.unmut.7
# GSM879127 MCL_IGHV_UNMUT（268-01-5TR）  mcl.ighv.unmut.8
# GSM879128 MCL_IGHV_MUT（269-01-5TR）    mcl.ighv.mut.3
# GSM879129 MCL_IGHV_UNMUT（043-01-4TR）  mcl.ighv.unmut.9
# GSM879130 MCL_IGHV_MUT（R1302）         mcl.ighv.mut.4
# GSM879131 MCL_IGHV_MUT（R1304）         mcl.ighv.mut.5
# GSM879132 MCL_IGHV_MUT（R1305）         mcl.ighv.mut.6
# GSM879133 MCL_IGHV_MUT（R1628）         mcl.ighv.mut.7
# GSM879134 MCL_IGHV_MUT（R1629）         mcl.ighv.mut.8
# GSM879135 MCL_IGHV_MUT（R1341）         mcl.ighv.mut.9
# GSM879136 MCL_IGHV_MUT（R1333）         mcl.ighv.mut.10
# GSM879137 MCL_IGHV_MUT（R1585）         mcl.ighv.mut.11
# GSM879138 MCL_IGHV_MUT（R1589）         mcl.ighv.mut.12
# GSM879139 MCL_IGHV_MUT（R1591）         mcl.ighv.mut.13
# GSM879140 MCL_IGHV_UNMUT（R1595）       mcl.ighv.unmut.10
# GSM879141 MCL_IGHV_MUT（R1640）         mcl.ighv.mut.14
# GSM879142 MCL_IGHV_MUT（R1626）         mcl.ighv.mut.15
# GSM879143 MCL_IGHV_MUT（R1610）         mcl.ighv.mut.16
# GSM879144 MCL_IGHV_MUT（02_201）        mcl.ighv.mut.17
# GSM879145 MCL_IGHV_MUT（R1339）         mcl.ighv.mut.18
# GSM879146 MCL_IGHV_MUT（R1301）         mcl.ighv.mut.19
# GSM879147 MCL_IGHV_UNMUT（R1338）       mcl.ighv.unmut.11
# GSM879148 MCL_IGHV_UNMUT（R1762）       mcl.ighv.unmut.12
# GSM879149 MCL_IGHV_MUT（R1938）         mcl.ighv.mut.20
# GSM879150 MCL_IGHV_MUT（R1940）         mcl.ighv.mut.21
# GSM879151 MCL_IGHV_MUT（R1942）         mcl.ighv.mut.22
# GSM879152 MCL_IGHV_MUT（R1941）         mcl.ighv.mut.23
# GSM879153 MCL_IGHV_MUT（R1970）         mcl.ighv.mut.24
# GSM879154 MCL_IGHV_UNMUT（R1897）       mcl.ighv.unmut.13
# GSM879155 MCL_IGHV_UNMUT（R1388）       mcl.ighv.unmut.14

c.all <- c("mcl.ighv.unmut.1","mcl.ighv.unmut.2","mcl.ighv.unmut.3","mcl.ighv.unmut.4","mcl.ighv.unmut.5",
           "mcl.ighv.mut.1","mcl.ighv.mut.2",
           "mcl.ighv.unmut.6","mcl.ighv.unmut.7","mcl.ighv.unmut.8",
           "mcl.ighv.mut.3",
           "mcl.ighv.unmut.9",
           "mcl.ighv.mut.4","mcl.ighv.mut.5", "mcl.ighv.mut.6","mcl.ighv.mut.7", "mcl.ighv.mut.8","mcl.ighv.mut.9", "mcl.ighv.mut.10","mcl.ighv.mut.11", "mcl.ighv.mut.12","mcl.ighv.mut.13",
           "mcl.ighv.unmut.10",
           "mcl.ighv.mut.14","mcl.ighv.mut.15", "mcl.ighv.mut.16","mcl.ighv.mut.17","mcl.ighv.mut.18", "mcl.ighv.mut.19",
           "mcl.ighv.unmut.11","mcl.ighv.unmut.12",
           "mcl.ighv.mut.20","mcl.ighv.mut.21", "mcl.ighv.mut.22","mcl.ighv.mut.23","mcl.ighv.mut.24",
           "mcl.ighv.unmut.13","mcl.ighv.unmut.14")

c.all.status <- c("unmut","unmut","unmut","unmut","unmut",
                  "mut","mut",
                  "unmut","unmut","unmut",
                  "mut",
                  "unmut",
                  "mut","mut", "mut","mut", "mut","mut", "mut","mut", "mut","mut",
                  "unmut",
                  "mut","mut", "mut","mut","mut", "mut",
                  "unmut","unmut",
                  "mut","mut", "mut","mut","mut",
                  "unmut","unmut")
sampleNames(affy.data) <- c.all
sampleNames(affy.data)

3、 質(zhì)量分析

3.1 芯片灰度圖

Affymetrix芯片在印刷時(shí)，會(huì)在在左上角印制芯片的名稱，在四個(gè)角印制實(shí)心的花紋，我們可以通過這個(gè)特征來了解芯片數(shù)據(jù)是否可靠。比如說，圖像特別黑，說明信號(hào)強(qiáng)度低；圖像特別亮，說明信號(hào)可能過飽和，圖像亮度適中，說明信號(hào)比較可靠。

從GSE36000芯片組中，38個(gè)芯片，以其中的一個(gè)樣本GSM879118為例構(gòu)造芯片灰度圖，可以看出，四角實(shí)心的花紋名下，并且也有芯片名稱，圖像亮度適中，我們可以認(rèn)為：該樣本信號(hào)比較可靠。

GSE36000芯片組-GSM879118芯片灰度圖

GSE36000芯片組-GSM879118芯片灰度圖2

代碼參考

# 2.5讀取芯片灰度圖
# 2.5.1 設(shè)置工作目錄——放置結(jié)果
setwd("C:/libtobrain/R/R_Package/workSpaces/Version5/result2/")
# 2.5.2 將基因芯片進(jìn)行解析讀取，存放為result目錄下的pdf文件，即芯片灰度圖
pdf(file = "芯片灰度圖.pdf")
image(affy.data[,1])
dev.off()
3.2 質(zhì)量控制

采取了以上兩種方式，來評(píng)估一個(gè)芯片是否質(zhì)量可靠，還是不夠，那么我們就需要采用更加綜合的方式來采取分析。

首先我們采取了對(duì)數(shù)據(jù)進(jìn)行線性擬合，并將權(quán)重圖、殘差圖、殘差符號(hào)圖保存到pdf上（也可以直接在RStudio上觀察，通過筆記本來跑較大數(shù)據(jù)的分析，例如本次下載的數(shù)據(jù)就高到165MB，還是相對(duì)比較吃力的，個(gè)人建議，為了減少資源的占用，直接將圖片輸出會(huì)高效一些）

以本次實(shí)驗(yàn)GSE36000為例，可以看出權(quán)重圖、殘差圖、殘差符號(hào)圖，色彩平均，那么我們可以認(rèn)定該芯片質(zhì)量均勻，不會(huì)因?yàn)槲恢玫牟煌霈F(xiàn)色彩過輕過重。如果圖片中出現(xiàn)明顯的色差變化，那么就認(rèn)定這組數(shù)據(jù)質(zhì)量不達(dá)標(biāo)。

1、 權(quán)重圖


GSE36000-權(quán)重圖

2、 殘差圖


GSE36000-殘差圖

3、 殘差符號(hào)圖
GSE36000-殘差符號(hào)圖

參考代碼：

# 2.6 質(zhì)量分析與控制
# 2.6.1 質(zhì)量控制，化權(quán)重圖、殘差圖等
library(preprocessCore)
library(gcrma)
library(affyPLM)

# 對(duì)探針數(shù)據(jù)集做線性擬合
Pset <- fitPLM(affy.data)
# 根據(jù)計(jì)算結(jié)果，畫權(quán)重圖
pdf(file = "權(quán)重圖.pdf")
image(Pset, type="weights", which=1, main="Weights")
dev.off()
# 根據(jù)計(jì)算結(jié)果，畫殘差圖
pdf(file = "殘差圖.pdf")
image(Pset,  type="resids", which=1, main="Residuals")
dev.off()
# 根據(jù)計(jì)算結(jié)果，畫殘差符號(hào)圖
pdf(file = "殘差符號(hào)圖.pdf")
image(Pset, type="sign.resids", which=1, main="Residuals.sign")
dev.off()
3.3 質(zhì)量分析報(bào)告

在獲取到了基因芯片的數(shù)據(jù)，那么為了保證數(shù)據(jù)質(zhì)量達(dá)標(biāo)，我們還需要使用到Bioconductor的simpleaffy包，使用qc的方法，對(duì)Affy芯片的數(shù)據(jù)，進(jìn)行質(zhì)量評(píng)估，如果質(zhì)量評(píng)估太差，那么我們也需要放棄這一組數(shù)據(jù)。
例如，本次實(shí)驗(yàn)使用到的數(shù)據(jù)（GSE36000)，我們進(jìn)行質(zhì)量評(píng)估，可以發(fā)現(xiàn)，第一列，是我們?cè)?讀取下載部分的數(shù)據(jù)"章節(jié)，進(jìn)行重命名的樣本名稱，第二列的數(shù)據(jù)是檢出率與平均背景噪音，第三列（QC Stats），由actin（空心三角形）、（空心圓形）gapdh、（實(shí)心圓形）尺度因子組成。
從質(zhì)量評(píng)估圖中，我們可以看出，這一芯片質(zhì)量整體良好。指標(biāo)出現(xiàn)藍(lán)色表示正常，指標(biāo)紅色表示存在質(zhì)量問題，比如unmut.8芯片相對(duì)來說actin遠(yuǎn)離尺度因子指標(biāo)，那么也可以考慮棄用這類數(shù)據(jù)，如果標(biāo)識(shí)有紅色bioB表示該樣品未檢測(cè)到，那么這類數(shù)據(jù)建議不采用。
另外，在采用GSE36000之前，本次實(shí)驗(yàn)原本是想通過基因芯片，GSE47811（已廢棄），來研究健康小鼠與患癌小鼠的唾液間的關(guān)系。但是因?yàn)椴捎觅|(zhì)量評(píng)估，發(fā)現(xiàn)大量數(shù)據(jù)不可用，所以棄用了GSE47811這組芯片。
實(shí)驗(yàn)的芯片：這里我們可以剔除數(shù)據(jù)：mut3/unmut3,unmut8


質(zhì)量分析報(bào)告

這里有個(gè)反例，是前實(shí)驗(yàn)棄用的芯片（已廢棄——研究健康小鼠與患癌小鼠的唾液間，GSE47811）
PS：這也是一個(gè)篩選數(shù)據(jù)的過程，如果不合適那么最好就拋棄，重新找更加合適的數(shù)據(jù)。


反例-質(zhì)量分析報(bào)告

參考代碼：

 2.6.2 將基因芯片進(jìn)行解析讀取，獲取質(zhì)量分析報(bào)告
library(genefilter)
library(gcrma)
library(simpleaffy)
affy.qc <- qc(affy.data)
pdf(file = "質(zhì)量分析報(bào)告.pdf", width = 10, height =20)
# 圖形化顯示分析報(bào)告
plot(affy.qc)
# 剔除這些質(zhì)量不合格的數(shù)據(jù)：mut3/unmut3,unmut8
dev.off()
3.4 RLE和NUSE箱線圖

RLE和NUSE，是質(zhì)量檢測(cè)的兩種手段。
RLE表達(dá)的是一個(gè)探針組在某個(gè)樣品在探針組所有樣品表達(dá)值的中位數(shù)取對(duì)數(shù)，而所有芯片如果都趨近于0值，那么表示這個(gè)芯片質(zhì)量比較可靠；
NUSE表達(dá)的是一個(gè)探針組在某個(gè)樣品表達(dá)值的PM值標(biāo)準(zhǔn)差除以這個(gè)探針組所有樣品表達(dá)值的PM值標(biāo)準(zhǔn)差的中位數(shù)，而所有芯片如果都趨近于1值，那么表示這個(gè)芯片質(zhì)量比較可靠。
以本次實(shí)驗(yàn)實(shí)驗(yàn)對(duì)象（GSE36000）為例，RLE趨向于0，NUSE趨向于1，芯片質(zhì)量達(dá)標(biāo)。但是如果芯片質(zhì)量有問題，那么會(huì)嚴(yán)重偏離0值（RLE)、1值（NUSE），這個(gè)時(shí)候就需要進(jìn)行舍棄了。


RLE箱線圖
NUSE箱線圖

參考代碼：

# 2.6.3質(zhì)量控制，繪制RLE和NUSE箱線圖
# 載入一組顏色
library(RColorBrewer)
library(graph)
colors <- brewer.pal(12,"Set3")
# 繪制RLE箱線圖，如果數(shù)據(jù)有部分?jǐn)?shù)據(jù)遠(yuǎn)離0值，那么就需要剔除這些質(zhì)量不合格的數(shù)據(jù)，這里無須剔除
Mbox(Pset,ylim=c(-1,1),col=colors,main="RLE",las=2)
# 繪制NUSE箱線圖，如果數(shù)據(jù)有部分?jǐn)?shù)據(jù)遠(yuǎn)離1值，那么就需要剔除這些質(zhì)量不合格的數(shù)據(jù)，這里無須剔除
boxplot(Pset,ylim=c(0.95,1.2),col=colors,main="NUSE",las=2)
3.5 質(zhì)量分析（RNA降解）

芯片質(zhì)量分析，還有一種方式，那就是根據(jù)RNA降解，如果從左到右，保持一定的斜率（不與X值水平），那么表示RNA降解不嚴(yán)重，芯片質(zhì)量可靠。以本次實(shí)驗(yàn)（GSE36000）為例，RNA降解曲線斜率高，沒有呈現(xiàn)下降趨勢(shì)，那么就表示芯片質(zhì)量可靠。不需剔除數(shù)據(jù)。


芯片質(zhì)量分析-RNA降解

參考代碼：

# 2.6.4 RNA降解，質(zhì)量分析，如果圖形中5'趨向3'斜率較低，那么就需要把這些不夠明顯區(qū)別的剔除,這里看出都明顯有斜率，這里無須剔除
affy.deg <- AffyRNAdeg(affy.data)
plotAffyRNAdeg(affy.deg, col = colors)
legend("topleft", rownames(pData(affy.data)), col = colors, lwd = 1, inset = 0.05, cex = 0.5)
3.6 初步聚類分析

從聚類分析的角度來看，未突變與突變組基本上能夠比較好的分開，但這也不能判定這個(gè)實(shí)驗(yàn)是成立的。只能說明我們的研究對(duì)象，從未突變到突變有一些因素產(chǎn)生了作用，但是其中具體問題還是需要我們?nèi)ゾ唧w分析。這里我們可以簡(jiǎn)單得出一個(gè)小結(jié)論，未突變到突變的基因，組間差異大于組內(nèi)差異，基因差異較為明顯。但這里我們不考慮去除不歸類的樣品。


初步聚類分析

4、 數(shù)據(jù)預(yù)處理

數(shù)據(jù)的預(yù)處理，將質(zhì)量分析篩選下來合格的數(shù)據(jù)，進(jìn)行背景校正、標(biāo)準(zhǔn)化、匯總?cè)齻€(gè)部分的處理，獲取用來做實(shí)驗(yàn)分析的矩陣
可以通過bgcorrect.methods()、normalize.methods(affy.data)、express.summary.stat.methods()的方式，獲取能用于進(jìn)行預(yù)處理的方法，如下圖：


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

具體流程：

affy包提供了通過expresso方法來實(shí)現(xiàn)背景校正、標(biāo)準(zhǔn)化、匯總的預(yù)處理，這里需要設(shè)定不同的參數(shù)進(jìn)行具體方法實(shí)現(xiàn)。
另外，如果想更快的實(shí)現(xiàn)同樣的功能，我們也可以直接通過采取mas5，gcrma，rma方法來進(jìn)行一體化的預(yù)處理。這里mas5需要進(jìn)行對(duì)數(shù)的轉(zhuǎn)換，獲取數(shù)字化的表達(dá)譜矩陣。

經(jīng)過上面那么多步驟，這里需要停下來，對(duì)數(shù)據(jù)進(jìn)行一次處理（篩選時(shí)刻），參考代碼：

# 2.6.5 將以上的篩選重新組成實(shí)驗(yàn)組
# GSM879118 MCL_IGHV_UNMUT（R1394）       mcl.ighv.unmut.1 
# GSM879119 MCL_IGHV_UNMUT（R1328）       mcl.ighv.unmut.2
# GSM879120 MCL_IGHV_UNMUT（R1329）       mcl.ighv.unmut.3    x
# GSM879121 MCL_IGHV_UNMUT（R1306）       mcl.ighv.unmut.4
# GSM879122 MCL_IGHV_UNMUT（R1332）       mcl.ighv.unmut.5
# GSM879123 MCL_IGHV_MUT（R1399_M）       mcl.ighv.mut.1
# GSM879124 MCL_IGHV_MUT（R1400）         mcl.ighv.mut.2
# GSM879125 MCL_IGHV_UNMUT（R1587）       mcl.ighv.unmut.6
# GSM879126 MCL_IGHV_UNMUT（R1680）       mcl.ighv.unmut.7 
# GSM879127 MCL_IGHV_UNMUT（268-01-5TR）  mcl.ighv.unmut.8    x
# GSM879128 MCL_IGHV_MUT（269-01-5TR）    mcl.ighv.mut.3      x
# GSM879129 MCL_IGHV_UNMUT（043-01-4TR）  mcl.ighv.unmut.9
# GSM879130 MCL_IGHV_MUT（R1302）         mcl.ighv.mut.4 
# GSM879131 MCL_IGHV_MUT（R1304）         mcl.ighv.mut.5
# GSM879132 MCL_IGHV_MUT（R1305）         mcl.ighv.mut.6
# GSM879133 MCL_IGHV_MUT（R1628）         mcl.ighv.mut.7
# GSM879134 MCL_IGHV_MUT（R1629）         mcl.ighv.mut.8
# GSM879135 MCL_IGHV_MUT（R1341）         mcl.ighv.mut.9
# GSM879136 MCL_IGHV_MUT（R1333）         mcl.ighv.mut.10
# GSM879137 MCL_IGHV_MUT（R1585）         mcl.ighv.mut.11
# GSM879138 MCL_IGHV_MUT（R1589）         mcl.ighv.mut.12
# GSM879139 MCL_IGHV_MUT（R1591）         mcl.ighv.mut.13
# GSM879140 MCL_IGHV_UNMUT（R1595）       mcl.ighv.unmut.10
# GSM879141 MCL_IGHV_MUT（R1640）         mcl.ighv.mut.14
# GSM879142 MCL_IGHV_MUT（R1626）         mcl.ighv.mut.15
# GSM879143 MCL_IGHV_MUT（R1610）         mcl.ighv.mut.16
# GSM879144 MCL_IGHV_MUT（02_201）        mcl.ighv.mut.17
# GSM879145 MCL_IGHV_MUT（R1339）         mcl.ighv.mut.18
# GSM879146 MCL_IGHV_MUT（R1301）         mcl.ighv.mut.19
# GSM879147 MCL_IGHV_UNMUT（R1338）       mcl.ighv.unmut.11
# GSM879148 MCL_IGHV_UNMUT（R1762）       mcl.ighv.unmut.12
# GSM879149 MCL_IGHV_MUT（R1938）         mcl.ighv.mut.20
# GSM879150 MCL_IGHV_MUT（R1940）         mcl.ighv.mut.21
# GSM879151 MCL_IGHV_MUT（R1942）         mcl.ighv.mut.22
# GSM879152 MCL_IGHV_MUT（R1941）         mcl.ighv.mut.23
# GSM879153 MCL_IGHV_MUT（R1970）         mcl.ighv.mut.24
# GSM879154 MCL_IGHV_UNMUT（R1897）       mcl.ighv.unmut.13
# GSM879155 MCL_IGHV_UNMUT（R1388）       mcl.ighv.unmut.14


c.experiment  <- c("mcl.ighv.unmut.1","mcl.ighv.unmut.2","mcl.ighv.unmut.4","mcl.ighv.unmut.5","mcl.ighv.unmut.6","mcl.ighv.unmut.7","mcl.ighv.unmut.9","mcl.ighv.unmut.10",
                   "mcl.ighv.unmut.11","mcl.ighv.unmut.12", "mcl.ighv.unmut.13","mcl.ighv.unmut.14",
                   "mcl.ighv.mut.1","mcl.ighv.mut.2","mcl.ighv.mut.4","mcl.ighv.mut.5", "mcl.ighv.mut.6","mcl.ighv.mut.7", "mcl.ighv.mut.8","mcl.ighv.mut.9", "mcl.ighv.mut.10",
                   "mcl.ighv.mut.11", "mcl.ighv.mut.12","mcl.ighv.mut.13","mcl.ighv.mut.14","mcl.ighv.mut.15", "mcl.ighv.mut.16","mcl.ighv.mut.17","mcl.ighv.mut.18", "mcl.ighv.mut.19",
                   "mcl.ighv.mut.20","mcl.ighv.mut.21", "mcl.ighv.mut.22","mcl.ighv.mut.23","mcl.ighv.mut.24"
                   )

c.experiment.status <- c("unmut","unmut","unmut","unmut","unmut","unmut",
                  "unmut","unmut","unmut","unmut","unmut","unmut",
                  "mut","mut", "mut","mut", "mut","mut", "mut","mut", "mut","mut",
                  "mut","mut", "mut","mut", "mut","mut", "mut","mut", "mut","mut",
                  "mut","mut", "mut"
                  )
length(c.experiment)
length(c.experiment.status)
affy.experiment <- affy.data[,c.experiment]
colnames(affy.experiment)

# 2.7 將原始數(shù)據(jù)組成數(shù)據(jù)框，對(duì)數(shù)據(jù)進(jìn)行一個(gè)分類的整理（為差異分析做準(zhǔn)備）
affy.experiment.frame <- data.frame(SampleID = c.experiment, Disease = c.experiment.status)

sampleNames(affy.experiment)

c.unmut <- c("mcl.ighv.unmut.1","mcl.ighv.unmut.2","mcl.ighv.unmut.4","mcl.ighv.unmut.5","mcl.ighv.unmut.6","mcl.ighv.unmut.7","mcl.ighv.unmut.9","mcl.ighv.unmut.10",
             "mcl.ighv.unmut.11","mcl.ighv.unmut.12", "mcl.ighv.unmut.13","mcl.ighv.unmut.14")
c.mut <- c("mcl.ighv.mut.1","mcl.ighv.mut.2","mcl.ighv.mut.4","mcl.ighv.mut.5", "mcl.ighv.mut.6","mcl.ighv.mut.7", "mcl.ighv.mut.8","mcl.ighv.mut.9", "mcl.ighv.mut.10",
           "mcl.ighv.mut.11", "mcl.ighv.mut.12","mcl.ighv.mut.13","mcl.ighv.mut.14","mcl.ighv.mut.15", "mcl.ighv.mut.16","mcl.ighv.mut.17","mcl.ighv.mut.18", "mcl.ighv.mut.19",
           "mcl.ighv.mut.20","mcl.ighv.mut.21", "mcl.ighv.mut.22","mcl.ighv.mut.23","mcl.ighv.mut.24")
length(c.unmut)
length(c.mut)


affy.experiment.unmut <- affy.experiment[,c.unmut]
affy.experiment.mut <- affy.experiment[,c.mut]

# 2.8 預(yù)處理，（預(yù)處理）背景校正、標(biāo)準(zhǔn)化、匯總
# 預(yù)處理的方式有ms、expresso、gcrma
bgcorrect.methods()
normalize.methods(affy.experiment)
express.summary.stat.methods()

# 2.8.1 通過聚類分析，查看數(shù)據(jù)質(zhì)量，預(yù)處理
# 使用gcrma算法來預(yù)處理數(shù)據(jù)，也可以通過mas5、rma算法
affy.experiment.gcrma <- gcrma(affy.experiment)
affy.experiment.gcrma.exprs <- exprs(affy.experiment.gcrma)

# 計(jì)算樣品兩兩之間的Pearson相關(guān)系數(shù)
pearson_cor <- cor(affy.experiment.gcrma.exprs)
# 得到Pearson距離的下三角矩陣
dist.lower <- as.dist(1 - pearson_cor)
# 聚類分析、畫圖
hc <- hclust(dist.lower,"ave")

pdf(file = "聚類分析.pdf", width = 10, height =10)
plot(hc)
dev.off()

# PCA
# samplenames <- sub(pattern = "\\.CEL", replacement = "", colnames(affy.data.gcrma.eset))
groups <- factor(affy.experiment.frame[,2])
# BiocManager::install("affycoretools")
# BiocManager::install("affycoretools")
library("affycoretools")
affycoretools::plotPCA(affy.experiment.gcrma.exprs, addtext=c.experiment, groups=groups, groupnames=levels(c.experiment.status))

# 2.8.2 使用expresso，進(jìn)行背景校正、標(biāo)準(zhǔn)化處理和匯總
affy.experiment.expresso.eset <- expresso(affy.experiment, bgcorrect.method = "mas", normalize.method = "constant", pmcorrect.method = "mas",
                                          summary.method = "mas")
affy.experiment.expresso.exprs = exprs(affy.experiment.expresso.eset) 

# 2.8.3 直接采用MAS5算法進(jìn)行數(shù)據(jù)預(yù)處理
affy.experiment.mas5.eset = mas5(affy.experiment)  
# 用exprs()從eset中獲取數(shù)字化的表達(dá)譜矩陣
affy.experiment.mas5.exprs = exprs(affy.experiment.mas5.eset)

5、 基因芯片差異分析

5.1 limma處理差異分析

limma是一個(gè)比較全面的，通過用于芯片、RNA-Seq、PCR進(jìn)行線性模型和差異表達(dá)函數(shù)的工具包，它通過貝葉斯方法提供穩(wěn)定的分析結(jié)果，limma支持多重復(fù)雜的設(shè)計(jì)實(shí)驗(yàn)，我們可以進(jìn)行多個(gè)差異條件進(jìn)行分析，差異條件越多，就越復(fù)雜。除了limma進(jìn)行差異分析，也可以通過DESeq來集成處理。

例如，本次實(shí)驗(yàn)（GES36000）主要是對(duì)未突變和突變兩種差異條件進(jìn)行差異分析，是比較簡(jiǎn)單的情況。

使用limma進(jìn)行差異分析，主要需要準(zhǔn)備好：分組矩陣，差異比較矩陣，表達(dá)矩陣，進(jìn)行線性模擬擬合與差值計(jì)算、貝葉斯檢驗(yàn)、生成檢驗(yàn)結(jié)果。

具體流程：
1）獲取實(shí)驗(yàn)設(shè)計(jì)矩陣，即未突變、突變兩種情況

實(shí)驗(yàn)設(shè)計(jì)矩陣-突變、未突變

實(shí)驗(yàn)設(shè)計(jì)矩陣-突變、未突變2

2）構(gòu)建對(duì)比模型，比較兩個(gè)實(shí)驗(yàn)條件下表達(dá)數(shù)據(jù)，本次實(shí)驗(yàn)只有兩種差異條件，如果需要讀個(gè)差異條件，通過makeContrasts進(jìn)行多對(duì)多特征交互比較，得到差異


構(gòu)建對(duì)比模型

3）獲取表達(dá)矩陣，矩陣要求進(jìn)行過預(yù)處理，并為對(duì)數(shù)轉(zhuǎn)換的的表達(dá)譜矩陣
4）對(duì)表達(dá)矩陣與實(shí)驗(yàn)設(shè)計(jì)矩陣進(jìn)行線性模擬擬合
5）根據(jù)對(duì)比模型進(jìn)行差值計(jì)算
6）對(duì)差值計(jì)算的結(jié)果進(jìn)行貝葉斯檢驗(yàn)
7）生成檢驗(yàn)結(jié)果的報(bào)告
8）對(duì)P.value進(jìn)行篩選，得到比較理想的差異表達(dá)基因


未完全處理的差異表達(dá)基因

PS：當(dāng)然也可以直接使用GEO2R做一些處理，后面再做分享吧。

參考代碼：

# 基因芯片差異分析
# 3.1 差異分析方法：limma
# 3.1.1 選取差異表達(dá)基因，使用Bioconductor中的limma包
# 將（未突變、突變）狀態(tài)數(shù)據(jù)框數(shù)據(jù)讀取為因子
library(limma)
status <- factor(affy.experiment.frame[, "Disease"])
# 構(gòu)建實(shí)驗(yàn)設(shè)計(jì)矩陣（即未突變、突變兩種差異對(duì)比）
design <- model.matrix(~-1+status)
head(design)
# 構(gòu)建對(duì)比模型，比較兩個(gè)實(shí)驗(yàn)條件下表達(dá)數(shù)據(jù),其中contrasts可以參考design的列名
contrast.matrix <- makeContrasts(contrasts = "statusmut - statusunmut", levels=design)

# 3.1.2 線性模型擬合
fit <- lmFit(affy.experiment.gcrma.exprs , design)
colnames(affy.experiment.gcrma.exprs)
# 根據(jù)對(duì)比模型進(jìn)行差值計(jì)算 
fit1 <- contrasts.fit(fit, contrast.matrix)
# 3.1.3 貝葉斯檢驗(yàn)
fit2 <- eBayes(fit1)
View(fit2)

# 3.1.4 生成所有基因的檢驗(yàn)結(jié)果報(bào)告
result <- topTable(fit2, coef="statusmut - statusunmut", n=nrow(fit2), lfc=log2(2))
nrow(result)  # 1446
# 用P.Value進(jìn)行篩選，得到全部差異表達(dá)基因
result <- result[result[, "P.Value"]<0.01,]
result2 <- result[result[, "P.Value"]<0.001,]
# 顯示一部分報(bào)告結(jié)果
head(result)
head(result2)
# 分析不同p.value值下，解析出來匹配的數(shù)據(jù)行數(shù)
nrow(result)  # 909
nrow(result2) # 701
write.table(result, file="statusmut_statusunmut+0.01.txt", quote=F, sep="\t")
write.table(result2, file="statusmut_statusunmut+0.001.txt", quote=F, sep="\t")
5.2 尋找在不同條件下的差異表達(dá)

另外我們可以采用倍數(shù)法、t-test方法來進(jìn)行基因的差異表達(dá)。
affy芯片中對(duì)探針雜交質(zhì)量評(píng)定有三個(gè)值P/M/A（有/臨界值/沒有），分別是Present/Marginal/Absent。如果探針集被標(biāo)記位P或者A，表示一個(gè)芯片中，基因有表達(dá)，PM數(shù)值與MM相比，沒有顯著差異。

具體流程：
1）對(duì)表達(dá)譜矩陣進(jìn)行對(duì)數(shù)化處理（針對(duì)mas5進(jìn)行預(yù)處理的情況）
2）計(jì)算樣品組與實(shí)驗(yàn)組的平均表達(dá)值，計(jì)算差異的相對(duì)表達(dá)量
3）使用倍數(shù)法，篩選出2倍以上差異的基因
4）使用t-test檢驗(yàn)差異基因（注意：根據(jù) ?test方式，讀取文檔中提及，兩個(gè)實(shí)驗(yàn)組需要保證，相同對(duì)比組，即樣品組與實(shí)驗(yàn)組需要同數(shù)量）
5）計(jì)算p值
6）檢驗(yàn)每張芯片中的三個(gè)值P/M/A（Present/Marginal/Absent），篩選出至少一個(gè)有表達(dá)的部分
7）設(shè)置FDR（False Discovery Rate）方法，來校正P值
8）篩選出2倍以上差異，P值小于0.05的，至少一個(gè)有表達(dá)的基因
9）構(gòu)造差異表達(dá)基因的熱圖

差異表達(dá)基因的熱圖

參考代碼

# 3.2 尋找在不同條件下存在差異表達(dá)的基因
# 3.2.1 對(duì)表達(dá)譜矩陣進(jìn)行對(duì)數(shù)化處理
affy.experiment.mas5.exprs.log <- log(affy.experiment.mas5.exprs, 2)
probeset.id <- row.names(affy.experiment.mas5.exprs.log)
# 保存數(shù)據(jù)
write.table(affy.experiment.mas5.exprs.log, file="affy.experiment.mas5.exprs.log.txt", quote=F, sep="\t")

# 3.2.2 計(jì)算每個(gè)基因在樣本中的平均表達(dá)值
unmut.mean <- apply(affy.experiment.mas5.exprs.log[, c.unmut], 1, mean)
mut.mean <- apply(affy.experiment.mas5.exprs.log[, c.mut], 1, mean)
# 計(jì)算差異基因的相對(duì)表達(dá)量 (log(fold change))
unmut.to.mut.mean <- unmut.mean - mut.mean

# 使用cbind，將數(shù)據(jù)表拼接在一起
unmut.to.mut.mean.data <- cbind(affy.experiment.mas5.exprs.log, unmut.mean, mut.mean, unmut.to.mut.mean)
head(unmut.to.mut.mean.data)
# 保存數(shù)據(jù)
write.table(unmut.to.mut.mean.data, file="unmut.to.mut.mean.data.txt", quote=F, sep="\t")

# 3.2.3 倍數(shù)法，尋找在不同條件下表達(dá)量存在2倍以上差異的基因(探針組) (log(fold change)>1 or <-1)
unmut.to.mut.mean.fc.probesets <- names(unmut.to.mut.mean[unmut.to.mut.mean >1 | unmut.to.mut.mean<(-1)])

# 3.2.4 用t test檢驗(yàn)?zāi)硞€(gè)基因在突變與未突變基因中的表達(dá)是否存在差異
dataset.unmut <- affy.experiment.mas5.exprs.log[1, c.unmut]
dataset.mut <- affy.experiment.mas5.exprs.log[1, c.mut]
test.gene <- t.test(dataset.unmut, dataset.mut, "two.sided")
test.gene
test.gene$p.value
length(dataset.mut)
length(dataset.unmut)
?t.test

# 3.2.5 使用apply函數(shù)，對(duì)t test檢驗(yàn)每個(gè)基因在未突變與突變基因中的表達(dá)是否存在差異,計(jì)算p值
dim(affy.experiment.mas5.exprs.log)
p.value.experiment.genes <- apply(affy.experiment.mas5.exprs.log, 1, function(x) { t.test(x[1:12], x[13:35]) $p.value } )

# 3.2.6 用mas5calls函數(shù)來檢測(cè)每張芯片中每個(gè)基因是否表達(dá),P/M/A（有/臨界值/沒有）
affy.experiment.mas5calls <- mas5calls(affy.experiment)
affy.experiment.mas5calls.exprs = exprs(affy.experiment.mas5calls)
# 保存數(shù)據(jù)
write.table(affy.experiment.mas5calls.exprs, file="affy.experiment.mas5calls.exprs.txt", quote=F, sep="\t")
affy.experiment.PMA <- apply(affy.experiment.mas5calls.exprs, 1, paste, collapse="")
head(affy.experiment.PMA)

# 3.2.7 把至少在一個(gè)芯片中有表達(dá)的基因選出來
genes.present <- names(affy.experiment.PMA[affy.experiment.PMA != "AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA"])

# 54675
length(affy.experiment.PMA)
# 40673
length(genes.present)
affy.experiment.mas5.exprs.log.present <- affy.experiment.mas5.exprs.log[genes.present,]

# 3.2.8 設(shè)置FDR（False Discovery Rate）方法，來校正P值
p.adjust.methods
raw.pvals.present <- p.value.experiment.genes[genes.present]
fdr.pvals.present <- p.adjust(raw.pvals.present, method="fdr")

# 3.2.9 對(duì)FDR p值按從小到大排序
fdr.pvals.present.sorted <- fdr.pvals.present[order(fdr.pvals.present)]
length(fdr.pvals.present.sorted)
# 輸出FDR p值最小的10個(gè)基因
fdr.pvals.present.sorted[1:10]
fdr.pvals.present.sorted[40663:40673]
# 將結(jié)果整理成表格
expression.plus.pvals <- cbind(affy.experiment.mas5.exprs.log.present, raw.pvals.present, fdr.pvals.present)
write.table(expression.plus.pvals, "mas5_expression.plus.pvals.txt", sep="\t", quote=F)

# 3.2.10 找出那些至少在一個(gè)芯片中有表達(dá)，且在不同組織中表達(dá)有差異（FDR p值<0.23)的基因(探針組)
DE.fdr.probesets <- names(fdr.pvals.present[fdr.pvals.present < 0.05])
DE.fc.probesets <- names(unmut.to.mut.mean[unmut.to.mut.mean >1 | unmut.to.mut.mean<(-1)])
DE.probesets <- intersect(DE.fdr.probesets,DE.fc.probesets)

# 獲取那些至少在一個(gè)芯片中有表達(dá)且在不同組織中表達(dá)有差異的基因的表達(dá)量
unmut.to.mut.mean.mas5_DE <- unmut.to.mut.mean.data[DE.probesets, c.experiment]

# 3.2.11 構(gòu)造熱圖
library(pheatmap)
pdf(file = "差異表達(dá)基因熱圖.pdf",width = 14,height = 35)
# pheatmap(unmut.to.mut.mean.mas5_DE,color=colorRampPalette(c("green","black","red"))(100),clustering_distance_rows = "correlation",clustering_distance_cols = "euclidean",clustering_method="complete")
pheatmap(unmut.to.mut.mean.mas5_DE,color=colorRampPalette(c("green","black","red"))(100),
         clustering_distance_rows = "correlation",clustering_distance_cols = "euclidean",clustering_method="complete")
dev.off()

rnames.mas5_DE <- as.matrix(rownames(unmut.to.mut.mean.mas5_DE))
colnames(rnames.mas5_DE) <- "probeset.id"

affy.experiment.mas5.DE.mat <- cbind(rnames.mas5_DE,unmut.to.mut.mean.mas5_DE)
#保存為txt文件
write.table(affy.experiment.mas5.DE.mat, "unmut.to.mut.mean.mas5_DE.txt", sep="\t", quote=F,row.names=FALSE)
5.3 小結(jié)

本次實(shí)驗(yàn)中，對(duì)P.value的一些感想：
P.value表示假設(shè)成立的條件下，重復(fù)n次試驗(yàn)，獲得現(xiàn)有統(tǒng)計(jì)量以及極端情況下的概率。在樣本量較低，使用P.value值可靠性會(huì)比較低，因?yàn)椴町惙譃榻M內(nèi)差異與組間差異（樣品組與實(shí)驗(yàn)組），如果樣本量少，那么組內(nèi)差異會(huì)比較大，那么就會(huì)影響到我們需要關(guān)注的組間差異。
P.value代表著閾值，我們可以提高閾值的方式（即將P.value<0.05轉(zhuǎn)為P.value<0.01），雖然這降低了假陽(yáng)性的概率，不過通過這種方式設(shè)定標(biāo)準(zhǔn)，會(huì)讓原本可能真的存在表達(dá)差異的數(shù)據(jù)達(dá)不到我們?cè)O(shè)定的0.01的閾值標(biāo)準(zhǔn)，這就會(huì)提高了假陰性率。
芯片數(shù)據(jù)分析，熱圖，是直觀展示出基因表達(dá)變化的二維土層，其中數(shù)值以紅綠黑為基本色調(diào)，能夠讓繁雜的數(shù)據(jù)一目了然，畢竟圖片比文字更能夠直觀的表達(dá)。
在本次差異分析中，GES36000樣本，右側(cè)文本為基因編號(hào)，熱圖的每一行都是代表這一個(gè)基因，下方文本為樣本編號(hào)，熱圖的每一列都是代表著一個(gè)樣本，左側(cè)為樹狀結(jié)構(gòu)，代表著基因之間的相似度聚落關(guān)系，頂部也為樹狀結(jié)構(gòu)，表示樣本之間的相似度聚落關(guān)系。其中，紅色表示上調(diào)，綠色表示下調(diào)
這里我們可以看出，GES36000樣本，突變與未突變的樣本在一些基因上的差異表達(dá)還是比較明顯的。

6、 基因功能分析

6.1 注釋工具包

通過bioconductor收錄的芯片注釋包，將對(duì)應(yīng)ID值進(jìn)行映射，得到Gene symbol和Entrez ID兩種探針注釋信息，有了這些注釋，能夠直觀的了解每個(gè)ID值對(duì)應(yīng)的是什么術(shù)語(yǔ)？確定每個(gè)探針對(duì)應(yīng)檢測(cè)哪個(gè)基因的表達(dá)。

在添加注釋工具包時(shí)，可以先通過show的方式來自動(dòng)下載對(duì)應(yīng)的注釋包，例如，本次實(shí)驗(yàn)中是引用了hgu133plus2.db注釋包

參考代碼：

# 4.基因功能分析
# 4.1 注釋工具包
# 4.1.1 加載注釋工具包
# BiocManager::install(affydb)
library(XML)
library(annotate)
library(org.Hs.eg.db)

# 4.1.2 獲得基因芯片注釋包名稱
affy.experiment@annotation
# 4.1.3 加載注釋包"hgu133plus2.db"
affy.experiment.affydb <- annPkgName(affy.experiment@annotation, type="db")
library(affy.experiment.affydb, character.only=TRUE)
# 根據(jù)每個(gè)探針組的ID,獲取那些至少在一個(gè)芯片中有表達(dá)且在不同組織中表達(dá)有差異的基因名、Entrez ID
unmut.to.mut.mean.symbols <- getSYMBOL(rownames(unmut.to.mut.mean.mas5_DE),affy.experiment.affydb)
unmut.to.mut.mean.EntrezID <- getEG(rownames(unmut.to.mut.mean.mas5_DE),affy.experiment.affydb)
unmut.to.mut.mean.mas5.DE.anno=cbind(unmut.to.mut.mean.EntrezID,unmut.to.mut.mean.symbols,unmut.to.mut.mean.mas5_DE)
# 保存為txt文件
write.table(unmut.to.mut.mean.mas5.DE.anno, "unmut.to.mut.mean.mas5.DE.anno.txt", sep="\t", quote=F,row.names = FALSE)

根據(jù)挑選出來的差異基因，GO富集分析會(huì)對(duì)實(shí)驗(yàn)結(jié)果有補(bǔ)充提示的作用，我們可以通過GO富集分析差異基因，找到富集差異基因的GO分類目錄，其中數(shù)據(jù)結(jié)果可以保存為本文，同時(shí)也可以保存為html的方式，這提供了較大的便利，我們可以通過這些方式，更加直觀的尋找不同樣品的差異基因可能和哪一些基因功能的改變有關(guān)系


GO富集分析

參考代碼：

# 5 GO富集分析
# 加載所需R包
library(Matrix)
library(Category)
library(graph)
library(GOstats)
library("GO.db")
# 5.1 獲取基因芯片所有探針組與差異表達(dá)基因的EntrezID
# 提取芯片中affy.data所有探針組對(duì)應(yīng)的EntrezID，注意保證uniq
entrezUniverse <- unique(getEG(rownames(affy.experiment),affy.experiment.affydb));
# 提取所有差異表達(dá)基因及其對(duì)應(yīng)的EntrezID，去除na值，注意保證uniq
entrezSelected <- unique(unmut.to.mut.mean.EntrezID[!is.na(unmut.to.mut.mean.EntrezID)]);
# 設(shè)置GO富集分析的所有參數(shù)
params <- new("GOHyperGParams", geneIds = entrezSelected, universeGeneIds = entrezUniverse, 
              annotation = affy.experiment.affydb, ontology = "BP", pvalueCutoff = 0.001, conditional = FALSE, testDirection = "over");
# 5.2 對(duì)所有的GO term根據(jù)params參數(shù)做超幾何檢驗(yàn)，并進(jìn)行匯總
hgOver <- hyperGTest(params);
bp <- summary(hgOver) ;

# 5.3 同時(shí)生成所有GO term的檢驗(yàn)結(jié)果文件，每個(gè)GOterm都有指向官方網(wǎng)站的鏈接，可以獲得其詳細(xì)信息
htmlReport(hgOver,  file='unmut.to.mut.mean.mas5_DE_go.html') ;
head(bp)
# 保存數(shù)據(jù)
write.table(bp, "unmut.to.mut.mean.mas5_DE_go.txt", sep="\t", quote=F,row.names = FALSE)

# 根據(jù)每個(gè)探針組的ID獲取對(duì)應(yīng)基因Gene Symbol，并作為新的一列
result$symbols <- getSYMBOL(rownames(result), affy.experiment.affydb)
# 根據(jù)探針I(yè)D獲取對(duì)應(yīng)基因Entrez ID
result$EntrezID <- getEG(rownames(result), affy.experiment.affydb)
# 顯示前幾行
head(result)
nrow(result)

四、結(jié)論

通過在NCBI下載一套人類基因（未突變和突變淋巴），其中包含了38組樣品基因，其中未突變14組，突變24組，前前后后對(duì)這38組數(shù)據(jù)編輯了5套實(shí)驗(yàn)?zāi)_本，進(jìn)行數(shù)據(jù)的載入、讀取、質(zhì)量控制、質(zhì)量分析、聚類分析、預(yù)處理、差異分析、注解、生成結(jié)果分析報(bào)告等。

整一套流程走了幾遍，主要收獲如下：基因芯片的差異分析，主要分為數(shù)據(jù)的處理與分析；大多數(shù)人會(huì)側(cè)重在分析上，也為了得到研究想要獲得的結(jié)果；但是數(shù)據(jù)的篩選與質(zhì)量分析也是一個(gè)很重要的部分。

在本次研究之前（已廢棄），是以研究健康小鼠與患癌小鼠的唾液間的關(guān)系，但因?yàn)榛蛸|(zhì)量不達(dá)標(biāo)，選擇放棄這組實(shí)驗(yàn)，另外在基因篩選上發(fā)現(xiàn)，如果樣品組與實(shí)驗(yàn)組（突變與未突變）統(tǒng)計(jì)量少（38組數(shù)據(jù)，只人工采用20組數(shù)據(jù)），亦或是在p值設(shè)置閾值偏低，那么組內(nèi)的差異很可能因?yàn)槲覀儗?shí)驗(yàn)設(shè)計(jì)的不合理，而升高，導(dǎo)致實(shí)驗(yàn)結(jié)果十分不理想（前后5組實(shí)驗(yàn)，對(duì)比分析結(jié)果）。

所以，我們?cè)谔幚頂?shù)據(jù)，進(jìn)行基因差異分析，需要保證我們的結(jié)果的可靠性，那么就需要我們注重我們?cè)诜治鰰r(shí)的嚴(yán)謹(jǐn)性，并反復(fù)測(cè)試不同實(shí)驗(yàn)變量，對(duì)比結(jié)果，得出合理的結(jié)論。

五、參考資料

[1] NCBI NCBI介紹[EB/OL]https://www.ncbi.nlm.nih.gov/home/about/
[2] W3Cschool R介紹[EB/OL]https://www.w3cschool.cn/r/
[3] 朱猛進(jìn).差異表達(dá)基因的分析[CP/OL]http://blog.sciencenet.cn/blog-295006-403640.html
[4] 顯著檢驗(yàn)與多重檢驗(yàn)校正[S/OL] http://www.omicshare.com/forum/thread-260-1-12.html
[5] tommyhechina芯片數(shù)據(jù)分析（探針注釋）[S/OL] https://blog.csdn.net/tommyhechina/article/details/80409983

六、本文源碼

github的repo地址：https://github.com/PinkSmallFan/Pbioinformation

對(duì)應(yīng)本文的R源碼：https://github.com/PinkSmallFan/Pbioinformation/tree/master/source_code/differential_analysis/Human_unmut_mut

PS：本文產(chǎn)生的源數(shù)據(jù)（例如芯片灰度圖、箱線圖、差異分析、聚集分析、PCA等等的數(shù)據(jù)與對(duì)應(yīng)的原圖這里就不共享啦~ 需要的話可以留言，酌情考慮）
PS2：如果這篇筆記有什么不足，或者疑惑不解的地方（可能對(duì)小白不太友好），可以提出來，看到會(huì)回答（留個(gè)TODO給大家~）
如果感覺對(duì)你有幫助，可以關(guān)注：專欄-生物信息學(xué)-小白成長(zhǎng)記

最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡(jiǎn)書系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容