TCGA-数据整理（counts =》表达矩阵）
2020-03-28 本文已影响0人
养猪场小老板
-
0、R包安装
rm(list=ls())
if (!requireNamespace("BiocManager", quietly=TRUE))
install.packages("BiocManager")
BiocManager::install("rtracklayer")
BiocManager::install("clusterProfiler")
BiocManager::install("org.Hs.eg.db")
BiocManager::install("rtracklayer")
BiocManager::install("DESeq2")
BiocManager::install("DOSE")
BiocManager::install("enrichplot")
BiocManager::install("edgeR")
BiocManager::install("BiocGenerics")
BiocManager::install("ggplot2")
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)
library(BiocGenerics)
library(ggplot2)
#1、下载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)
}
-
1、差异表达矩阵
#############################################################
###下载好的数据压缩包全部提取到gunziped的文件夹里#############
#4.1##############
metadata <- jsonlite::fromJSON("metadata.cart.2020-03-25.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"
#4.2
#方法：
###### 提取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 做成一个数据框
将工作目录转到 "gunziped" 中
##把filename和 associated_entities中的 entity_submitter_id提取出来，做成一个数据框，然后我批量对应转换
naid_df <- data.frame()
###
##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)#(407行才对)
###换工作目录
###提取一个，将数据框
colnames(naid_df) <- c("filename","TCGA_id")#添加列名
test <- data.table::fread(paste0("./",naid_df$filename[1]))#先读取一个，查看效果
View(test)
expr_df <- data.frame(matrix(NA,nrow(test),nrow(naid_df)))#先建一个数据框，行ensemble_id号，列为样本名TCGA_id号
View(expr_df)
for (i in 1:nrow(naid_df)) {
print(i)
expr_df[,i]<- data.table::fread(paste0("./",naid_df$filename[i]))[,2]#把每次读取出来的第二列数据变成一列
}
View(expr_df)
#用对应的TCGA id 给获得的数据命名
colnames(expr_df) <- naid_df$TCGA_id
View(expr_df)
#读取第一个文件，将其第一列合并到大数据框上，作为基因id
gene_id <- data.table::fread(paste0("./",naid_df$filename[1]))$V1
gene_id
expr_df <- cbind(gene_id=gene_id,expr_df)#cbind() 把矩阵横向合并成一个大矩阵（列方式），而rbind()是纵向合并（行方式）
View(expr_df)
#去除最后5行，看下面就知道了
tail(expr_df$gene_id,10)
expr_df <- expr_df[1:(nrow(expr_df)-5),]
View(expr_df)
save(expr_df,file = 'expr_df.Rdata')
load(file = 'expr_df.Rdata')
#去掉版本号
library(tidyr)
expr_df_nopoint <- expr_df %>%
tidyr::separate(gene_id,into = c("gene_id"),sep="\\.")
View(expr_df_nopoint)
{
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")
load(file = "expr_id.Rda")#加载保存的数据
write.csv(expr_id,file = "expr_id.csv",row.names = F)
expr_gena <- expr_id[,(-1)]#去掉ensemble_id列
expr_gena <- expr_gena[!duplicated(expr_id$SYMBOL),]#或者函数uniqe()
rownames(expr_gena) <- expr_gena[,1] #将symbol列设置为行名
exprSet <- expr_gena[,(-1)]#去掉symbol列
exprSet <- exprSet[!apply(exprSet,1,function(x){sum(floor(x)==0)>3}),]#有三个表达为0，即去除。
exprSet <- na.omit(exprSet)#去NA
#View(expr_gena)
#View(exprSet)
save(exprSet,file="exprSet,Rda")
load(file="exprSet,Rda")
################################
group_list=ifelse(as.numeric(substr(colnames(exprSet),14,15)) < 10,'tumor','normal')#分组
table(group_list)#查看
-
2、基因差异分析(三选一就可以进行后续分析)
测序的数据，DESeq2 ，edgeR一般用的比较多，接受的是RNAseq的count数据，非log的，如果数据已经log则要取log。limma是最经典的，但是limma是必须接受log之后的值，才能正确算出差异，一般芯片数据用limma包，（有些从GEO数据库下载的数据，经过标准化处理的时候是已经log过了的就不用log了）
1)DESeq2
if(T){
library(DESeq2)
group_list<-as.factor(group_list)#把group_list变成因子，因为它接受的是因子
colData <- data.frame(row.names=colnames(exprSet, group_list=group_list) )#得到样本和分组对应的一个表
dds <- DESeqDataSetFromMatrix(countData = exprSet,
colData = colData,
design = ~ group_list)#构建模型
dds <- DESeq(dds)#差异分析前处理
res <- results(dds,
contrast=c("group_list","tumor","normal"))
#进行对比后提取数据，用results提取是因为看了说明书此函数的差异分析结果就是用这个提取
resOrdered <- res[order(res$padj),]#按照差异性大小排序
head(resOrdered)
DEG =as.data.frame(resOrdered)#变为数据框
DESeq2_DEG = na.omit(DEG) #删除NA值
nrDEG=DESeq2_DEG[,c(2,6)]#2为log2FoldChange 6为 pvalue 列，取出这两个
colnames(nrDEG)=c('log2FoldChange','pvalue')
}
2)edgeR
if(T){
library(edgeR)
d <- DGEList(counts=exprSet,group=factor(group_list))#构建模型，只需要两个因素
keep <- rowSums(cpm(d)>1) >= 2#做筛选，此方法需要这样做
table(keep)
d <- d[keep, , keep.lib.sizes=FALSE]
d$samples$lib.size <- colSums(d$counts)
d <- calcNormFactors(d)
d$samples
dge=d
design <- model.matrix(~0+factor(group_list))
rownames(design)<-colnames(dge)
colnames(design)<-levels(factor(group_list))
dge=d
dge <- estimateGLMCommonDisp(dge,design)
dge <- estimateGLMTrendedDisp(dge, design)
dge <- estimateGLMTagwiseDisp(dge, design)
fit <- glmFit(dge, design)
# https://www.biostars.org/p/110861/
lrt <- glmLRT(fit, contrast=c(-1,1)) #因为是两组所以这么做，按照说明书
nrDEG=topTags(lrt, n=nrow(dge))
nrDEG=as.data.frame(nrDEG)
head(nrDEG)
edgeR_DEG =nrDEG
nrDEG=edgeR_DEG
}
need_DEG=nrDEG#View(need_DEG),colnames(need_DEG)
-
3）voom from limma
if(T){
suppressMessages(library(limma))
design <- model.matrix(~0+factor(group_list))
colnames(design)=levels(factor(group_list))
rownames(design)=colnames(exprSet)
design
dge <- DGEList(counts=exprSet)
dge <- calcNormFactors(dge)
logCPM <- cpm(dge, log=TRUE, prior.count=3)#注意limma包的表达数据是需要经过log的
v <- voom(dge,design,plot=TRUE, normalize="quantile")#说明书需要 normalize
fit <- lmFit(v, design)
group_list
cont.matrix=makeContrasts(contrasts=c('tumor-normal'),levels = design)#注意这里不要反了，反了就是相反的结果
fit2=contrasts.fit(fit,cont.matrix)
fit2=eBayes(fit2)
tempOutput = topTable(fit2, coef='tumor-normal', n=Inf)
DEG_limma_voom = na.omit(tempOutput)
head(DEG_limma_voom)
nrDEG=DEG_limma_voom[,c(1,4)]
colnames(nrDEG)=c('log2FoldChange','pvalue')
draw_h_v(exprSet,nrDEG,'limma',group_list,1)
}
#一路跑下来最终得到结果
tmp_f=file.path(Rdata_dir,'TCGA-STAD-DEG_results.Rdata')
if(file.exists(tmp_f)){
save(DEG_limma_voom,DESeq2_DEG,edgeR_DEG, file = tmp_f)
}else{
load(file = tmp_f)
}
nrDEG1=DEG_limma_voom[,c(1,4)]
colnames(nrDEG1)=c('log2FoldChange','pvalue')
nrDEG2=edgeR_DEG[,c(1,5)]
colnames(nrDEG2)=c('log2FoldChange','pvalue')
nrDEG3=DESeq2_DEG[,c(2,6)]
colnames(nrDEG3)=c('log2FoldChange','pvalue')
mi=unique(c(rownames(nrDEG1),rownames(nrDEG1),rownames(nrDEG1)))
lf=data.frame(lf1=nrDEG1[mi,1],
lf2=nrDEG2[mi,1],
lf3=nrDEG3[mi,1])
cor(na.omit(lf))
-
3、接下来，画热图和火山图
draw_h_v <- function(exprSet,need_DEG,n='DEseq2',group_list,logFC_cutoff){
## we only need two columns of DEG, which are log2FoldChange and pvalue
## heatmap
need_DEG=nrDEG
library(pheatmap)
choose_gene=head(rownames(need_DEG),50) ## 50 maybe better
choose_matrix=exprSet[choose_gene,]#提取前50的表达矩阵
choose_matrix[1:4,1:4]
choose_matrix=t(scale(t(log2(choose_matrix+1)))) #画热图要先归一化
## http://www.bio-info-trainee.com/1980.html
annotation_col = data.frame( group_list=group_list )#变成一个表
rownames(annotation_col)=colnames(exprSet)#再取个名字
pheatmap(choose_matrix,show_colnames = F,annotation_col = annotation_col,
filename = paste0(1,'_need_DEG_top50_heatmap.png'))
library(ggfortify)
df=as.data.frame(t(choose_matrix))#行列转换
choose_matrix
df$group=group_list
png(paste0(2,'_DEG_top50_pca.png'),res=120)#存为什么文件
p=autoplot(prcomp( df[,1:(ncol(df)-1)] ), data=df,colour = "group")+theme_bw()#画图
print(p)
dev.off()
#画火山图
if(! logFC_cutoff){
logFC_cutoff <- with(need_DEG,mean(abs( log2FoldChange)) + 2*sd(abs( log2FoldChange)) )
}#如果没有给定阈值（差异倍数认为多少是上下调），则用平均值加上两倍的方差
logFC_cutoff=1
need_DEG$change = as.factor(ifelse(need_DEG$pvalue < 0.05 & abs(need_DEG$log2FoldChange) > logFC_cutoff,
ifelse(need_DEG$log2FoldChange > logFC_cutoff ,'UP','DOWN'),'NOT')
)
this_tile <- paste0('Cutoff for logFC is ',round(logFC_cutoff,3),
'\nThe number of up gene is ',nrow(need_DEG[need_DEG$change =='UP',]) ,
'\nThe number of down gene is ',nrow(need_DEG[need_DEG$change =='DOWN',])
)
library(ggplot2)
g = ggplot(data=need_DEG,
aes(x=log2FoldChange, y=-log10(pvalue),
color=change)) +
geom_point(alpha=0.4, size=1.75) +
theme_set(theme_set(theme_bw(base_size=20)))+
xlab("log2 fold change") + ylab("-log10 p-value") +
ggtitle( this_tile ) + theme(plot.title = element_text(size=15,hjust = 0.5))+
scale_colour_manual(values = c('blue','black','red')) ## corresponding to the levels(res$change)
print(g)
ggsave(g,filename = paste0(3,'_volcano.png'))
dev.off()
}
