R語言介紹

R語言是用于統(tǒng)計(jì)分析，圖形表示和報(bào)告的編程語言和軟件環(huán)境。 R語言由Ross Ihaka和Robert Gentleman在新西蘭奧克蘭大學(xué)創(chuàng)建，目前由R語言開發(fā)核心團(tuán)隊(duì)開發(fā)。

R語言的核心是解釋計(jì)算機(jī)語言，其允許分支和循環(huán)以及使用函數(shù)的模塊化編程。 R語言允許與以C，C ++，.Net，Python或FORTRAN語言編寫的過程集成以提高效率。

R語言在GNU通用公共許可證下免費(fèi)提供，并為各種操作系統(tǒng)（如Linux，Windows和Mac）提供預(yù)編譯的二進(jìn)制版本。
R是一個(gè)在GNU風(fēng)格的副本左側(cè)的自由軟件，GNU項(xiàng)目的官方部分叫做GNU S.

R的演變

R語言最初是由新西蘭奧克蘭奧克蘭大學(xué)統(tǒng)計(jì)系的Ross Ihaka和Robert Gentleman寫的。 R語言于1993年首次亮相。
一大群人通過發(fā)送代碼和錯(cuò)誤報(bào)告對R做出了貢獻(xiàn)。

自1997年年中以來，已經(jīng)有一個(gè)核心組（“R核心團(tuán)隊(duì)”）可以修改R源代碼歸檔。

R的特點(diǎn)

如前所述，R語言是用于統(tǒng)計(jì)分析，圖形表示和報(bào)告的編程語言和軟件環(huán)境。 以下是R語言的重要特點(diǎn)：

• R語言是一種開發(fā)良好，簡單有效的編程語言，包括條件，循環(huán)，用戶定義的遞歸函數(shù)以及輸入和輸出設(shè)施。
• R語言具有有效的數(shù)據(jù)處理和存儲(chǔ)設(shè)施，
• R語言提供了一套用于數(shù)組，列表，向量和矩陣計(jì)算的運(yùn)算符。
• R語言為數(shù)據(jù)分析提供了大型，一致和集成的工具集合。
• R語言提供直接在計(jì)算機(jī)上或在紙張上打印的圖形設(shè)施用于數(shù)據(jù)分析和顯示。

作為結(jié)論，R語言是世界上最廣泛使用的統(tǒng)計(jì)編程語言。 它是數(shù)據(jù)科學(xué)家的第一選擇，并由一個(gè)充滿活力和有才華的貢獻(xiàn)者社區(qū)支持。 R語言在大學(xué)教授并部署在關(guān)鍵業(yè)務(wù)應(yīng)用程序中。 本教程將教您R編程與適當(dāng)?shù)睦釉诤唵魏腿菀椎牟襟E。

