為什麼我一行代碼就可以完成3個R包的RNA-seq差異分析呢

收錄於話題 #RNA 36個

在教師節收到學生提問，刷我B站74小時視頻的時候看到我演示了RNA-seq差異分析只用了一行代碼就完成了3大R包的全部分析，並且輸出了對應的圖表結果，覺得很神奇，但是B站視頻並沒有配套講義和代碼還有測試數據。

首先我一直使用airway數據集做測試

airway數據集這裡我就不多說了，搜索生信技能樹早期教程可以看到很多介紹，使用下面代碼就可以簡單探索。

if(F){
  library(airway)
  data(airway)
  exprSet=assay(airway)
  group_list=colData(airway)[,3]
  save(exprSet,group_list,file = 'airway_exprSet.Rdata')
}

load(file = 'airway_exprSet.Rdata')

if(T){
  colnames(exprSet)
  pheatmap::pheatmap(cor(exprSet))
  group_list
  tmp=data.frame(g=group_list)
  rownames(tmp)=colnames(exprSet)
  
  pheatmap::pheatmap(cor(exprSet),annotation_col = tmp)
  dim(exprSet)
  exprSet=exprSet[apply(exprSet,1, function(x) sum(x>1) > 5),]
  dim(exprSet)

  exprSet=log(edgeR::cpm(exprSet)+1)
  dim(exprSet)
  exprSet=exprSet[names(sort(apply(exprSet, 1,mad),decreasing = T)[1:500]),]
  dim(exprSet)
  M=cor(log2(exprSet+1))
  tmp=data.frame(g=group_list)
  rownames(tmp)=colnames(M)
  pheatmap::pheatmap(M,annotation_col = tmp)
  pheatmap::pheatmap(M,annotation_col = tmp,filename = 'cor.png')

  library(pheatmap)
  pheatmap(scale(cor(log2(exprSet+1))))

}

很明顯可以看到， 組內的樣本的相似性應該是要高於組間的！

而且為了顯示這個規律，我還做了一個統計學技巧展示，當然了，很多人非常的不用心，所以把視頻聽10遍也看不懂，get不到我的點，需要批評！

使用我包裝好的函數即可

可以看到，下面的代碼非常簡潔，因為僅僅是使用了 run_DEG_RNAseq 函數，就根據表達矩陣和分組信息，完成了全部的分析！

rm(list = ls())
options(stringsAsFactors = F)
load(file = 'airway_exprSet.Rdata')
group_list
group_list=relevel(group_list,ref = 'untrt')
source('run_DEG_RNA-seq.R')
run_DEG_RNAseq(exprSet,group_list,
               g1="untrt",g2="trt",
               pro='airway')

這就是大家看視頻後提的問題，為什麼這麼神奇呢？下面的圖表是如何自動出來的呢？

因為這個 run_DEG_RNAseq 函數的代碼非常長，這裡我就不貼在公眾號了哈，大家可以在我的GitHub的GEO項目找到它！

https://github.com/jmzeng1314/GEO/blob/master/airway_RNAseq/run_DEG_RNA-seq.R

GEO傳奇代碼

一不留神，這個GEO項目就成為了點讚數最多的，直接孵化出12篇數據挖掘類SCI文章，至於間接的那些就不計其數了，因為大家都是偷偷的使用，也不告訴我，甚至某些別有用心者還不告訴身邊的人，要一個人獨享這些代碼。

https://github.com/jmzeng1314/GEO

https://github.com/jmzeng1314/GEO/tree/master/airway_RNAseq

既然是多個R包，結果該如何取捨呢？

這個時候是沒有標準答案的，因為每個R包都非常熱門，引用量都是好幾千，你選擇哪個都符合市場規律，不過，我這裡有一個代碼，對3個結果根據閾值篩選交集。

https://github.com/jmzeng1314/GEO/tree/master/airway_RNAseq

差異基因後是不是也可以批量GO/KEGG資料庫注釋呢？

當然是啊，都會寫代碼了，還有什麼是不能為所欲為的呢？

同樣的，代碼也是在GitHub，需要你仔細理解，不過我有一個小小的要求，請不要把我的代碼雪藏，或者刻意隱瞞。

https://github.com/jmzeng1314/GEO/tree/master/airway_RNAseq

值得一提的是這裡面的一行代碼是需要格外注意的哦：

group_list=relevel(group_list,ref = 'untrt')

本來到這裡應該是要貼上我們全國巡講的宣傳，但是中秋節臨時加開的廣州特別班，我們沒想過會有很多人報名，畢竟只有區區十天不到的報名時間，但是很快就滿20人了，而且非常多迫切想學習的小夥伴找到我們，即使拒絕後仍然是發生如下的對話：

因為加人就要換大一點會議室，成本就增加幾千塊錢，所以只能是看有沒有3個以上的人報名，至少把成本cover掉！

大家仍然是可以嘗試報名，廣州今年就這一場了，還有很多其他城市的粉絲嗷嗷待哺者等著我們！

號外：中秋節廣州3天入門課程報名馬上截止：（中秋節一起來學習！）全國巡講第16站-廣州（生信入門課加量不加價）