TGCA數據的標準化以及差別分析--轉載

轉載果子學生信  https://mp.weixin.qq.com/s/Ph1O6V5RkxkyrKpVmB5ODA

前面咱們從GDC下載了TCGA腫瘤數據庫的數據,也可以把GDC下載的多個TCGA文件批量讀入R

今天咱們講一下TCGA數據的標準化，以及差別分析，獲得了標準化後的數據，咱們就能夠按照之前的帖子，作一系列操做

Y叔推薦的這個圖有毒！

圖有毒系列之2

多個基因在多亞組疾病中的展現

在獲得了差別分析的結果後，咱們能夠完成熱圖，火山圖，GO分析，KEGG分析，GSEA分析，就跟這個帖子中的同樣。
來完成你的生信做業，這是最有誠意的GEO數據庫教程

下面開始今天的教程：
首先加載上一次課得到的數據；

### 加載數據
load("expr_df.Rdata")

如今的數據是這個樣子的

處理前

去掉ensemble ID的點號

library(tidyr)
expr_df_nopoint <- expr_df %>% 
  tidyr::separate(gene_id,into = c("gene_id"),sep="\\.") 

如今的數據是這個樣子的

處理後

去掉點號，是爲了用gtf文件。
gtf文件的獲取和做用在這裏
GTF文件有什麼用啊？別的不談，最起碼能提lncRNA

加載gtf文件,這是目前咱們能接觸的最大文件，有260萬行。

load(file = "gtf_df.Rda")

提取mRNA

mRNA_exprSet <- gtf_df %>% 
  dplyr::filter(type=="gene",gene_biotype=="protein_coding") %>% #篩選gene,和編碼指標
  dplyr::select(c(gene_name,gene_id,gene_biotype)) %>% 
  dplyr::inner_join(expr_df_nopoint,by ="gene_id") %>% 
  tidyr::unite(gene_id,gene_name,gene_id,gene_biotype,sep = " | ")

最終獲得19668行，這是編碼基因的個數，如今的數據是這個樣子的

編碼RNA

提取lncRNA

這裏頗有爭議，而個人理由是，即便是編碼基因，也會出現非編碼轉錄本，而長鏈非編碼RNA，指的是轉錄本，因此不能用gene的編碼與否來界定

ncRNA <- c("sense_overlapping","lincRNA","3prime_overlapping_ncRNA",
           "processed_transcript","sense_intronic",
           "bidirectional_promoter_lncRNA","non_coding",
           "antisense_RNA")

LncRNA_exprSet <- gtf_df %>% 
  dplyr::filter(type=="transcript",transcript_biotype %in% ncRNA) %>% #注意這裏是transcript_biotype
  dplyr::select(c(gene_name,gene_id,transcript_biotype)) %>% 
  dplyr::distinct() %>% #刪除多餘行
  dplyr::inner_join(expr_df_nopoint,by ="gene_id") %>% 
  tidyr::unite(gene_id,gene_name,gene_id,transcript_biotype,sep = " | ")

最終獲得25530個非編碼轉錄本，數據是這個樣子的

非編碼RNA

數據標準化

標準化和差別分析都是用Deseq2這個包來完成，首先要構建dds對象，構建這個對象須要兩個文件，第一是輸入數據，咱們已經有了，第二個是分組文件metadata，他至少由兩列構成，一列是樣本名稱，一列是分組信息。

首先把樣本名稱變成數據框格式

metadata <- data.frame(TCGA_id =colnames(expr_df)[-1])

分組信息包含在TCAG_id的第14,15字符頗有用，他指示了樣本是癌症仍是癌旁或者是轉移 病竈

官網解釋以下,01-09是癌症，10-19是正常，20-29是癌旁

Tumor types range from 01 - 09, normal types from 10 - 19 and control samples from 20 - 29

TCGA barcode的詳細信息以下：

同時咱們要注意，即便是腫瘤組織，01-09意義各不相同，好比，01表明原發竈，02表明轉移竈，詳細信息以下：

咱們用table這個函數統計一下腦膠質瘤GBM樣本的分類

table(substring(metadata$TCGA_id,14,15))

有154個是原發竈，有13個是轉移竈，很奇怪是吧，沒有癌旁。可是這個是能理解的，人的大腦正常組織是有用的，不一樣於肝臟這類奇怪多一塊少一塊無所謂，切取大腦正常組織是沒有倫理的。實際上TCGA裏面還有一部分腫瘤是沒有癌旁的，好比，淋巴瘤。

這一部分沒有正常對照的腫瘤如何進行差別分析呢，一種方法是，使用GTEx數據庫中的正常組織，這個咱們留一個坑，之後再講。

可是，今天咱們的活仍是要作，咱們就用復發和非復發來區分便可。

sample <- ifelse(substring(metadata$TCGA_id,14,15)=="01","cancer","recur")
## 這裏的factor是爲了dds的須要
metadata$sample <- as.factor(sample)

此時metadata是這個樣子的

構建dds的兩個文件所有準備好，咱們開始下一步

mRNA標準化

這一步是爲了代碼複用，把counts文件統一命名

mycounts <- mRNA_exprSet

構建dds對象，若是mycounts中的TCGA_id是行名，tidy這個參數設置爲FASLE

dds <-DESeqDataSetFromMatrix(countData=mycounts, 
                             colData=metadata, 
                             design=~sample,
                             tidy=TRUE)

Deseq2分析，這裏面有不少步驟都本身運行了，這一步十分耗時，取決於樣本數以及電腦內存大小，個人16g內存電腦運行5分鐘，而個人學員們有的人要運行20個小時。甚至，若是，你分析的是乳腺癌，1000多個樣本，小電腦根本過不去，此時，你能夠考慮升級一下裝備。

dds <- DESeq(dds)

這個數據很重要，並且有些人得到也不容易，因此，須要保存一下,方便之後使用。

save(dds,file="mRNA_exprSet_dds_sample.Rdata")

vst標準化,這一步跟上一步同樣，速度取決於樣本量和電腦

vsd <- vst(dds, blind = FALSE)

爲何選擇vst呢？看這個
轉錄組的高級分析前該如何標準化數據？
Deseq2標準化的原理是什麼，youtube上的StatQuest小哥視頻說的特別好，能夠看看這個帖子
DESeq2的標準化方法

這時候，Deseq2還內置了主成分分析來看一下樣本分佈

plotPCA(vsd, "sample")

從圖上咱們能夠看出，原發竈和轉移竈，並不能完美分開，生物學意義就是，轉移竈不是新的類型的腫瘤，他實際上仍是腦膠質瘤，後續可能發生的結果是，下游額差別分析接結果很差，可能的解決方法是，找出配對的原發竈和轉移竈來分析。咱們看結果來講話。

獲取標準化後的數據,這一步還會自動過濾掉不符合規定的基因，這時候，數據明顯被標準化了

mRNA_exprSet_vst <- as.data.frame(assay(vsd))

標準化以後

保存一下這個數據，調整一下格式，就能夠用於本文開頭說的那一系列操做。

save(ncRNA_exprSet_vst,file = "ncRNA_exprSet_vst.Rda")

差別分析

這裏用到前面保存的dds,使用results函數提取

res <- results(dds, tidy=TRUE) 

咱們看到這個數據，有foldchange值，有pvlaue，那麼篩選差別基因，熱圖，火山圖，GO，KEGG分析，GSEA分析就瓜熟蒂落啦。