前言

  當(dāng)我們想研究不同sample的某個(gè)變量A之間的差異時(shí)，往往會(huì)因?yàn)槠渌恍┳兞緽對該變量的固有影響，而影響不同sample變量A的比較，這個(gè)時(shí)候需要對sample變量A進(jìn)行標(biāo)準(zhǔn)化之后才能進(jìn)行比較。標(biāo)準(zhǔn)化的方法是對sample 的 A變量和B變量進(jìn)行l(wèi)oess回歸，擬合變量A關(guān)于變量B的函數(shù) f(b)，f(b)則表示在B的影響下A的理論取值，A-f(B)（A對f(b）殘差）就可以去掉B變量對A變量的影響,此時(shí)殘差值就可以作為標(biāo)準(zhǔn)化的A值在不同sample之間進(jìn)行比較。

Loess局部加權(quán)多項(xiàng)式回歸

  LOWESS最初由Cleveland 提出,后又被Cleveland&Devlin及其他許多人發(fā)展。在R中l(wèi)oess 函數(shù)是以lowess函數(shù)為基礎(chǔ)的更復(fù)雜功能更強(qiáng)大的函數(shù)。主要思想為:在數(shù)據(jù)集合的每一點(diǎn)用低維多項(xiàng)式擬合數(shù)據(jù)點(diǎn)的一個(gè)子集,并估計(jì)該點(diǎn)附近自變量數(shù)據(jù)點(diǎn)所對應(yīng)的因變量值,該多項(xiàng)式是用加權(quán)最小二乘法來擬合;離該點(diǎn)越遠(yuǎn),權(quán)重越小,該點(diǎn)的回歸函數(shù)值就是這個(gè)局部多項(xiàng)式來得到,而用于加權(quán)最小二乘回歸的數(shù)據(jù)子集是由最近鄰方法確定。

  最大優(yōu)點(diǎn):不需要事先設(shè)定一個(gè)函數(shù)來對所有數(shù)據(jù)擬合一個(gè)模型。并且可以對同一數(shù)據(jù)進(jìn)行多次不同的擬合，先對某個(gè)變量進(jìn)行擬合，再對另一變量進(jìn)行擬合，以探索數(shù)據(jù)中可能存在的某種關(guān)系，這是普通的回歸擬合無法做到的。

LOESS平滑方法

  1. 以x0為中心確定一個(gè)區(qū)間，區(qū)間的寬度可以靈活掌握。具體來說，區(qū)間的寬度取決于q=fn。其中q是參與局部回歸觀察值的個(gè)數(shù)，f是參加局部回歸觀察值的個(gè)數(shù)占觀察值個(gè)數(shù)的比例,n是觀察值的個(gè)數(shù)。在實(shí)際應(yīng)用中，往往先選定f值，再根據(jù)f和n確定q的取值，一般情況下f的取值在1/3到2/3之間。q與f的取值一般沒有確定的準(zhǔn)則。增大q值或f值，會(huì)導(dǎo)致平滑值平滑程度增加，對于數(shù)據(jù)中前在的細(xì)微變化模式則分辨率低，但噪聲小，而對數(shù)據(jù)中大的變化模式的表現(xiàn)則比較好；小的q值或f值，曲線粗糙，分辨率高，但噪聲大。沒有一個(gè)標(biāo)準(zhǔn)的f值，比較明智的做法是不斷的調(diào)試比較。

  2. 定義區(qū)間內(nèi)所有點(diǎn)的權(quán)數(shù)，權(quán)數(shù)由權(quán)數(shù)函數(shù)來確定，比如立方加權(quán)函數(shù)weight = (1 - (dist/maxdist)^3)^3),dist為距離x的距離，maxdist為區(qū)間內(nèi)距離x的最大距離。任一點(diǎn)(x0,y0)的權(quán)數(shù)是權(quán)數(shù)函數(shù)曲線的高度。權(quán)數(shù)函數(shù)應(yīng)包括以下三個(gè)方面特性：(1)加權(quán)函數(shù)上的點(diǎn)(x0,y0)具有最大權(quán)數(shù)。(2)當(dāng)x離開x0(時(shí)，權(quán)數(shù)逐漸減少。(3)加權(quán)函數(shù)以x0為中心對稱。

  3. 對區(qū)間內(nèi)的散點(diǎn)擬合一條曲線y=f(x)。擬合的直線反映直線關(guān)系，接近x0的點(diǎn)在直線的擬合中起到主要的作用，區(qū)間外的點(diǎn)它們的權(quán)數(shù)為零。

  4. x0的平滑點(diǎn)就是x0在擬合出來的直線上的擬合點(diǎn)(y0,f( x0))。

  5. 對所有的點(diǎn)求出平滑點(diǎn)，將平滑點(diǎn)連接就得到Loess回歸曲線。

R語言代碼

 loess(formula, data, weights, subset, na.action, model = FALSE,  span = 0.75, enp.target, degree = 2,  parametric = FALSE, drop.square = FALSE, normalize = TRUE,  family = c("gaussian", "symmetric"),  method = c("loess", "model.frame"),  control = loess.control(...), ...)

  formula是公式，比如y~x,可以輸入1到4個(gè)變量;

  data是放著變量的數(shù)據(jù)框，如果data為空，則在環(huán)境中尋找;

  na.action指定對NA數(shù)據(jù)的處理，默認(rèn)是getOption("na.action");

  model是否返回模型框；

  span是alpha參數(shù)，可以控制平滑度,相當(dāng)于上面所述的f，對于alpha小于1的時(shí)候，區(qū)間包含alpha的點(diǎn)，加權(quán)函數(shù)為立方加權(quán)，大于1時(shí)，使用所有的點(diǎn)，最大距離為alpha^(1/p)，p 為解釋變量;

  anp.target，定義span的備選方法；

  normalize，對多變量normalize到同一scale；

  family，如果是gaussian則使用最小二乘法，如果是symmetric則使用雙權(quán)函數(shù)進(jìn)行再下降的M估計(jì)；

  method，是適應(yīng)模型或者僅僅提取模型框架；

  control進(jìn)一步更高級的控制，使用loess.control的參數(shù)；

  其它參數(shù)請自己參見manual并且查找資料

loess.control(surface = c("interpolate", "direct"),   statistics = c("approximate", "exact"),   trace.hat = c("exact", "approximate"),   cell = 0.2, iterations = 4, ...)

  surface，擬合表面是從kd數(shù)進(jìn)行插值還是進(jìn)行精確計(jì)算；

  statistics,統(tǒng)計(jì)數(shù)據(jù)是精確計(jì)算還是近似，精確計(jì)算很慢

  trace.hat,要跟蹤的平滑的矩陣精確計(jì)算或近似？建議使用超過1000個(gè)數(shù)據(jù)點(diǎn)逼近，

  cell,如果通過kd樹最大的點(diǎn)進(jìn)行插值的近似。大于cell floor(nspancell)的點(diǎn)被細(xì)分。

  robust fitting使用的迭代次數(shù)。

predict(object, newdata = NULL, se = FALSE, na.action = na.pass, ...)

  object，使用loess擬合出來的對象；

  newdata,可選數(shù)據(jù)框，在里面尋找變量并進(jìn)行預(yù)測；

  se,是否計(jì)算標(biāo)準(zhǔn)誤差；

  對NA值的處理

實(shí)例

  生物數(shù)據(jù)分析中，我們想查看PCR擴(kuò)增出來的擴(kuò)增子的測序深度之間的差異，但不同的擴(kuò)增子的擴(kuò)增效率受到GC含量的影響，因此我們首先應(yīng)該排除掉GC含量對擴(kuò)增子深度的影響。

數(shù)據(jù)

amplicon 測序數(shù)據(jù)，處理后得到的每個(gè)amplicon的深度，每個(gè)amplicon的GC含量，每個(gè)amplicon的長度

先用loess進(jìn)行曲線的擬合

gcCount.loess <- loess(log(RC+0.01)~GC,data=RC_DT,control = loess.control(surface = "direct"),degree=2)

畫出擬合出來的曲線

predictions1<- predict (gcCount.loess,RC_DT$GC)#plot scatter and line plot(RC_DT$GC,log(RC_DT$RC+0.01),cex=0.1,xlab="GC Content",ylab=expression(paste("log(NRC"["lib"],"+0.01)",sep="")))lines(RC_DT$GC,predictions1,col = "red")

取殘差，去除GC含量對深度的影響

#sustract the influence of GCresi <- log(RC_DT$RC+0.01)-predictions1RC_DT$RC <- resisetkey(RC_DT,GC)

此時(shí)RC_DT$RC就是normalize之后的RC

畫圖顯示nomalize之后的RC,并將擬合的loess曲線和normalize之后的數(shù)據(jù)保存

#plot scatter and line using Norm GC dataplot(RC_DT$GC,RC_DT$RC,cex=0.1,xlab="GC Content",ylab=expression("NRC"["GC"]))gcCount.loess <- loess(RC~GC,data=RC_DT,control = loess.control(surface = "direct"),degree=2)save(gcCount.loess,file="/home/ywliao/project/Gengyan/gcCount.loess.Robject")predictions2 <- predict(gcCount.loess,RC_DT$GC)lines(RC_DT$GC,predictions2,col="red")save(RC_DT,file="/home/ywliao/project/Gengyan/RC_DT.Rdata")

當(dāng)然，也想看一下amplicon 長度len 對RC的影響，不過影響不大

全部代碼如下(經(jīng)過修改，可能與上面完全匹配)：

library(data.table)load("/home/ywliao/project/Gengyan/RC_DT.Rdata")RRC_DT <- RC_DT[Type=="WBC" & !is.na(RC),]lst <- list()for (Samp in unique(RC_DT$Sample)){RC_DT <- RRC_DT[Sample==Samp]####loess GC vs RC####gcCount.loess <- loess(log(RC+0.01)~GC,data=RC_DT,control = loess.control(surface = "direct"),degree=2)predictions1<- predict (gcCount.loess,RC_DT$GC)#plot scatter and line #plot(RC_DT$GC,log(RC_DT$RC+0.01),cex=0.1,xlab="GC Content",ylab=expression(paste("log(NRC"["lib"],"+0.01)",sep="")))#lines(RC_DT$GC,predictions1,col = "red")#sustract the influence of GCresi <- log(RC_DT$RC+0.01)-predictions1RC_DT$NRC <- resisetkey(RC_DT,GC)#plot scatter and line using Norm GC data#plot(RC_DT$GC,RC_DT$NRC,cex=0.1,xlab="GC Content",ylab=expression("NRC"["GC"]))gcCount.loess <- loess(NRC~GC,data=RC_DT,control = loess.control(surface = "direct"),degree=2)predictions2 <- predict(gcCount.loess,RC_DT$GC)#lines(RC_DT$GC,predictions2,col="red")lst[[Samp]] <- RC_DT}NRC_DT <- rbindlist(lst)save(RC_DT,file="/home/ywliao/project/Gengyan/NRC_DT.Rdata")####loess len vs RC###setkey(RC_DT,Len)len.loess <- loess(RC_DT$NRC~RC_DT$Len, control = loess.control(surface = "direct"),degree=2)predictions2<- predict (len.loess,RC_DT$Len)#plot scatter and line plot(RC_DT$Len,RC_DT$NRC,cex=0.1,xlab="Length",ylab=expression(paste("log(RC"["GC"],"+0.01)",sep="")))lines(RC_DT$Len,predictions2,col = "red")

總結(jié)

以上就是這篇文章的全部內(nèi)容了，希望本文的內(nèi)容對大家的學(xué)習(xí)或者工作具有一定的參考學(xué)習(xí)價(jià)值，如果有疑問大家可以留言交流，謝謝大家對武林網(wǎng)的支持。