今天继续我们的TCGA数据分析系列。
TCGA数据下载的方式有很多,本次我们利用UCSC Xena数据库下载数据,网址是:https:///。该平台内置了一些公共数据集,比如来自TCGA, ICGC等大型癌症研究项目的数据,不仅可以对数据进行分析,而且还提供了对应文件的下载功能。 打开后界面是这样的 点击DATA SETS,里面有很多数据集,我们选择TCGA肝癌数据 接着我们选择HTSeq-Count 这里可以看到值log2(count+1),为什么加一呢,因为很多基因的表达值是0,无法取log。 下载下来,解压后打开是这个样子 接下来我们进行差异分析 首先加载R包 rm(list = ls())#一键清空 #安装并加载R包 if(length(getOption("CRAN"))==0) options(CRAN="https://mirrors.tuna./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") 读取数据 #导入表达谱数据 LIHCdata=read.table("TCGA-LIHC.htseq_counts.tsv",header=T,sep='\t') LIHCdata[1:4,1:4] 利用我们之前讲到的方法去掉ensemble ID的点号R语言学习系列之separate {tidyr} LIHCdata1<-separate(LIHCdata,Ensembl_ID,into = c("Ensembl_ID"),sep="\\.") LIHCdata1[1:4,1:4] 接下来我们需要对ID进行转换,转换的方法也有很多,有R包,在线数据库。小工具等,这里我们通过下载最新版的GTF文件来进行转换。 首先,打开ensembl网址:http://asia./index.html 点击Download 再点击Download database 再点击human的GTF文件 下载Homo_sapiens.GRCh38.99.chr.gtf.gz文件 然后放到我们R语言的工作目录内,打开文件 gtf <- rtracklayer::import('Homo_sapiens.GRCh38.99.chr.gtf.gz') #转换为数据框 gtf <- as.data.frame(gtf) 查看文件,保存文件为Rdata,将来方便我们直接打开 dim(gtf) save(gtf,file = "Homo_sapiens.GRCh38.99基因组注释文件.Rda") 可见文件有290万行,27列,由于GTF文件中,基因ID的列名是gene_id,所以我们把LIHCdata1中的基因列名改成一样的,方便后续合并 colnames(LIHCdata1)[1]<-'gene_id' 通过浏览文件看到我们需要的主要信息在 1 type,这一列我们需要选择gene 2 gene_biotype,这一列我们需要选择protein_coding,当然你也可以选择其他的种类,比如miRNA,长链非编码等。 所以我们首先把蛋白编码的基因的行都筛选出来 a=dplyr::filter(gtf,type=="gene",gene_biotype=="protein_coding") dim(a) 这个时候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) 把行名由数字换成基因 rownames(mRNAdata)<- mRNAdata[,1] mRNAdata<-mRNAdata[,-1] 我们下载的数据取过了log2(count+1),这里我们再返回count mRNAdata <- 2^mRNAdata -1 保存文件,大功告成! write.csv(mRNAdata,"mRNAdata.csv",quote = F,row.names = T) save(mRNAdata,file = "mRNAdata.Rda") 好了,今天先讲到这,下回我们来进行后续的差异分析。 |
|