R

相比于網(wǎng)頁工具，使用編程語言處理科研數(shù)據(jù)的一大優(yōu)勢，在于高度的定制化，以及批量處理數(shù)據(jù)的快捷性和高效性

目錄
批量處理——for循環(huán)批量計(jì)算組間差異
批量處理——apply批量計(jì)算組間差異
批量處理——for循環(huán)畫圖
批量處理——for循環(huán)遷移文件

對于批量處理數(shù)據(jù)的方法，之前使用for循環(huán)和apply語句進(jìn)行處理過，但是不夠系統(tǒng)，學(xué)習(xí)果子生信課程后有一個清晰的認(rèn)識，寫下來，一是可以調(diào)用方便，二是自己寫過之后，才能算是完全掌握。當(dāng)然一切以解決問題為主，不陷于技術(shù)深究。

批量計(jì)算基因和基因之間的相關(guān)性，也是一項(xiàng)很好的應(yīng)用。

場景

計(jì)算RNA-seq得到的基因與DDR1基因相關(guān)系數(shù)，并挑選出顯著相關(guān)基因（p < 0.05）

數(shù)據(jù)準(zhǔn)備

得到的數(shù)據(jù)是TCGA結(jié)直腸癌的數(shù)據(jù)

為上圖的清潔數(shù)據(jù)

for循環(huán)計(jì)算

首先要明確單次操作的情況

cor_data <- data.frame()

yourgene <- 'DDR1'

genelist <-  setdiff(colnames(COAD_data),c("sample","TCGA_id","type","subtype",yourgene))


dd <- cor.test(COAD_data[,yourgene], COAD_data[,genelist[1]], method = 'spearman')
dd$statistic
dd$p.value
cor_data[genelist[1],1] = dd$statistic
cor_data[genelist[1],2] = dd$p.value
colnames(cor_data) = c('statistic', 'p.value')

批量處理時，針對for循環(huán)就是提前做好一個容器，然后把變量替換成一般表達(dá)

cor_data <- data.frame()

yourgene <- 'DDR1'

genelist <-  setdiff(colnames(COAD_data),c("sample","TCGA_id","type","subtype",yourgene))

system.time(
  for (gene in genelist) {
  dd <- cor.test(COAD_data[,yourgene], COAD_data[,gene], method = 'spearman')
  statistic = dd$statistic
  p.value = dd$p.value
  cor_data[gene,1] = statistic
  cor_data[gene,2] = p.value 
}
)
# 用戶  系統(tǒng)  流逝 
# 34.62  0.66 36.62 

colnames(cor_data) = c('statistic', 'p.value')

lapply計(jì)算

yourgene <- 'DDR1'
genelist <-  setdiff(colnames(COAD_data),c("sample","TCGA_id","type","subtype",yourgene))

# 寫一個函數(shù)
mycor <- function(x){
  dd = cor.test(COAD_data[,yourgene], COAD_data[,x], method = 'spearman')
  data.frame(gene1 = yourgene, gene2 = x, R = dd$estimate, p.value = dd$p.value)
}

# 測試一下
mycor(genelist[1])

# 批量
lapplylist <- lapply(genelist, mycor)

# do.call
cor_data <- do.call(rbind, lapplylist)

# 整理成一個函數(shù)
cor_data <- do.call(rbind, lapply(genelist, function(x){
  dd = cor.test(COAD_data[,yourgene], COAD_data[,x], method = 'spearman')
  data.frame(gene1 = yourgene, gene2 = x, R = dd$estimate, p.value = dd$p.value)
}))

再測試一下時間

system.time(
  cor_data <- do.call(rbind, lapply(genelist, function(x){
    dd = cor.test(COAD_data[,yourgene], COAD_data[,x], method = 'spearman')
    data.frame(gene1 = yourgene, gene2 = x, R = dd$estimate, p.value = dd$p.value)
  }))
)

# 用戶  系統(tǒng)  流逝 
# 31.29  0.25 32.15 

結(jié)果

兩種計(jì)算方法的得到結(jié)果是一致的

關(guān)于篩選和排序

那就只是filter函數(shù)和arrange函數(shù)的操作，可參閱R for Data Science（筆記） ---數(shù)據(jù)變換（filter使用） ，R for Data Science（筆記） ---數(shù)據(jù)變換（行排序）操作。

后記：本篇是寫了如何批量求與某一個基因的基因的相關(guān)系數(shù)。后續(xù)設(shè)計(jì)兩個基因在不同腫瘤的相關(guān)性。