|
致 歉 信 本文早在2020年3月13日就已經(jīng)在本公眾號發(fā)布，截止到2020年5月4日，已經(jīng)有60人付費，由于本人的操作不當(dāng)，意外把文章刪除了。對已經(jīng)付費的粉絲們，深感抱歉。其實，我早已將視頻上傳到B站(未提供代碼)： https://www.bilibili.com/video/BV117411Z7Cq 相信付費的您已經(jīng)拿到代碼和相關(guān)文件，如果有什么問題可以聯(lián)系我。 郵箱：bioinfocloud@aliyun.com DoubleHelix 2020.05.04   關(guān)于TCGA數(shù)據(jù)庫的教程，前期我們已經(jīng)推出了一些文章： 【2】R語言TCGA-Assembler包下載TCGA數(shù)據(jù)前面的這些文章雖然介紹了TCGA數(shù)據(jù)庫的差異表達(dá)，但部分用了perl腳本，如果不懂perl的同學(xué)，數(shù)據(jù)庫更新，數(shù)據(jù)格式可能會變換，運行腳本就得不到想要的結(jié)果，所以這里我們就只用R語言進行系統(tǒng)性的講解，包括數(shù)據(jù)下載、整理、融合、基因ID轉(zhuǎn)換以及表達(dá)差異分析，最后通過火山圖進行可視化。所以你需要有一定的R基礎(chǔ)，能看的懂代碼，能改一些參數(shù)。 關(guān)于差異表達(dá)分析我們利用DESeq2和EdgeR包，其實在我們前面基因芯片數(shù)據(jù)挖掘序列文章中都已介紹： 基因芯片數(shù)據(jù)分析（一）：芯片數(shù)據(jù)初探 基因芯片數(shù)據(jù)分析（二）：讀取芯片數(shù)據(jù) 基因芯片數(shù)據(jù)分析（三）：數(shù)據(jù)質(zhì)控 基因芯片數(shù)據(jù)分析（四）：獲取差異表達(dá)基因 基因芯片數(shù)據(jù)分析（五）：edgeR包的基本原理 基因芯片數(shù)據(jù)分析（六）：DESeq2包的基本原理 基因芯片數(shù)據(jù)分析（七）：edgeR差異分析實戰(zhàn)案例 基因芯片數(shù)據(jù)分析（八）：DESeq2差異分析實戰(zhàn)案例 可以先通過五~八這4篇文章了解這2個包的原理和使用教程。其實，你只要能從TCGA數(shù)據(jù)庫下載的數(shù)據(jù)整理得到表達(dá)矩陣，參照七、八這2篇文章就可以得到差異表達(dá)結(jié)果。下面我們步入正題........... 一.數(shù)據(jù)下載 數(shù)據(jù)下載有3中方式，官網(wǎng)在線下載；官放下載工具下載；R語言包下載，比如：TCGAbiolinks。TCGA-Assembler等。我們推薦使用TCGAbiolinks，個人覺得這是挖掘TCGA數(shù)據(jù)比較好用的包。 ########################## 這里下載的是Counts數(shù)據(jù)，不是FPKM數(shù)據(jù)####################setwd('.')##包的安裝if (!requireNamespace('BiocManager', quietly = TRUE)) install.packages('BiocManager')BiocManager::install('TCGAbiolinks')# 加載相應(yīng)的包，可能會需要其他包，提示錯誤就安裝缺少的包。# 因為每個人已經(jīng)安裝的包不一樣。library(TCGAbiolinks)# 請求數(shù)據(jù)。query <- GDCquery(project = 'TCGA-LUAD', data.category = 'Transcriptome Profiling', data.type = 'Gene Expression Quantification', workflow.type = 'HTSeq - Counts')# 從query中獲取結(jié)果表，它可以選擇帶有cols參數(shù)的列，并使用rows參數(shù)返回若干行。# 594個barcode samplesDown <- getResults(query,cols=c('cases'))# 533個barcodedataSmTP <- TCGAquery_SampleTypes(barcode = samplesDown, typesample = 'TP')# 59個barcodedataSmNT <- TCGAquery_SampleTypes(barcode = samplesDown, typesample = 'NT')# 59個正常組織和533個腫瘤組織樣本作為研究dataSmTP_short <- dataSmTP[1:533]dataSmNT_short <- dataSmNT[1:59]# 根據(jù)前面的篩選，再次請求數(shù)據(jù)queryDown <- GDCquery(project = 'TCGA-LUAD', data.category = 'Transcriptome Profiling', data.type = 'Gene Expression Quantification', workflow.type = 'HTSeq - Counts', barcode = c(dataSmTP_short, dataSmNT_short)) 網(wǎng)頁篩選后添加購物車就可以直接下載了。最后一種下載方式是通過官方工具Data Transfer Tool： https://gdc./access-data/gdc-data-transfer-tool。 關(guān)于網(wǎng)頁下載需要下載幾個文件。  如果直接使用TCGAbiolinks包下載，并用該包進行數(shù)據(jù)分析，就沒有這么麻煩。我們這里主要是講解通過網(wǎng)頁下載的數(shù)據(jù)，直接整理，分析。讓大家知道更加詳細(xì)的流程。只是下載數(shù)據(jù)的話，上面三種方式都可以，不過現(xiàn)在TCGAbiolinks包下載速度可慢了，而且這個包現(xiàn)在好像因為R版本的更新，部分函數(shù)好像會出錯。網(wǎng)頁版下載數(shù)據(jù)是把篩選好的樣本數(shù)據(jù)打包下載，文件太大，網(wǎng)速不好可能會中斷。但一般也沒什么問題。 二. 數(shù)據(jù)的整理 我們網(wǎng)頁上下載的數(shù)據(jù)是一個壓縮包，需要解壓。  解壓后我們會看到很多文件夾：   解壓以后是一個txt的文本文件，打開文本文件就是該文件對于病人的RNA-Seq數(shù)據(jù)，第一列是Ensembl ID ，第二列就是Counts數(shù)。   三 .ID轉(zhuǎn)換 得到表達(dá)矩陣后，Ensembl ID我們看不懂，需要轉(zhuǎn)換成我們能看的懂的基因名稱，比如TP53。我們就需要通過gtf文件進行基因轉(zhuǎn)換，人的基因注釋文件下載地址可通過gencode：https://www./human/ 下載，NCBI和Ensembl 也可以下載。  但是我們需要注意的是多個Ensembl ID可能對于一個基因。所以我們需要去重,去重可以使用下面函數(shù)，取平均值。 當(dāng)然也可以選擇其中一個，刪掉重復(fù)的數(shù)據(jù)，建議取平均值。 四，表達(dá)差異分析 得到這樣的表達(dá)矩陣，我們就可以進行差異表達(dá)分析了。  DESeq2和EdgeR包進行表達(dá)差異分析需要原始的Counts表達(dá)矩陣文件，我們前面下載的也是HTSeq - Counts數(shù)據(jù)，如果是HTSeq - FPKM的數(shù)據(jù)可以用limma包直接分析，這2類數(shù)據(jù)在下載和處理得到矩陣文件過程是一樣的。后面也會附代碼。 進行差異分析以后，我們繪制火山圖進行可視化。 |
|
|
