0.安装包
rm(list=ls())
if (!requireNamespace("BiocManager", quietly=TRUE))
install.packages("BiocManager")
BiocManager::install("dplyr")
BiocManager::install("rtracklayer")
BiocManager::install("SummarizedExperiment")
BiocManager::install("DESeq2")
BiocManager::install("tidyr")
library(dplyr)
library(rtracklayer)
library(SummarizedExperiment)
library(DESeq2)
library(tidyr)
library(R.utils)###解压缩函数gunzip()
library(jsonlite)
library(clusterProfiler)
library(org.Hs.eg.db)
library(AnnotationDbi)
library(edgeR)
1、下载好HTseq-FPKM的TCGA数据
#查看目录下files文件数目是否和网页上一直
length(list.files())
2、将所有文件夹下的压缩包复制到一个文件夹下
######创建新的文件夹,确保文件夹排在第一位
dir.create("00_data_read_in_one_file")
#####确认第一个文件夹是刚刚所创
dir()[1]
####批量循环操作将所有.gz格式的压缩文件导入到新文件夹中
for (dirname in dir()[2:length(dir())]){
file <- list.files(dirname,pattern = "*.gz") #找到对应文件夹中的内容，pattern可以是正则表达式
file.copy(paste0(dirname,"/",file),"00_data_read_in_one_file") #复制内容到新的文件夹
}
3、解压缩到同一文件夹#gunzip()函数
####修改工作目录到新创文件夹"00_data_read_in_one_file"
####批量解压缩00_data_read_in_one_file文件夹下的压缩包,然后复制到gunziped文件夹下
dir.create("gunziped")
####循环批量操作
for (i in 1:length(dir())) {
foo <- list.files(pattern = "*.gz")[i]
faa <- gunzip(foo)
file.copy(faa,"gunziped",overwrite = TRUE)
}
4、提取藏在metadata文件中的转换信息
metadata <- jsonlite::fromJSON("metadata.cart.2020-03-20.json")###网页上下载
dim(metadata)#读入json格式的文件,该文件是一个407行，14列的数据框
#[1] 407 14
colnames(metadata)#有如下列
#[1] "experimental_strategy" "associated_entities" "access"
#[4] "data_type" "file_name" "file_id"
#[7] "state" "data_format" "data_category"
#[10] "file_size" "md5sum" "analysis"
#[13] "submitter_id" "annotations"
#方法：
###### 提取filename和 associated_entities中的 entity_submitter_id出来
######做成一个数据框
######然后批量对应转换
##转换的信息就是两列filename和associated_entities，我们把它选出来
metadata$associated_entities[1][[1]]#有四项
##entity_submitter_id entity_type
##TCGA-BR-8366-01A-11R-2343-13 aliquot
##entity_id case_id
##f461ac7f-a73e-4291-9391-10ef93ef4240 66c04b8f-42fd-44c5-945b-ee149bcf5dad
#4.2.1 提取filename和样本名称(样本名就是associated_entities中的entity_submitter_id)
metadata_id <- metadata %>% dplyr::select(c(file_name,associated_entities))
View(metadata_id )
#在associated_entities列表中里面包括了 entity_id, case_id, entity_submitter_id, entity_type这四个项目
#4.2.2 做成一个数据框
##把filename和 associated_entities中的 entity_submitter_id提取出来，做成一个数据框，然后我批量对应转换
naid_df <- data.frame()
###将工作目录转到 "gunziped" 中
##4.2.3 批量对应转换
for (i in 1:407){
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
}
#View(naid_df)
#head(naid_df)
save(naid_df,file = "naid_df.Rda")
load(file = "naid_df.Rda")
5、提取解压后407个txt文件中的信息，并制成表达矩阵
-
5.1读取该目录下的407个文件,先用一个运行一遍，然后建数据框，接着通过gene_id批量关联起来
#转工作目录
###提取一个，将数据框
nameList <- list.files("./")#确认目录下文件数
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("./",nameList[1]),stringsAsFactors = F, header = F)
#读入第一个文件，保存为data.frame
names(expr_df) <- c("gene_id",TCGA_id) #给刚才数据库命名
#View(expr_df)
for (i in 2:length(nameList)){
location <- which(naid_df==nameList[i],arr.ind = TRUE)
###根据文件名，返回该文件在naid_df中的坐标
TCGA_id <- as.character(naid_df[location[1],2])
###根据返回坐标的行数为行，列为2，再字符串处理，提取TCGA-id号
dfnew <- read.table(paste0("./",nameList[i]),stringsAsFactors = F,header = F)
###读取每一个txt文件下ENSG号和表达量
names(dfnew) <- c("gene_id",TCGA_id)
###将dfnew的两行分别命名为基因名和TCGA_id号
expr_df <- inner_join(expr_df,dfnew,by="gene_id")
#通过上面的与建库的expr_df关联起来
}
View(expr_df)
head(expr_df)#tail()
save(expr_df,file = "expr_df.Rda")
load(file = "expr_df.Rda")
6、 id转换
-
6.1 去掉gene_id的ensg_id中小数点及后面的版本数；
##因为expr_df表达矩阵里面的gene_id有小数点，而注释文件中没有,调整一下，先以“.”分列，再去掉小数点后的列
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")
#View(expr_df_nopoint)
load(file = "expr_df_nopoint.Rda")
keytypes(org.Hs.eg.db)#查看ID转换的类型
char <- expr_df_nopoint$gene_id
gen_ids <- bitr(char,fromType = "ENSEMBL",toType = c("SYMBOL","GENENAME"),OrgDb = "org.Hs.eg.db")#Y叔的clusterProfiler包
colnames(gen_ids)[1] <-c("gene_id")
gen_ids <- gen_ids[,1:2]
expr_id <- merge(gen_ids,expr_df_nopoint,by="gene_id")
#View(expr_id)#注释结果
save(expr_id,file = "expr_id.Rda")
7.差异分析
-
7.1 去掉有ensg号的列，保存基因名symbol列
load(file = "expr_id.Rda")#加载保存的数据
expr_gena <- expr_id[,(-1)]
-
7.2 去掉重复基因名行（或者取中位数、平均数等处理）
expr_gena <- expr_gena[!duplicated(expr_id$SYMBOL),]#或者函数uniqe()
rownames(expr_gena) <- expr_gena[,1] #将symbol列设置为行名
exprSet <- expr_gena[,(-1)]
#View(exprSet)
7.4 将FPKM改成TPM格式，方便后续差异分析
###FPKM确实存在不准确性，推荐使用TPM
###链接地址：//www.greatytc.com/p/9dfb65e405e8
fpkmToTpm <- function(fpkm)
{
exp(log(fpkm) - log(sum(fpkm)) + log(1e6))
}
exprSet_tpms <- apply(exprSet,2,fpkmToTpm)#如果要计算TPM值，只需要用一下apply函数
save(exprSet_tpms,file ="exprSet_tpms.Rda")
load(file ="exprSet_tpms.Rda")
colSums(exprSet_tpms)#根据TPM的特征进行检查，看每列加和是否一致
View(exprSet_tpms)
