如何用前面的步骤进行其他数据的应用，现在是应用于TCGA 的RNA count数据进行的可视化
“#”为改动的地方
1.新建文件夹，对下述文件进行等同迁移

image.png

2.从网页下载数据
示例miRNA-Seq，clinical数据和joson没有区别
gdc tcga网站——files——data type☑️miRNA Expression qualification——experimental strategy ☑️miRNA-Seq——workflow type☑️BCGSC miRNA Profiling
case——program☑️TCGA——project☑️所需癌症类型——
manifest进行下载，储存至新建文件夹
3.打开tcga.1.Rproj
复制需要文件，至新建文件夹

image.png

options(stringsAsFactors = F)
library(stringr)
cancer_type="TCGA-CHOL"  #癌症类型改动
if(!dir.exists("clinical"))dir.create("clinical")    ##新建空文件夹
if(!dir.exists("expdata"))dir.create("expdata")
dir()  
#下面两行命令在terminal完成，根据前面dir看所有文件进行命令改动，改为该数据clinical/expdata文件名
#./gdc-client download -m gdc_manifest_cl.2020-03-23.txt -d clinical
#./gdc-client download -m gdc_manifest_expdata.2020-03-23.txt -d expdata
#gdc_manifest_cl.2020-03-23.txt 
#gdc_manifest_expdata.2020-03-23.txt
###将网站下载的清单数据下载至新建文件夹
length(dir("./clinical/"))
length(dir("./expdata/"))

可以看到，下载的文件是按样本存放的，我们需要得到的是表格，需要将他们批量读入R语言并整理。

3.整理临床信息

library(XML)  
result <- xmlParse("./clinical/142aea0e-7a7b-4ac4-9dbb-0f62e2379599/nationwidechildrens.org_clinical.TCGA-W5-AA2O.xml")  
#对应更改clinical里面的文件名

rootnode <- xmlRoot(result)  
xmlSize(rootnode)

print(rootnode[1]) 
#print(rootnode[2])  #患者有用的临床信息大多都在第二个节点
xmldataframe <- xmlToDataFrame(rootnode[2])  
head(t(xmlToDataFrame(rootnode[2]))) 

xmls = dir("clinical/",pattern = "*.xml$",recursive = T)
cl = list()   
for(i in 1:length(xmls)){
  result <- xmlParse(paste0("clinical/",xmls[[i]]))
  rootnode <- xmlRoot(result)
  cl[[i]] = xmlToDataFrame(rootnode[2])
}
clinical <- do.call(rbind,cl)
clinical[1:3,1:3]

4.整理表达矩阵

探索数据：先任选两个counts文件读取，并观察geneid的顺序是否一致。

options(stringsAsFactors = F)
x = read.table("expdata/03aee74e-4e37-4a58-a720-c90e807d2f40/be5bc6a0-9720-47ac-953e-fa8d0c32cd82.htseq.counts.gz",
               row.names = 1,sep = "\t") 
 #更改expdata文件中数据文件名（自动补齐）
x2 = read.table("expdata/10d08172-48d2-49e7-b760-721163492cc1/c1071bcd-5a0c-4e09-a578-fc4b6dbe26ad.htseq.counts.gz", 
                row.names = 1,sep = "\t") 
 #更改expdata文件中数据文件名（自动补齐）
identical(rownames(x),rownames(x2)) 

看一眼读取出来的数据是否符合格式，并进行处理，只需要count这一列：[1]，首行是否为数据，：header=T

image.png

由此可知，他们的geneid顺序是一致的，可以直接cbind，不会导致顺序错乱。

批量读取所有的counts.gz文件（将不同患者的表达矩阵所处文件名汇总到一个列表）

count_files = dir("expdata/",pattern = "*.htseq.counts.gz$",recursive = T)    
#观察expdata文件夹中数据格式pattern = "*.htseq.counts.gz$"并对应改变，找出数据共同表达后缀，进行替换

前面是文件夹的名字/后面是文件的名字

image.png

exp = list()
for(i in 1:length(count_files)){
  exp[[i]] <- read.table(paste0("expdata/",count_files[[i]]),row.names = 1,sep = "\t")  
#如果前面x,x2那里数据有改动，这里也要改动
}
exp <- do.call(cbind,exp)
dim(exp) #检查exp
exp[1:4,1:4]

image.png

image.png

发现问题：这样产生出来的表达矩阵没有列名（也就是不同患者id信息无法显示，之后就与临床信息无法对应）

解决办法：找到一个文件名与样本ID一一对应的文件。cart-json文件。
返回gdc网页，跟前面一样的步骤，case选种类，files选基因表达矩阵，勾选内容一致，✅加入购物车add all files to chart ✅右上角购物车图标☑️metadata自动下载数据，格式为cart-json。

