• 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()
}

参考mayoneday