TCGA-数据下载整理（二）
加入我们用TCGAbiolinks下载之后我们可以看一下另外一种操作
image.png
data_download_from_gdc
打开后发现是1208个文件夹，为什么是1208个，这跟我当时学的数据有关系，可以不用管让 每个文件夹里面还有个压缩文件，我们现在的任务就是，如何把每个文件夹里面的压缩文件放置在统一的文件夹中 假设我们的文件夹是这个,
首先设置工作目录到这一文件夹下：
00_data_read_in_one_file
并且他就是创建在 data_download_from_gdc中，所以现在文件夹的总体数目是应该是1208加1就是1209个 我们在R语言中输入
length(list.files())
这个命令中list.files()的意思就是列举当前工作目录中文件，length的意思就是有有多少个，一看，果真是1209个 其实
length(dir())
这个命令也能做这个事情。 为什么我的文件夹要起一个这么诡异的名字？ 以00开头？？
这是因为我需要让他处于文件中第一个位置，我们来验证一下：
dir()[1]
发现第一个文件夹就是他，这样我就每次循环访问后面的1208个文件夹，每次都把看到的东西复制过来就可以了
你说我为什么这么事无巨细，因为我可以很简洁，也可以很啰嗦，取决于我是否想要别人听的懂。
在R语言里面是可以直接创建文件夹的，鼠标右击也可以
dir.create("00_data_read_in_one_file") #创建新的文件夹,确保文件夹排在第一位
遍历和复制,为什么从2开始，因为第一个文件夹你已经知道了
for (dirname in dir()[2:length(dir())]){
file <- list.files(dirname,pattern = "*.counts") #找到对应文件夹中的内容，pattern可以是正则表达式
file.copy(paste0(dirname,"/",file),"00_data_read_in_one_file") #复制内容到新的文件夹
}
运行完了之后，我们可以打开文件夹看一下，确实在里面，这时候可以全部选择，
将文件解压到data_unzip文件夹,解压数据1.42GB，
我们发现即使这个样子，文件的名称也是怪怪的
那我们就来转换，转换的信息藏在metadata文件中，这个要去看开始的那个帖子下载
切换一下工作目录
setwd("~/skill_practice/BRCA_999")
注意一开始我就建立了 BRCA_999这个文件夹， data_download_from_gdc是他的子文件夹， 00_data_read_in_one_file又是 data_download_from_gdc的子文件夹 metadata是json格式的
读入json格式的文件,他是一个1208行，15列的数据框
metadata <- jsonlite::fromJSON("metadata.cart.2017-11-15T09_56_59.722935.json")
我们转换的信息就是两列filename和associatedentities，我们把它选出来
require(dplyr)
metadata_id <- metadata %>%
dplyr::select(c(file_name,associated_entities))
我需要的是 file_name和样本名称，样本名称藏在了 associated_entities 列表中里面包括了 entity_id, case_id, entity_submitter_id, entity_type这四个项目，查看第一个了解一下
metadata$associatedentities[1] [[1]]
entityid caseid 1 52033f64-1e6f-4657-a4fb-7cfeffc61951 39de7761-e762-4811-b95c-8216b79ae06b entitysubmitterid entitytype 1 TCGA-AN-A0XW-01A-11R-A109-07 aliquot
现在的想法是我把filename和 associated_entities中的 entity_submitter_id提取出来，做成一个数据框，然后我批量对应转换
naid_df <- data.frame()
for (i in 1:1208){
naid_df[i,1] <- substr(metadata_id$file_name[i],1,nchar(metadata_id$file_name[i])-3)
naid_df[i,2] <- metadata_id$associated_entities[i][[1]]$entity_submitter_id
}
现在把1208个小文件读入一个矩阵文件，并且给每一个文件加上 filename和 entity_submitter_id 论坛有道题目就是处理的这个问题
生信编程直播第四题：多个同样的行列式文件合并起来
http://www.biotrainee.com:8080/thread-603-1-1.html
我自己也给出了R语言的解法
在R语言中将多个同样的行列式文件合并起来
http://guoshipeng.com/2017/11/10/14-R-for-binding-colums/
但是当时不知道，TCGA的单个文件是没有列名的，导致无法合并，所以本次要复杂一点
读入所有解压的文件 1208个
nameList <- list.files("data_unzip/")
location <- which(naid_df==nameList[1],arr.ind = TRUE) ##which函数有一个已知value返回坐标的功能
TCGA_id <- as.character(naid_df[location[1],2]) ##通过坐标，获取TCGA_id
expr_df<- read.table(paste0("data_unzip/",nameList[1]),stringsAsFactors = F, header = F) #读入第一个文件，保存为data.frame
names(expr_df) <- c("gene_id",TCGA_id) #给刚才数据库命名
这边开始批量作业
for (i in 2:length(nameList)){
location <- which(naid_df==nameList[i],arr.ind = TRUE)
TCGA_id <- as.character(naid_df[location[1],2])
dfnew <- read.table(paste0("data_unzip/",nameList[i]),stringsAsFactors = F,header = F)
names(dfnew) <- c("gene_id",TCGA_id)
expr_df <- inner_join(expr_df,dfnew,by="gene_id")
}
晚上走的时候没运行完，早上来的时候已经完毕，限速环节应该是read.table，早上再来尝试运行一次总是说内存不够 我尝试了一下fread来解决这个问题：
require(data.table)
nameList <- list.files("data_unzip/")
location <- which(naid_df==nameList[1],arr.ind = TRUE)
TCGA_id <- as.character(naid_df[location[1],2])
expr_df<- fread(paste0("data_unzip/",nameList[1]))
names(expr_df) <- c("gene_id",TCGA_id)
for (i in 2:length(nameList)){
location <- which(naid_df==nameList[i],arr.ind = TRUE)
TCGA_id <- as.character(naid_df[location[1],2])
dfnew <- fread(paste0("data_unzip/",nameList[i]))
names(dfnew) <- c("gene_id",TCGA_id)
expr_df <- inner_join(expr_df,dfnew,by="gene_id")
}
结果大概2分钟搞定，速度喜人！！！！ 总共60488行,查看最后几行发现有5行不是我们要的
tail(expr_df$gene_id,10)
去掉最后五行
expr_df <- expr_df[1:(length(expr_df$gene_id)-5),]
保存数据，大概是3个G左右
save(expr_df,file = "expr_df.Rda")
下面开始id转换，信息在GTF文件中表达矩阵里面的gene_id有小数点， 而GTF文件中没有,调整一下，先以“.”分列，在去掉小数点后的列
require(dplyr)
require(tidyr)
expr_df_nopoint <- expr_df %>%
tidyr::separate(gene_id,into = c("gene_id","drop"),sep="\\.") %>%
dplyr::select(-drop)
save(expr_df_nopoint,file = "expr_df_nopoint.Rda")
load(file = "expr_df_nopoint.Rda")
下载GTF文件来注释 ftp://ftp.ensembl.org/pub/release-90/gtf/homo_sapiens
安装包：
source("https://bioconductor.org/biocLite.R")
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
biocLite("rtracklayer")
biocLite("SummarizedExperiment")
读入GTF数据
gtf1 <- rtracklayer::import('Homo_sapiens.GRCh38.90.chr.gtf')
gtf_df <- as.data.frame(gtf1)
保存数据
save(gtf_df,file = "gtf_df.Rda")
读入27个变量，2612129个观测，测试一下显示的不错
test <- gtf_df[1:5,]
View(test)
进展的很快，我们现在可以提取mRNA的表达矩阵啦， 以gtf文件中的 gene_biotype为标准，里面写 protein_coding的就是编码基因
首先要把这些基因提取出来，然后与表达谱融合，我在这个例子还提取了 gene_name, gene_id，所以最后的时候，我把三种表达方式合在了一起 这样，以后我无论用什么方式都可以选出我要的基因了
require(dplyr)
require(tidyr)
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 = " | ")
得到19688行，跟我们的认知很吻合
保存数据
save(mRNA_exprSet,file = "mRNA_exprSet.Rda")
下面进行标准化
suppressPackageStartupMessages(library(DESeq2))
suppressPackageStartupMessages(library(edgeR))
本次使用的方法是DESeq2中vst法 至于有哪几种方法，还有该选什么方法，明天会有一个帖子介绍
制作metadata，区别肿瘤和正常组，这时候需要对TCGA的id有个了解，从左往右数数，第14,15上的数字很重要 其中01-09是tumor，也就是癌症样本；其中10-29是normal，也就是癌旁
metadata <- data.frame(names(mRNA_exprSet)[-1])
for (i in 1:length(metadata[,1])) {
num <- as.numeric(substring(metadata[i,1],14,15))
if (num %in% seq(1,9)) {metadata[i,2] <- "T"}
if (num %in% seq(10,29)) {metadata[i,2] <- "N"}
}
加个名称
names(metadata) <- c("TCGA_id","sample")
将sample转化成factor
metadata$sample <- as.factor(metadata$sample)
我们可总结一下,有112例normal，1096例肿瘤组织
metadata %>% dplyr::group_by(sample) %>% summarise(n())
保存数据，我想问一下，为什么你老是要保存数据，因为我眼里常含泪花
save(metadata,file = "metadata.Rda")
制作countData
mycounts <- mRNA_exprSet
keepGene=rowSums(edgeR::cpm(mycounts[-1])>0) >=2
table(keepGene);dim(mycounts)
这是数据维度19668 和1209
dim(mycounts[keepGene,])
剩余19546 1209
mycounts <-mycounts[keepGene,]
构建 DESeq对象
dds <-DataSetFromMatrix(countData=mycounts,
colData=metadata,
design=~sample,
tidy=TRUE)
标准化vst法，5分钟
vsd <- vst(dds, blind = FALSE)
mRNA_exprSet_vst <- as.data.frame(assay(vsd))
save(mRNA_exprSet_vst,file = "mRNA_exprSet_vst.Rda")
load(file = "dds.Rda")
下面就是完整的过程
dds <- DESeq(dds)
save(dds,file = "dds_DEseq.Rda")
res <- results(dds, tidy=TRUE) #获取结果
res <- as_tibble(res)
require(dplyr)
res <- res %>%
separate(row,c("symbol","ensemble","genetype"),sep = " \\| ") %>%
dplyr::select(- c(ensemble,genetype)) %>%
arrange(desc(abs(log2FoldChange))) %>% #排序。为了去重
distinct(symbol,.keep_all = TRUE) %>%
arrange(desc(log2FoldChange))#再次按照log2FoldChange从大到小排序
save(res,file = "result.Rda")
既然做了差异分析，那对下游基因注释就是顺手 的事情啦，我们没有选择，就用Y叔的包，clusterprofiler 在这之前我们需要获得差异基因列表 这一步，就是基本数据框的基本操作，按照pvalue和fc来排序，选择自己一定数量的基因，数量你来定，最终得到基因列表gene(我没有演示，需要自己做)
加载包，加载之前确定自己装了
library(DOSE)
library(GO.db)
library(org.Hs.eg.db)
library(topGO)
library(GSEABase)
library(clusterProfiler)
基因名称转换，返回的是数据框
gene = bitr(gene, fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db")
head(gene)
GO分析 细胞组分
ego_CC <- enrichGO(gene = gene$ENTREZID,
OrgDb= org.Hs.eg.db,
ont = "CC",
pAdjustMethod = "BH",
minGSSize = 1,
pvalueCutoff = 0.01,
qvalueCutoff = 0.01,
readable = TRUE)
生物过程
ego_BP <- enrichGO(gene = gene$ENTREZID,
OrgDb= org.Hs.eg.db,
ont = "BP",
pAdjustMethod = "BH",
minGSSize = 1,
pvalueCutoff = 0.01,
qvalueCutoff = 0.01,
readable = TRUE)
分子功能：
ego_MF <- enrichGO(gene = gene$ENTREZID,
OrgDb= org.Hs.eg.db,
ont = "MF",
pAdjustMethod = "BH",
minGSSize = 1,
pvalueCutoff = 0.01,
qvalueCutoff = 0.01,
readable = TRUE)
作图
barplot(ego_CC, showCategory=20,title="EnrichmentGO_CC")#条状图，按p从小到大排的
dotplot(ego_BP,title="EnrichmentGO_BP_dot")#点图，按富集的数从大到小的
KEGG分析
kk <- enrichKEGG(gene = gene$ENTREZID,
organism ="human",
pvalueCutoff = 0.01,
qvalueCutoff = 0.01,
minGSSize = 1,
#readable = TRUE ,
use_internal_data =FALSE)
作图
barplot(kk)
dotplot(kk)