meta <- jsonlite::fromJSON("metadata.cart.2020-03-23.json") #“xxx”里面的内容根据下载的JSON文件改动
colnames(meta)
ids <- meta$associated_entities;class(ids)
ids[[1]]
ids[[1]][,1] #检查tcga ID所处第几列，并改变后面的取值ids[[1]][,？] 

可以看到，meta$associated_entities是个列表，这个列表里包含数据框，数据框的第一列内容就是tcga样本id。

image.png

注意，换了数据需要自己探索存放在哪一列。不一定是完全一样的，需要确认清楚。
批量提取associated_entities（ids）列表元素里面的基因名

ID = sapply(ids, function(x){x[,1]})   #如果前面ids列位发生改变，此处对应修改[,?]
file2id = data.frame(file_name = meta$file_name,
                     ID = ID) 

image.png

image.png

文件名与TCGA样本ID的对应关系已经得到，接下来是将其添加到表达矩阵中，成为行名。需要找到读取文件的顺序，一一对应修改。

head(file2id$file_name,2)
head(count_files,2)
count_files2 = stringr::str_split(count_files,"/",simplify = T)[,2]
table(count_files2 %in% file2id$file_name)
#检查基因表达矩阵文件名与临床信息文件名是否完全一致

文件名格式并不完全一样，截取count_files“/”后面的名字组成新向量

image.png

count_files2（表达矩阵）的顺序就是列名的顺序，根据它来调整file2id（样本）的顺序。此处需要再次理解一下match函数。

file2id = file2id[match(count_files2,file2id$file_name),]
colnames(exp) = file2id$ID
exp[1:4,1:4]

表达矩阵整理完成,需要过滤一下那些在很多样本里表达量都为0的基因。过滤标准不唯一。

dim(exp)
#exp = exp[rowSums(exp)>0,]  
 #筛选条件所有样本A基因表达量总和>0则保留此行
exp = exp[apply(exp, 1, function(x) sum(x > 1) > 9), ]
#对exp这个矩阵1（行）/2（列）执行后面的函数，A基因表达了>1的样本总数>9的行保留，>9这个条件可根据情况改变
dim(exp)
exp[1:4,1:4]
exp = as.matrix(exp)

分组信息

image.png

根据样本ID的第14-15位，给样本分组（tumor和normal）

table(str_sub(colnames(exp),14,15))  
Group = ifelse(as.numeric(str_sub(colnames(exp),14,15)) < 10,'tumor','normal')
Group = factor(Group,levels = c("normal","tumor"))
table(Group)
save(exp,clinical,Group,cancer_type,file = paste0(cancer_type,"_gdc.Rdata"))

三大R包差异分析

if(!require(stringr))install.packages('stringr')
if(!require(ggplotify))install.packages("ggplotify")
if(!require(patchwork))install.packages("patchwork")
if(!require(cowplot))install.packages("cowplot")
if(!require(DESeq2))BiocManager::install('DESeq2')
if(!require(edgeR))BiocManager::install('edgeR')
if(!require(limma))BiocManager::install('limma')

rm(list = ls())
load("TCGA-CHOL_gdc.Rdata") #癌症类型及名字,后面内容无需改变！！！！#全为前面的注释
table(Group)

deseq2----

library(DESeq2)
colData <- data.frame(row.names =colnames(exp), 
                      condition=Group)
if(!file.exists(paste0(cancer_type,"_dd.Rdata"))){  #第一次运行需要将数据转换为deseq2所需要的类型，并保存_dd.Rdata格式，再次运行将跳过此步骤
  dds <- DESeqDataSetFromMatrix(
  countData = exp,
  colData = colData,
  design = ~ condition)
  dds <- DESeq(dds)
  save(dds,file = paste0(cancer_type,"_dd.Rdata"))
}
load(file = paste0(cancer_type,"_dd.Rdata"))

# 两两比较
res <- results(dds, contrast = c("condition",rev(levels(Group))))  #三包都需要提供分组信息，默认因子normol在前，tumor在后，加上rev函数逆转顺序
resOrdered <- res[order(res$pvalue),] # 按照P值排序
DEG <- as.data.frame(resOrdered)  #转换为数据框
DEG = na.omit(DEG)  #当取log及其他操作时，数据太小无法识别则返回NA，去除NA操作
head(DEG)

