我們在上幾期已經介紹完R語言的基礎及統計學方法，今天小編將給大家整理了用R語言進行芯片數據預處理、數據分析、可視化分析。
01芯片數據預處理（affy包）
基因芯片技術是使用寡聚核苷酸探針檢測基因ab 結果分析存儲到哪，比方說我們用函數讀取CEL文件獲得的數據是探針水平的，即雜交信號，而芯片數據預處理的目的是將雜交信號轉換成表達水平數據?；蛐酒结標綌祿幚淼腞包有affy、、、gcrma等。
library(BiocInstaller)biocLite(c("affy","gcrma","affyPLM","affycomp"))
affy包芯片數據預處理有三個步驟：
背景處理
歸一化處理
匯總
Affy基因芯片的探針長度是25個堿基，每個mRNA用11-20個探針檢測，檢測同一個mRNA的一組探針稱為probe sets。因為探針較短，為保證雜交的特異性，affy公司為每個基因設計了兩類探針，一類探針的序列與基因完全匹配，稱為 match (PM) ；另一類為不匹配探針，稱為 (MM) 。PM和MM探針序列除了第13個堿基外其他都一樣，MM中把PM第13個堿基換為互補堿基。PM和MM探針成對出現。
1.1
載入芯片數據、修改芯片名稱
library(affy)library(tcltk)filters <- matrix(c("CEL file", ".[Cc][Ee][Ll]", "All", ".*"), ncol = 2, byrow = T)cel.files <- tk_choose.files(caption = "Select CELs", multi = TRUE, filters = filters, index = 1)##這里查看一下文件名稱basename(col.files)## [1] "NRID9780_Zarka_2-1_MT-0HCA(SOIL)_Rep1_ATH1.CEL"## [2] "NRID9781_Zarka_2-2_MT-0HCB(SOIL)_Rep2_ATH1.CEL"## [3] "NRID9782_Zarka_2-3_MT-1HCA(SOIL)_Rep1_ATH1.CEL"## [4] "NRID9783_Zarka_2-4_MT-1HCB(SOIL)_Rep2_ATH1.CEL"## [5] "NRID9784_Zarka_2-5_MT-24HCA(SOIL)_Rep1_ATH1.CEL"## [6] "NRID9785_Zarka_2-6_MT-24HCB(SOIL)_Rep2_ATH1.CEL"## [7] "NRID9786_Zarka_2-7_MT-7DCA(SOIL)_Rep1_ATH1.CEL"## [8] "NRID9787_Zarka_2-8_MT-7DCB(SOIL)_Rep2_ATH1.CEL"data.raw <- ReadAffy(filenames = cel.files)sampleNames(data.raw) <-paste("CHIP",1:length(cel.files),sep="")##用pm和mm函數可查看每個探針的檢測情況：pm.data.raw <- pm(data.raw)head(pm.data.raw, 2)## CHIP1 CHIP2 CHIP3 CHIP4 CHIP5 CHIP6 CHIP7 CHIP8## 501131 127.0 166.3 112 139.8 111.3 85.5 126.3 102.8## 251604 118.5 105.0 82 101.5 94.0 81.3 103.8 103.0mm.data.raw <- mm(data.raw)head(mm.data.raw, 2)## CHIP1 CHIP2 CHIP3 CHIP4 CHIP5 CHIP6 CHIP7 CHIP8## 501843 89.0 88.0 80.5 91.0 77.0 75 79.0 72.0##?252316?134.3??77.3??77.0?107.8??98.5????75??99.5??71.3
1.2
背景處理
這一步要處理背景值和噪聲信號。注意把背景和噪聲區分開。
Affy包中用于消減背景噪聲的函數是bg.()，可用MAS和RMA方法。
MAS方法將芯片分為k個網格區域，用每個區域使用信號強度最低的2%探針計算背景值和噪聲。
data.rma <- bg.correct(data.raw, method="rma")data.mas?<-?bg.correct(data.raw,?method="mas")class(data.rma)## [1] "AffyBatch"## attr(,"package")## [1] "affy"class(data.mas)## [1] "AffyBatch"## attr(,"package")## [1] "affy"class(data.raw)## [1] "AffyBatch"## attr(,"package")## [1] "affy"
()讀入的CEL芯片數據以類數據形式存儲，而背景消減后得到的是類數據。
MAS方法應用后PM和MM的信號強度被重新計算。而RMA方法只用在PM探針數據，背景調整后MM信號值不變。
head(pm(data.raw)-pm(data.mas),2)## CHIP1 CHIP2 CHIP3 CHIP4 CHIP5 CHIP6 CHIP7 CHIP8## 501131 79.34 62.93 63.23 72.87 62.48 64.43 62.97 60.86## 251604 77.57 63.06 63.73 80.69 63.07 66.53 63.33 60.34head(pm(data.raw)-pm(data.rma),2)## CHIP1 CHIP2 CHIP3 CHIP4 CHIP5 CHIP6 CHIP7 CHIP8## 501131 111.2 102.21 93.23 116.36 93.75 76.15 102.82 85.21## 251604 103.9 88.69 72.57 90.75 82.23 72.64 87.75 85.32head(mm(data.raw)-mm(data.mas),2)## CHIP1 CHIP2 CHIP3 CHIP4 CHIP5 CHIP6 CHIP7 CHIP8## 501843 79.34 62.93 63.24 72.90 62.48 64.44 62.97 60.85##?252316?77.56?63.06?63.73?80.66?63.07?66.51?63.33?60.34##差值為全部為0，說明rma方法沒有對mm數值進行處理head(mm(data.raw)-mm(data.rma),2)## CHIP1 CHIP2 CHIP3 CHIP4 CHIP5 CHIP6 CHIP7 CHIP8## 501843 0 0 0 0 0 0 0 0## 252316 0 0 0 0 0 0 0 0identical(mm(data.raw), mm(data.rma))## [1] TRUE
1.3
歸一化處理
同一個RNA樣品用相同類型的芯片進行雜交所得的結果不可能完全一樣。為了使不同芯片所得到的結果具有可比性，必須進行歸一化處理。歸一化處理的affy函數為()，以對象和處理方法為參數。
Affy公司在affy軟件中使用的是線性縮放方法，先選一個芯片作為參考，將其他芯片和參考芯片逐一做線性回歸分析，用沒有截距的回歸直線對其他芯片的信號值做縮放。Affy軟件做回歸分析前去除了2%最強和最弱的信號。
Mas方法做背景消減后進行歸一化分析的R語句是：
data.mas.ln <- normalize(data.mas, method = "constant")head(pm(data.mas.ln)/pm(data.mas), 5)head(mm(data.mas.ln)/mm(data.mas), 5)
1.4
匯總
這一步使用合適的統計方法，通過（包含多個探針）的雜交信號計算出表達量。Affy包的函數是。
兩個參數可設定的值通過下面函數獲得：
pmcorrect.methods()generateExprSet.methods()
常用的匯總方法包括、和mas。方法只用PM做背景校正，
pmcorrect.method=”pmonly”02基因表達譜芯片數據分析
基因芯片數據分析是對從基因芯片高密度雜交點陣圖中提取的雜交點熒光強度信號進行的定量分析，通過有效數據的篩選和相關基因表達譜的聚類，最終整合雜交點的生物學信息，發現基因的表達譜與功能可能存在的練習。
基因芯片數據分析方法分為
差異基因表達分析：主要篩選出不同條件下表達明顯差異的基因?？捎肎EO2R在線工具幫助尋找差異基因。用標準的統計學方法如t檢驗、方差分析、回歸分析等方法進行差異基因表達分析。（操作步驟可參閱上一期的。）
聚類分析：分析基因或樣本之間的相互關系。常用聚類方法包括分層聚類和K-均值聚類。
判別分析：以某些在不同樣品中表達差異顯著的基因作為模板，通過判別分析建立有效的疾病診斷方法。目前使用的判別分析方法主要包括：支持向量機、決策樹、貝葉斯分類、神經網絡法。
03分析數據結果的展示方式
3.1
Venn圖構建（包）
包提供了2到5個集合的繪圖函數：包括繪制兩個集合的韋恩圖的draw..venn，三個集合的draw..venn，四個集合的draw.quad.venn，五個集合的函數draw..venn。
安裝并加載包
Install.packages(“VennDiagram”)library(VennDiagram)
包中畫維恩圖的函數是venn.(). 下面我們來看一下venn.()函數的使用及參數說明。
venn.diagram(x, filename, height = 3000, width = 3000, resolution =500, imagetype = "tiff", units = "px", compression ="lzw", na = "stop", main = NULL, sub = NULL, main.pos= c(0.5, 1.05), main.fontface = "plain",main.fontfamily = "serif", main.col = "black",main.cex = 1, main.just = c(0.5, 1), sub.pos = c(0.5,1.05), sub.fontface = "plain", sub.fontfamily ="serif", sub.col = "black", sub.cex = 1, sub.just =c(0.5, 1), category.names = names(x), force.unique =TRUE, print.mode = "raw", sigdigs = 3, direct.area =FALSE, area.vector = 0, hyper.test = FALSE, total.population = NULL,????lower.tail?=?TRUE,?...)
參數說明
x: a list of , e.g:list(A=1:10, B=3:8, C=5:13)
:設置圖形輸出文件名
:輸出圖形的清晰度，DPI數值
:輸出圖形的格式，tiff, png, svg 等
alpha:設置每個區塊的透明度
main:圖形標題
main.:字體樣式，比如斜體，粗體等
main.:字體，比如Time New Roman等
利用函數draw..venn繪制兩個集合的韋恩圖
draw.pairwise.venn(area1=Length_A,area2=Length_B,cross.area=Length_AB,category=c('A','B'),參數說明
lwd=rep(1,1)##邊框線寬度
lty=rep(2,2)##lty設定圓弧的線形
col=c('blue','red') ##邊框顏色
alpha= 0.7 #透明度,
fill=c('red','green') #填充色 ,
cex= 1 #標簽字體大小,
= 0.01 #邊際距離,
cat.col=c('red','green') ##表示集合名稱的顯示顏色,r
.=180) ##調整圖形的旋轉角度
draw..venn函數繪制三個集合的韋恩圖
draw.triple.venn(area1=Length_A,?area2=Length_B,?area3=Length_C???##area1、area2、area3分別指第一個、第二個、第三個集合的大小n12=Length_AB,?n23=Length_BC,?n13=Length_AC,?n123=Length_ABC???##n12表示第一個與第二個集合的交集大小category?=?c('A','B','C')???#指定集合名稱col=c('red','green','blue')??#邊框顏色fill=c('red','green','blue') #填充色cat.col=c('red','green','blue') #表示集合名稱的顯示顏色reverse = FALSE) #指是否對圖形進行反轉
利用draw.quad.venn繪制四個集合的韋恩圖
draw.quad.venn(area1=Length_A, area2=Length_B, area3=Length_C, area4=Length_D,n12=Length_AB, n13=Length_AC, n23=Length_BC, n14=Length_AD, n24=Length_BD , n34=Length_CD,n123=Length_ABC, n124=Length_ABD, n234=Length_BCD, n134=Length_ACD,n1234=Length_ABCD ,category = c('A','B','C','D') ,col=c('red','green','blue','orange'), fill=c('red','green','blue','orange') ,alpha = 0.7 ,lty=rep(4), cat.col=c('pink') ,reverse = FALSE)
3.2
GO/KEGG富集分析（包）
包主要做GO和KEGG的功能富集及其可視化。
但需要注意的是，GO的注釋信息來自，提供19個物種的org類型的GO注釋信息。
常用的分析方法包括：
過表征分析ORA：找出差異基因后，對差異基因進行GO/KEGG富集。
基因富集分析GSEA：對所有基因進行富集，綜合考慮每個基因對基因集的貢獻。
GO是基因本體論聯合會建立的一個數據庫，旨在建立一個適用于各種物種的、對基因和蛋白功能進行限定和描述的，并能隨著研究不斷深入而更新的語義詞匯標準。
GO注釋的三大類：
分子生物學功能（ ，MF）
生物學過程（ ，BP）
細胞學組分（ ，CC）
KEGG是一個整合了基因組、化學和系統功能信息的綜合數據庫。KEGG有4個大類和17和子數據庫，其中KEGG 數據庫最為常用。
3.2.1 GO分析
包中的()函數可以對差異基因進行GO分析。
常用的GO-基因轉換的ID：
：歐洲生物信息數據庫提供的ENSG+11位數字
：NCBI提供，純數字
：基因名
：NCBI提供的參考序列數據ab 結果分析存儲到哪，以NG/NM等開頭
做GO分析時最好使用格式的ID，通過bitr()函數轉換基因ID。
bitr(geneID, fromType=””, toType=””, OrgDb=, drop=TRUE)參數說明
：需要轉換的基因ID文件
：文件中ID類型
：輸出的ID類型
OrgDb：注釋對象
Drop：是否去除空值
data(geneList, package="DOSE")head(geneList)##選取差異基因gene <- names(geneList)[abs(geneList) >2]轉換IDgeneID <- bitr(gene, fromType = "ENTREZID",toType = c("ENSEMBL", "SYMBOL"),OrgDb = org.Hs.eg.db)head(geneID)
使用()函數可以對DEG基因進行GO富集分析
gene.GO <- enrichGO(gene = geneID$ENTREZID,OrgDb = org.Hs.eg.db,keyType = "ENTREZID",ont = "CC",pAdjustMethod = "BH",pvalueCutoff = 0.01,qvalueCutoff = 0.05,readable = T)gene.GO[1:2,]##可保存生成的 @result 文件write.table(gene.GO@result, file = "geneGO.txt",??????????????sep?=?"\t",?row.names?=?T,col.names?=?T)
3.2.2KEGG分析
只用()函數即可。
使用參數指定物種，需要查找物種對應的
earch_kegg_organism("human", by="common_name")gene.KEGG <- enrichKEGG(gene = geneID$ENTREZID,organism = "hsa",keyType = "kegg",pAdjustMethod = "BH",pvalueCutoff = 0.01,????????????????????????qvalueCutoff?=?0.05?)
3.2.3GO/KEGG分析結果的可視化
可使用包進行富集結果可視化分析。
/
library(enrichplot)p1 <- barplot(gene.GO, showCategory = 10)p2 <- dotplot(gene.GO, showCategory = 10)plot_grid(p1,p2,ncol = 2)
heatplot(gene.GO, foldChange = geneList)顯示GO基因集中所富集到的基因，則將該基因與GO集連線。
cnetplot(gene.GO, showCategory = 5)顯示富集到的GO基因集之間的關系。若不同基因集見存在重疊基因，則基因集間將會連線。
emapplot(gene.GO, showCategory = 20)/
goplot(gene.ID)3.3
生存曲線分析（包）
生存分析：研究各個因素與生存時間有無關系及關聯程度大小。
生存分析的變量包括：
生存曲線是以生存時間為橫軸、生存率為縱軸繪制一條生存曲線。
生存分析的部分研究內容包括：
描述生存過程：研究生存時間的分布特點，估算生存率和標準誤，繪制生存曲線。常用乘積極限法和壽命表法。如果生存曲線為單因素分析，用中位生存時間表示生存時間的平均水平。
比較生存過程：常用log-rank檢驗
影響生存時間的因素分析：常用多因素生存分析方法包括Cox比例風險回歸模型。
在這里，我們介紹如何用R-包進行生存分析，并繪制KM曲線圖。
##載入R包，讀取數據library(survival)dat <- read.table('Test.txt',header=TRUE)##估計生存函數，觀察不同組間的區別fit <- survfit(Surv(dat$OS_time,as.numeric(dat$OS_Status))~dat$COV1,data=dat)##獲得的survial列就是生存率summary(fit)##比較不同因子分組的生存效果，檢驗顯著性survdiff(Surv(dat$OS_time,as.numeric(dat$OS_Status))~dat$COV1,data=dat)##繪制KM曲線圖plot(fit,xlab="Time(Days)",ylab="Survival",main="",col=c("blue","red"),lty=2,lwd=2)legend("topright",c("A","B"),col=c("blue","red"),lty=2,lwd=2,cex=0.7)
橫軸表示生存時間，縱軸表示生存概率，為一條梯形下降的曲線。下降幅度越大，表示生存率越低或生存時間越短。
好了，今天就介紹到這里，小編在下一期將給大家分享生物軟件。
希望我們分享的文章能幫助你發表更多SCI論文。如果你有疑問，歡迎在下面評論區留言。
溫馨提示
如果你喜歡本文，請分享到朋友圈，想要獲得更多信息，請關注我。
南博屹相伴，科研不孤單