今天繼續(xù)我們的TCGA數(shù)據(jù)分析系列。
TCGA數(shù)據(jù)下載
TCGA數(shù)據(jù)下載的方式有很多，本次我們利用UCSC Xena數(shù)據(jù)庫下載數(shù)據(jù)，網(wǎng)址是：https://xenabrowser.net/。該平臺(tái)內(nèi)置了一些公共數(shù)據(jù)集，比如來自TCGA, ICGC等大型癌癥研究項(xiàng)目的數(shù)據(jù)，不僅可以對(duì)數(shù)據(jù)進(jìn)行分析，而且還提供了對(duì)應(yīng)文件的下載功能。
打開后界面是這樣的
點(diǎn)擊DATA SETS，里面有很多數(shù)據(jù)集，我們選擇TCGA肝癌數(shù)據(jù)
接著我們選擇HTSeq-Count
這里可以看到值log2（count+1），為什么加一呢，因?yàn)楹芏嗷虻谋磉_(dá)值是0，無法取log。
下載下來，解壓后打開是這個(gè)樣子
接下來我們進(jìn)行差異分析
首先加載R包
rm(list= ls())#一鍵清空#安裝并加載R包if(length(getOption("CRAN"))==0) options(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/")if(!require("rtracklayer")) BiocManager::install("rtracklayer")if(!require("tidyr")) BiocManager::install("tidyr")if(!require("dplyr")) BiocManager::install("dplyr")if(!require("DESeq2")) BiocManager::install("DESeq2")if(!require("limma")) BiocManager::install("limma")if(!require("edgeR")) BiocManager::install("edgeR")
讀取數(shù)據(jù)
#導(dǎo)入表達(dá)譜數(shù)據(jù)LIHCdata=read.table("TCGA-LIHC.htseq_counts.tsv",header=T,sep='\t')LIHCdata[1:4,1:4]
利用我們之前講到的方法去掉ensemble ID的點(diǎn)號(hào)R語言學(xué)習(xí)系列之separate {tidyr}
LIHCdata1<-separate(LIHCdata,Ensembl_ID,into= c("Ensembl_ID"),sep="\\.")LIHCdata1[1:4,1:4]
利用我們之前講到的方法去掉ensemble ID的點(diǎn)號(hào)R語言學(xué)習(xí)系列之separate {tidyr}
LIHCdata1<-separate(LIHCdata,Ensembl_ID,into= c("Ensembl_ID"),sep="\\.")LIHCdata1[1:4,1:4]
接下來我們需要對(duì)ID進(jìn)行轉(zhuǎn)換，轉(zhuǎn)換的方法也有很多，有R包，在線數(shù)據(jù)庫。小工具等，這里我們通過下載最新版的GTF文件來進(jìn)行轉(zhuǎn)換。
首先，打開ensembl網(wǎng)址：http://asia.ensembl.org/index.html
點(diǎn)擊Download
再點(diǎn)擊Download database
再點(diǎn)擊human的GTF文件
下載Homo_sapiens.GRCh38.99.chr.gtf.gz文件
然后放到我們R語言的工作目錄內(nèi)，打開文件
gtf <- rtracklayer::import('Homo_sapiens.GRCh38.99.chr.gtf.gz')#轉(zhuǎn)換為數(shù)據(jù)框gtf <-as.data.frame(gtf)
查看文件，保存文件為Rdata，將來方便我們直接打開
dim(gtf)save(gtf,file ="Homo_sapiens.GRCh38.99基因組注釋文件.Rda")
可見文件有290萬行，27列，由于GTF文件中，基因ID的列名是gene_id,所以我們把LIHCdata1中的基因列名改成一樣的，方便后續(xù)合并
colnames(LIHCdata1)[1]<-'gene_id'
通過瀏覽文件看到我們需要的主要信息在
1 type，這一列我們需要選擇gene
2 gene_biotype，這一列我們需要選擇protein_coding，當(dāng)然你也可以選擇其他的種類，比如miRNA，長鏈非編碼等。
所以我們首先把蛋白編碼的基因的行都篩選出來
a=dplyr::filter(gtf,type=="gene",gene_biotype=="protein_coding")dim(a)
這個(gè)時(shí)候a文件只有19939行了，列下來我們只選擇gene_name，gene_id和gene_biotype這三列，其他都不要了
b=dplyr::select(a,c(gene_name,gene_id,gene_biotype))b[1:4,]
接下來用LIHCdata1和b文件中共有的gene_id列還合并文件
c=dplyr::inner_join(b,LIHCdata1,by ="gene_id")c[1:5,1:5]
再去掉第2，3列，基因名再去重
d=select(c,-gene_id,-gene_biotype)mRNAdata=distinct(d,gene_name,.keep_all = T)
把行名由數(shù)字換成基因
rownames(mRNAdata)<-mRNAdata[,1]mRNAdata<-mRNAdata[,-1]
我們下載的數(shù)據(jù)取過了log2（count+1），這里我們?cè)俜祷豤ount
mRNAdata<-2^mRNAdata -1
保存文件，大功告成！
write.csv(mRNAdata,"mRNAdata.csv",quote = F,row.names = T)save(mRNAdata,file ="mRNAdata.Rda")
好了，今天先講到這，下回我們來進(jìn)行后續(xù)的差異分析。
需要今天下載的兩個(gè)文件的話，公眾號(hào)回復(fù)“大功告成”可獲取鏈接。
另外，最近收集了一些很好的資源，想分享給大家，順便能漲一些粉，主要有
1. 19年中標(biāo)的各門類國自然題目匯總，以及17年的國自然匯總，部分含摘要！
2. R語言學(xué)習(xí)書籍
R語言實(shí)戰(zhàn)(中文完整版)
R數(shù)據(jù)科學(xué)(中文完整版）
ggplot2：****數(shù)據(jù)分析與圖形藝術(shù)
30分鐘學(xué)會(huì)ggplot2
3. TCGA數(shù)據(jù)整理
前期從https://xenabrowser.net/datapages/ (UCSC Xena)數(shù)據(jù)庫下載的TCGA數(shù)據(jù)，傳到了百度云上備份。
感興趣的話，轉(zhuǎn)發(fā)朋友圈或者100人以上的微信群，截圖發(fā)到公眾號(hào)，即可獲取全部資源的百度云鏈接，鏈接7天有效，希望大家趕緊下載。你們的支持是我前進(jìn)的動(dòng)力，感謝。
