測(cè)序上游分析系列:
mRNA-seq轉(zhuǎn)錄組二代測(cè)序從raw reads到表達(dá)矩陣:上中游分析pipeline
miRNA-seq小RNA高通量測(cè)序pipeline:從raw reads,鑒定已知miRNA-預(yù)測(cè)新miRNA后专,到表達(dá)矩陣【一】
miRNA-seq小RNA高通量測(cè)序pipeline:從raw reads划鸽,鑒定已知miRNA-預(yù)測(cè)新miRNA,到表達(dá)矩陣【二】
其他文章系列:ggplot2作圖篇:
(1)基于ggplot2的RNA-seq轉(zhuǎn)錄組可視化:總述和分文目錄
(2)測(cè)序結(jié)果概覽:基因表達(dá)量rank瀑布圖戚哎,高密度表達(dá)相關(guān)性散點(diǎn)圖裸诽,注釋特定基因及errorbar的表達(dá)相關(guān)性散點(diǎn)圖繪制
(3)單/多個(gè)基因在組間同時(shí)展示的多種選擇:分組小提琴圖、分組/分面柱狀圖丈冬、單基因蜂群點(diǎn)圖拼圖的繪制
(4)配對(duì)樣本基因表達(dá)趨勢(shì):優(yōu)化前后的散點(diǎn)連線圖+拼圖繪制
TCGA是The Cancer Genome Atlas的簡(jiǎn)稱蓄氧,由美國(guó)NIH管理,是目前出現(xiàn)頻率最高,最權(quán)威的腫瘤公共數(shù)據(jù)庫(kù)。其中收錄了人體多個(gè)部位腫瘤的高通量測(cè)序/基因芯片mRNA局扶、miRNA表達(dá)譜數(shù)據(jù)畴蒲、基因組變異、甲基化等數(shù)據(jù),以及相應(yīng)樣本的人口學(xué)特征及臨床數(shù)據(jù)(包括暴露情況、預(yù)后、腫瘤分期骚腥、治療手段)。為腫瘤研究的數(shù)據(jù)挖掘和課題前期分析提供了豐富的資源。
目前課題非常常見(jiàn)的手(tao)段(lu)就是使用公共數(shù)據(jù)庫(kù)挖掘腫瘤高通量轉(zhuǎn)錄組數(shù)據(jù),發(fā)現(xiàn)對(duì)抑癌/促癌具有潛在價(jià)值的基因或基因模塊,隨后結(jié)合功能實(shí)驗(yàn)進(jìn)行驗(yàn)證和深入機(jī)制的探討。所以對(duì)沒(méi)有條件自己送臨床樣本/動(dòng)物實(shí)驗(yàn)測(cè)序的人來(lái)說(shuō),TCGA是妥妥的寶藏~
就從TCGA RNA-seq表達(dá)譜數(shù)據(jù)下載開(kāi)始吧泻轰。寫(xiě)下這個(gè)備忘錄泳挥,以免以后又出現(xiàn)代碼不用就忘記的尷尬筑煮。
TCGA官方網(wǎng)站為:https://portal.gdc.cancer.gov
1.點(diǎn)擊Repository框。
- 接下來(lái)就進(jìn)入數(shù)據(jù)類型選擇的頁(yè)面桑谍。
選擇下載轉(zhuǎn)錄組二代測(cè)序數(shù)據(jù)(芯片已經(jīng)基本成為過(guò)去式了)雹仿,在左側(cè)的File中需要選擇Data Category > Transcriptome Profiling, Data Type > Gene Expression Quantification, Experimental Strategy > RNA-seq.
Workflow type >: 關(guān)于檢測(cè)gene expression的RNA-seq的數(shù)據(jù)類型邑商,TCGA提供已經(jīng)將reads比對(duì)好并注釋后的Counts (未經(jīng)標(biāo)準(zhǔn)化的aligned reads數(shù)), FPKM(已標(biāo)準(zhǔn)化fragments per kilobase per million), FPKM-UQ(FPKM-upper quartile)下載含鳞。
如何選擇:(1)如果需要做差異表達(dá)分析熔吗,則需下載HTSeq-Counts咨堤,數(shù)據(jù)整合后使用R的EdgeR或DESeq2包分析凸克,它們均需輸入未經(jīng)標(biāo)準(zhǔn)化的原始counts matrix舆逃。(2)如需做WGCNA模組分析览祖,推薦下載HTSeq-Counts苔咪,經(jīng)DESeq2包標(biāo)準(zhǔn)化后輸入分析。(3)如需做CIBERSORT數(shù)字免疫細(xì)胞分析舔清,需下載HTSeq-FPKM作為輸入數(shù)據(jù)臼婆。(4)...
這里以FPKM數(shù)據(jù)下載為例,但Counts數(shù)據(jù)的流程也完全一樣!
- 選擇需要分析的Case即腫瘤類型和臨床樣本特征。
將左上角的File調(diào)到Case骄瓣,出現(xiàn)Primary site, Program, Gender, Age等信息需要勾選停巷。需根據(jù)課題需要進(jìn)行相關(guān)人口學(xué)特征的勾選。
注:為避免后續(xù)數(shù)據(jù)處理中ID轉(zhuǎn)換麻煩扒磁,Program需選擇TCGA(具有統(tǒng)一的ID特征,可批量處理),project和disease type根據(jù)需要的疾病類型勾選系吭。如需做生存分析敬扛,Vital status不要勾選某一特定生存狀態(tài)侮邀。
本次以TCGA項(xiàng)目中特定人種 (white) 的肺鱗癌 (LUSC, Lung squamous cell carcinoma) 為例坏怪。勾選條件如下:
- 添加數(shù)據(jù)到購(gòu)物車(chē)(Cart)。
此時(shí)樣本類型和特征都已經(jīng)選好绊茧。點(diǎn)擊頁(yè)面中上的Add all files to cart铝宵,然后點(diǎn)擊右上角Cart進(jìn)入購(gòu)物車(chē)下載。Cart上顯示已經(jīng)納入了392個(gè)樣本等待下載按傅。
- 下載樣本的測(cè)序數(shù)據(jù)捉超,臨床數(shù)據(jù)和樣本ID信息胧卤。
點(diǎn)擊Download > Manifest唯绍,該文件隨后將作為下載器的輸入信息文件下載測(cè)序數(shù)據(jù)。(不推薦通過(guò)下載Cart直接獲取測(cè)序數(shù)據(jù)枝誊,因?yàn)閿?shù)據(jù)庫(kù)網(wǎng)絡(luò)速度慢况芒,遇到斷網(wǎng)或無(wú)響應(yīng)就不能斷點(diǎn)續(xù)傳了)。然后同時(shí)需點(diǎn)擊下載Metadata以及Clinical > json。
注:第一次使用TCGA绝骚,需要點(diǎn)擊下載右上角的GDC Data Transfer Tool耐版,根據(jù)自己的平臺(tái)系統(tǒng)選擇合適的zip包下載。解壓后放入用于儲(chǔ)存TCGA下載數(shù)據(jù)的新建文件夾中压汪。
6.測(cè)序FPKM表達(dá)矩陣下載。
注:本例使用蘋(píng)果macbook OS系統(tǒng)操作止剖。我沒(méi)有windows電腦腺阳,所以windows的操作我搜一搜,再補(bǔ)上穿香!
打開(kāi)terminal終端亭引,把目錄轉(zhuǎn)換到儲(chǔ)存上述文件所在的文件夾(根據(jù)自己的文件夾路徑進(jìn)行更改),并以下載的manifest文件作為輸入皮获,使用gdc編譯器下載數(shù)據(jù):
cd /Users/zx/Desktop/TCGA/TCGA-LUSC-white/FPKM/tutorial
./gdc-client download -m gdc_manifest_20200308_044238.txt
接下來(lái)就是等待...
Tips:國(guó)內(nèi)最好在白天下載焙蚓,米國(guó)人睡了,傳輸速度會(huì)快一丟丟...但總體來(lái)說(shuō)還是非常感人洒宝。而且建議掛學(xué)校的學(xué)術(shù)梯子购公,有條件掛米國(guó)學(xué)校的梯子是最好了。
僅作參考待德,392個(gè)樣本君丁,國(guó)內(nèi)下午,我用國(guó)內(nèi)普通網(wǎng)絡(luò)下載了3個(gè)小時(shí)将宪;掛了Yale的梯子下載同一批數(shù)據(jù)绘闷,只用了15分鐘。:)))
7.數(shù)據(jù)全部下載完成:
數(shù)據(jù)下載完成后terminal報(bào)告success. 文件夾中會(huì)出現(xiàn)392個(gè)名字神似亂碼(然而不是)的新文件夾较坛,每個(gè)新文件夾中都包含有.FPKM.txt.gz(一種壓縮包形式)印蔗,這是我們需要的表達(dá)矩陣。我們隨后使用R將每個(gè)單樣本的FPKM表達(dá)矩陣進(jìn)行提取丑勤,整合华嘹,并與臨床數(shù)據(jù)進(jìn)行匹配。
接下來(lái)的部分使用Rstudio完成:
需要R包:dplyr, R.utils, tidyr, jsonlite, rtracklayer法竞,第一次使用需要先使用install.packages()下載安裝耙厚。
install.packages('dplyr')
install.packages('R.utils')
install.packages('tidyr')
install.packages('jsonlite')
#rtracklayer is based on bioconductor, so we need to change the way of installation.
source("https://bioconductor.org/biocLite.R")
biocLite("rtracklayer")
- 設(shè)置自定義工作路徑,加載數(shù)據(jù)岔霸,解壓薛躬,并將解壓后.txt表達(dá)矩陣統(tǒng)一放入新文件夾:
TIPS:gunzip函數(shù)是會(huì)默認(rèn)刪除壓縮包文件,僅保留解壓文件的呆细。數(shù)據(jù)文件下載非常不容易型宝,所以衷心建議各位在批量解壓前先整體備份一盤(pán)...如果因?yàn)槠ヅ鋚attern或者各種問(wèn)題導(dǎo)致代碼沒(méi)有循環(huán)完成,操作就比較麻煩了!建議找出問(wèn)題后直接用備份文件重來(lái)即可~
library('R.utils')
setwd('/Users/zx/Desktop/TCGA/TCGA-LUSC-white/FPKM/tutorial')
#to make sure this directory will be in the first place when performing the for loop.
dir.create('00000_extracted_FPKM')
#The char length of the name of each downloaded expr matrix directory is 36, which is different from other files downloaded from the TCGA website.
dir_FPKM <- dir()[nchar(dir()) == 36]
#decompress the separate files and integrate them into a expression matrix.
for (data in dir_FPKM) {
filename <- list.files(data, pattern = ".*FPKM")
R.utils::gunzip(paste0(data,"/",filename))
file_extracted <- list.files(data, pattern = ".*FPKM.txt$")
file.copy(paste0(data,"/",file_extracted),"00000_extracted_FPKM")
}
- 根據(jù)每個(gè)單獨(dú)表達(dá)矩陣文件的gene_id趴酣,將所有樣本合并為一個(gè)總表達(dá)矩陣梨树。
每個(gè)樣本的FPKM.txt由兩列組成,第一列為各基因的ENSEMBL GENE ID岖寞,第二列為FPKM表達(dá)值抡四,不含標(biāo)題行。
library('dplyr')
setwd('/Users/zx/Desktop/TCGA/TCGA-LUSC-white/FPKM/tutorial/00000_extracted_FPKM')
file_extracted_FPKM <- list.files()
first_file <- read.table(file_extracted_FPKM[1], header = F, sep = '\t', stringsAsFactors = F)
names(first_file) <- c('gene_id', substr(file_extracted_FPKM[1], 1, nchar(file_extracted_FPKM[1])-9))
for (extracted_FPKM in file_extracted_FPKM[2:length(file_extracted_FPKM)]) {
file_appended <- read.table(extracted_FPKM, header = F, sep = '\t', stringsAsFactors = F)
names(file_appended) <- c('gene_id', substr(extracted_FPKM, 1, nchar(file_extracted_FPKM[1])-9)) #delete the suffix .FPKM.txt which is made of 9 chars and just maintain the original file names of each sample for following ID conversion.
first_file <- dplyr::inner_join(first_file, file_appended, by='gene_id')
}
expr_FPKM <- first_file
- 暫停點(diǎn)(Optional):這一步不是必須仗谆。不少小伙伴反映后面expr_FPKM報(bào)錯(cuò)床嫌,然后發(fā)現(xiàn)這里文件名的處理可能有問(wèn)題。因?yàn)榇蠹蚁螺d的文件都是不同的胸私,需要靈活應(yīng)變處理輸入時(shí)文件名變更的情況厌处。
沒(méi)有信心解決的朋友可以直接跳過(guò)這一步。
此時(shí)已獲得最原始的總體FPKM表達(dá)矩陣岁疼,第一列為各基因的ENSEMBL GENE ID阔涉,第一行為樣本名稱(為txt文件中的file_names),可將其輸出備用:
#for original data backup:
setwd('/Users/zx/Desktop/TCGA/TCGA-LUSC-white/FPKM/tutorial/00000_extracted_FPKM')
write.csv(expr_FPKM, 'expr_FPKM.csv', quote = F, row.names = F)
若需將原始表達(dá)矩陣重新載入:
注:由于由36個(gè)字符組成的file_names很多以數(shù)字開(kāi)頭捷绒,而R不允許這樣的命名方式瑰排,所以載入后部分樣本名稱前會(huì)自動(dòng)加一個(gè)字符X。所以需要去除暖侨。這里合并操作:
#to load the expr_FKPM data exported before and remove the char 'X' automatically added in front of colnames starting with number:
setwd('/Users/zx/Desktop/TCGA/TCGA-LUSC-white/FPKM/tutorial/00000_extracted_FPKM')
expr_FPKM <- read.csv('expr_FPKM.csv', header = T, stringsAsFactors = F)
for_replacement_indices <- grep(names(expr_FPKM), pattern = '^X')
names(expr_FPKM)[for_replacement_indices] <- substr(names(expr_FPKM)[for_replacement_indices], 2, 37)
- 從網(wǎng)站下載的metadata文件中獲得各樣本的TCGA_IDs椭住,并將合并表達(dá)矩陣的樣本名轉(zhuǎn)換為完整的TCGA ID (eg: TCGA-33-4538-01A-01R-1201-07)。
#get the TCGA IDs of each sample.
library('jsonlite')
setwd('/Users/zx/Desktop/TCGA/TCGA-LUSC-white/FPKM/tutorial')
metadata <- jsonlite::fromJSON(txt = "metadata.cart.2020-03-08.json")
metadata_id <- dplyr::select(.data = metadata, c(file_name, associated_entities))
file_name <- c()
TCGA_ID <- c()
sub_TCGA_ID_func <- function(x){
x$entity_submitter_id
}
ID_conversion_table <- data.frame('file_names'=substr(metadata_id$file_name, 1, 36),
'TCGA_IDs'=as.character(lapply(metadata_id$associated_entities, FUN = sub_TCGA_ID_func)),
stringsAsFactors = F)
rownames(ID_conversion_table) <- ID_conversion_table$file_names
colnames(expr_FPKM) <- append('gene_id', ID_conversion_table[colnames(expr_FPKM)[2:length(expr_FPKM)],'TCGA_IDs'])
- 進(jìn)一步處理表達(dá)矩陣字逗,將第一列的ENSEMBL GENE ID (eg. ENSG00000273842.1)轉(zhuǎn)化為HUGO GENE SYMBOL (eg. DVL1)京郑,便于認(rèn)識(shí)和處理。并去除表達(dá)矩陣中可能存在的重復(fù)基因葫掉。最終生成的表達(dá)矩陣前三列為Gene symbol, gene id和bio_type.
TIPS:ENSEMBL數(shù)據(jù)庫(kù)下載注釋GTF文件速度真的很一言難盡些举。作為參考,掛了Yale的梯子俭厚,45mb左右的文件下載了2個(gè)小時(shí)户魏。建議不要著急,要微笑:)))
注:ENSEMBL ID與GENE SYMBOL轉(zhuǎn)換時(shí)不含其中的小數(shù)點(diǎn)和之后的數(shù)字挪挤,因此需要將其先去除后再使用GTF注釋文件轉(zhuǎn)換叼丑。
library('tidyr')
library('rtracklayer')
#trim the gene ID by decimals and translate them into gene symbols
expr_FPKM <- tidyr::separate(expr_FPKM, gene_id, into = c('gene_id', 'junk'), sep='\\.')
expr_FPKM <- expr_FPKM[,-2]
#download the annotation file containing ENS IDs, gene symbols, and biotypes.
download.file('ftp://ftp.ensembl.org/pub/release-99/gtf/homo_sapiens/Homo_sapiens.GRCh38.99.chr.gtf.gz',
'Homo_sapiens.GRCh38.99.chr.gtf.gz')
R.utils::gunzip('Homo_sapiens.GRCh38.99.chr.gtf.gz')
gtf1 <- rtracklayer::import('Homo_sapiens.GRCh38.99.chr.gtf')
gtf_df <- as.data.frame(gtf1)
#Only select the protein coding genes.
mRNA_exprSet <- gtf_df %>%
dplyr::filter(type=="gene",gene_biotype=="protein_coding") %>%
dplyr::select(c(gene_name,gene_id,gene_biotype)) %>%
dplyr::inner_join(expr_FPKM,by ="gene_id")
mRNA_exprSet <- mRNA_exprSet[!duplicated(mRNA_exprSet$gene_name),]
- 載入臨床特征數(shù)據(jù),并進(jìn)行其中缺失數(shù)據(jù)的清洗扛门。本例中需要保存的臨床數(shù)據(jù)包括病人的性別(gender)鸠信,腫瘤分期(tumor stage),生存時(shí)間(overall survival)尖飞,以及生存狀態(tài)(原vital status症副,數(shù)據(jù)清洗后改為censoring status)。
TIPS:之前下載的392個(gè)數(shù)據(jù)中既包含所有病人的癌組織也包含部分病人的癌旁組織政基,因此文件數(shù)量多于實(shí)際病人數(shù)贞铣。此時(shí)獲得的臨床特征數(shù)據(jù)以病人為單位(n=349),小于下載總樣本數(shù)沮明,等于總樣本中癌癥樣本數(shù)辕坝。
注:overall survival的計(jì)算方法:臨床數(shù)據(jù)中包含days_to_death以及days_to_last_follow_up兩行數(shù)據(jù),僅vital status為Dead的樣本會(huì)有days_to_death數(shù)據(jù)荐健,而Alive的樣本記錄有days_to_last_follow_up酱畅,以及有極少數(shù)樣本兩項(xiàng)數(shù)據(jù)均缺失。因此在overall_survival中江场,死亡樣本以days_to_death為結(jié)局發(fā)生天數(shù)纺酸,而存活樣本以days_to_last_follow_up作為刪失發(fā)生天數(shù)。overall_survival與censoring_status配合址否,在生存分析中會(huì)用到餐蔬。
clinical_data_original <- jsonlite::fromJSON('clinical.cart.2020-03-08.json')
filter_follow_up <- function(x){
x$days_to_last_follow_up
}
filter_tumor_stage <- function(x){
x$tumor_stage
}
days_to_last_follow_up <- as.numeric(lapply(clinical_data_original$diagnoses, filter_follow_up))
tumor_stage <- as.character(lapply(clinical_data_original$diagnoses, filter_tumor_stage))
clinical_trait <- clinical_data_original$demographic[,c('submitter_id', 'gender', 'days_to_death', 'vital_status')] %>%
cbind(days_to_last_follow_up, tumor_stage)
clinical_trait$tumor_stage <- as.character(clinical_trait$tumor_stage)
clinical_trait <- tidyr::separate(clinical_trait, col= 'submitter_id',
into = c('submitter_id', 'junk'), sep='_', remove = T)[,-2]
clinical_trait <- clinical_trait[!duplicated(clinical_trait$submitter_id),]
clinical_trait[is.na(clinical_trait$days_to_death), 'days_to_death'] <- clinical_trait[is.na(clinical_trait$days_to_death), 'days_to_last_follow_up']
clinical_trait <- clinical_trait[,-5]
names(clinical_trait) <- c('submitter_id', 'gender', 'overall_survival', 'censoring_status', 'tumor_stage')
clinical_trait <- clinical_trait %>% dplyr::filter(!(is.na(overall_survival))) %>% dplyr::filter(tumor_stage != 'not reported')
- 通過(guò)TCGA_ID將所有樣品分組。
思路是:TCGA ID由dash隔開(kāi)佑附,第四個(gè)模塊記錄了樣本的分組樊诺。若該模塊中編號(hào)為01-09則為腫瘤,若編號(hào)在10-29則為正常癌旁樣本音同。
#classify the samples by TCGA_IDs (TCGA-XX-XXXX-NN: NN < 10 are cancer samples, otherwise are normal samples).
condition_table <- tidyr::separate(data=data.frame('ids'=colnames(expr_FPKM)[2:length(expr_FPKM)]),
col='ids',
sep='-',
into=c('c1','c2','c3','c4','c5','c6','c7'))
condition_table <- cbind('ids'=colnames(expr_FPKM)[2:length(expr_FPKM)], condition_table) %>%
cbind('submitter_IDs'=substr(colnames(expr_FPKM)[2:length(expr_FPKM)], 1, 12))
condition_table$c4 <- gsub(condition_table$c4, pattern = '^0..', replacement = 'cancer')
condition_table$c4 <- gsub(condition_table$c4, pattern = '^1..', replacement = 'normal')
condition_table <- condition_table[,c('ids','c4','submitter_IDs')]
names(condition_table) <- c('TCGA_IDs', 'sample_type', 'submitter_id')
condition_table$submitter_id <- as.character(condition_table$submitter_id)
condition_table$TCGA_IDs <- as.character(condition_table$TCGA_IDs)
- 去除表達(dá)矩陣中腫瘤或正常組內(nèi)可能存在的重復(fù)樣本词爬,并將cancer樣本與臨床數(shù)據(jù)匹配。
由于一個(gè)樣本可能測(cè)序多次权均,所以可以通過(guò)submitter_ID (也就是TCGA_ID的前12位)將組內(nèi)重復(fù)去除顿膨。我這里采取的辦法是保留第一次在表達(dá)矩陣中出現(xiàn)的樣本,去除隨后出現(xiàn)的重復(fù)樣本叽赊。此外由于此前有部分病人存在信息缺失虽惭,已被剔除,現(xiàn)也只保留含有全部信息的病人cancer樣本蛇尚。(normal樣本不需匹配臨床信息芽唇,因?yàn)閮H用于差異表達(dá)分析,不會(huì)涉及后續(xù)生存分析或WCGNA分析)
注:由于normal組的癌旁組織也來(lái)源于病人取劫,所以需要組內(nèi)分開(kāi)去除重復(fù)匆笤,否則來(lái)源于統(tǒng)一病人的normal和cancer樣本會(huì)被去除其一。
#remove the duplicated samples within groups and only maintain the cancer samples with clinical data.
condition_table_cancer <- condition_table[condition_table$sample_type=='cancer',]
condition_table_cancer <- condition_table_cancer[!duplicated(condition_table_cancer$submitter_id),] %>%
dplyr::inner_join(clinical_trait, by='submitter_id')
condition_table_normal <- condition_table[condition_table$sample_type=='normal',]
condition_table_normal <- condition_table_normal[!duplicated(condition_table_normal$submitter_id),]
condition_table <- rbind(condition_table_cancer[1:3], condition_table_normal)
- 獲得最終可以進(jìn)行下游分析的規(guī)范FPKM表達(dá)矩陣谱邪。
該矩陣第一列為GENE SYMBOL炮捧,第一行為樣本TCGA_ID,其余為FPKM表達(dá)數(shù)據(jù)惦银。
#get the final total FPKM expr matrix and expr matrix containing cancer samples for downstream cibersort analyses.
mRNA_exprSet_cancer_only <- mRNA_exprSet[,c('gene_name', condition_table_cancer$TCGA_IDs)]
mRNA_exprSet <- mRNA_exprSet[,c('gene_name', condition_table$TCGA_IDs)]
write.csv(mRNA_exprSet, '/Users/zx/Desktop/TCGA/TCGA-LUSC-white/FPKM/tutorial/FPKM_matrix_LUSC_.csv', quote = F, row.names = F)
到此為止TCGA表達(dá)矩陣的數(shù)據(jù)下載和清洗整合就全部結(jié)束了E乜巍D┦摹!
開(kāi)不開(kāi)心J轵健@琛!(不)
最后輸出的表達(dá)矩陣就是長(zhǎng)這樣的:
說(shuō)在最后:關(guān)于counts數(shù)據(jù)的下載...
幾乎一模一樣殊校,鏡像操作即可晴玖!唯一的區(qū)別是在初始表達(dá)矩陣整合完成后,需要去除dataframe的最后5行为流,該5行是樣本的總結(jié)性數(shù)據(jù)呕屎,在后續(xù)分析中不會(huì)用到,直接去除即可敬察。