#添加change列标记基因上调下调
logFC_cutoff <- with(DEG,mean(abs(log2FoldChange)) + 2*sd(abs(log2FoldChange)) )#计算logFC_cutoff 95%置信区间； #with的作用是替代$将列变为局部变量进行操作
#logFC_cutoff <- 2/以2作为判断标准和前面置信区间计算出来的值都可以
k1 = (DEG$pvalue < 0.05)&(DEG$log2FoldChange < -logFC_cutoff)
k2 = (DEG$pvalue < 0.05)&(DEG$log2FoldChange > logFC_cutoff)
DEG$change = ifelse(k1,"DOWN",ifelse(k2,"UP","NOT"))
table(DEG$change)  #给数据框加上上调下调的列
head(DEG)

DESeq2_DEG <- DEG

edgeR----

library(edgeR)

dge <- DGEList(counts=exp,group=Group)  #提供表达数据及分组信息
dge$samples$lib.size <- colSums(dge$counts)
dge <- calcNormFactors(dge) 

design <- model.matrix(~0+Group)
rownames(design)<-colnames(dge)
colnames(design)<-levels(Group)

dge <- estimateGLMCommonDisp(dge, design)
dge <- estimateGLMTrendedDisp(dge, design)
dge <- estimateGLMTagwiseDisp(dge, design)

fit <- glmFit(dge, design)
fit2 <- glmLRT(fit, contrast=c(-1,1)) #之前内容都可看为一个函数，fit2是我们想要的东西

DEG=topTags(fit2, n=nrow(exp)) #赋值给DEG
DEG=as.data.frame(DEG)
logFC_cutoff <- with(DEG,mean(abs(logFC)) + 2*sd(abs(logFC)) )

#logFC_cutoff <- 2  #有需要就更换阈值条件
k1 = (DEG$PValue < 0.05)&(DEG$logFC < -logFC_cutoff)
k2 = (DEG$PValue < 0.05)&(DEG$logFC > logFC_cutoff)
DEG$change = ifelse(k1,"DOWN",ifelse(k2,"UP","NOT"))

head(DEG)
table(DEG$change)
edgeR_DEG <- DEG

limma----

library(limma)

design <- model.matrix(~0+Group)  #0是为了给后面的比较矩阵cont.matrix留下位置
colnames(design)=levels(Group)
rownames(design)=colnames(exp)

dge <- DGEList(counts=exp)
dge <- calcNormFactors(dge)

v <- voom(dge,design, normalize="quantile")
fit <- lmFit(v, design)

constrasts = paste(rev(levels(Group)),collapse = "-")
cont.matrix <- makeContrasts(contrasts=constrasts,levels = design) 
fit2=contrasts.fit(fit,cont.matrix)
fit2=eBayes(fit2)

DEG = topTable(fit2, coef=constrasts, n=Inf)
DEG = na.omit(DEG)

logFC_cutoff <- with(DEG,mean(abs(logFC)) + 2*sd(abs(logFC)) )

#logFC_cutoff <- 2
k1 = (DEG$P.Value < 0.05)&(DEG$logFC < -logFC_cutoff)
k2 = (DEG$P.Value < 0.05)&(DEG$logFC > logFC_cutoff)
DEG$change = ifelse(k1,"DOWN",ifelse(k2,"UP","NOT"))
table(DEG$change)
head(DEG)

检验我们的上下调基因是否比反

exp["ENSG00000042781.11",]  #看到第一行基因不同样本表达量情况
as.numeric(exp["ENSG00000042781.11",])
boxplot(as.numeric(exp["ENSG00000042781.11",])~Group)   #显示下调 #根据基因表达情况画箱线图，对比基因上调下调情况
DEG$logFC[1]  #结果返回值也是负数，表示下调，操作正确

不同方式算出来基因表达数量

limma_voom_DEG <- DEG
tj = data.frame(deseq2 = as.integer(table(DESeq2_DEG$change)),
           edgeR = as.integer(table(edgeR_DEG$change)),
           limma_voom = as.integer(table(limma_voom_DEG$change)),
           row.names = c("down","not","up")
          );tj    注意稳定的基因占大多数才是正常情况
save(DESeq2_DEG,edgeR_DEG,limma_voom_DEG,Group,tj,file = paste0(cancer_type,"_DEG.Rdata"))

如果最后画出来的火山图基因数目比较少，可通过改变阈值扩大范围，logFC_cutoff<- 1或者0.585，只要有出现都要进行更改

当热图上面的图标并未按照分组进行聚类比较时：

draw_heatmap(dat[c(up,down),],Group,n_cutoff = 2,
            cluster_cols = F) 
#cluster_cols = F取消聚类，聚类文件按表达矩阵原文件随机排序
m = dat[c(up,down),]  #挑选上下调基因所在行
m = m[,order(Group)]  #所选基因按照分组信息进行排序
Group = sort(Group) #分组信息随之更换，因为化热图需要
draw_heatmap(m,Group,n_cutoff = 2,
             cluster_cols = F)