TCGA-19-临床数据整理及人群分组
代码解释：
# 07 miRNA差异表达
rm(list = ls())
# miRNA数据下载及整理
library(TCGAbiolinks)
library(limma)
# 查看有哪些分类数据 project="TCGA-CHOL"
getProjectSummary("TCGA-CHOL")
query <- GDCquery(project = 'TCGA-CHOL',
data.category = "Transcriptome Profiling",
data.type = "miRNA Expression Quantification",
experimental.strategy = "miRNA-Seq",
workflow.type = "BCGSC miRNA Profiling")
GDCdownload(query)
expData<- GDCprepare(query = query)
#读取下载的query，用GDCprepare函数以SummarizedExperiment形式读入rstudio
expData[1:4,1:4] # counts第二列，rpm第三列
#这个代码的用处就是提取前4行和前4列展示出来，没用
dim(expData)
#查看维度，没用的代码
expr <- expData[,seq(2,136,3)]
#看OneNote解释
dim(expr)
colnames(expr) <- substring(colnames(expr),12)
# substr()需要指定stop值，substring()不需要
#提取列标题的12位及以后的字符串作为标题
expr <- data.frame(expData$miRNA_ID,expr)
#expData中提取出来的miRNA_ID列与expr数据框合并，重新命名为expr
names(expr)[1] <- 'id'
#把exper列名命名为id
#这样就得到未处理的counts矩阵
# 筛选
expr=as.matrix(expr)
#强制转为矩阵
rownames(expr)=expr[,1]
#第一列命名为行名
expr.matrix=expr[,2:ncol(expr)]
#提取除第一列以外的所有列，命名为expr.matrix
dimnames=list(rownames(expr.matrix),colnames(expr.matrix))
#定义一个列表 名为dimnames
data=matrix(as.numeric(as.matrix(expr.matrix)),nrow=nrow(expr.matrix),dimnames=dimnames)
#定义一个矩阵，名为data，强制转为数值型矩阵
data=avereps(data)
#avereps()函数：去重复：对于行上出现的相同基因，取平均
data=data[rowMeans(data)>0,]
#去除低表达基因，同一行上数据的平均值大于零的留下
dim(data)
data1 <- rbind(colnames(data),data)
#把列名作为第一行，命名为data1
write.table(data1, file="CHOL-miRNA-matrix.txt", sep="\t", quote=F, col.names=F)
#导出data1位txt文件，名称为CHOL-miRNA-matrix
rm(list = ls())
# 导入数据
rt <- read.table("CHOL-miRNA-matrix.txt",sep="\t",header=T,check.names=F)
#读取文件数据入R，check.names=TRUE检查变量名在R中是否有效
#数据处理
rownames(rt) <- rt[,1]
#把第一列作为行名
rt <- rt[,-1]
#删除第一列
dim(rt)
#返回维度信息
# 用substr函数在TCGA数据中提取样本信息
group <- factor(ifelse(as.integer(substr(colnames(rt),14,15))<10,'tumor','normal'),
levels = c('normal','tumor'))
#提取第14和15位，用来判断是正常组还是肿瘤组
table(group)
#table()函数——向量内元素频数统计，本例题是提取group中的频数统计
# 差异性分析
library(edgeR)
logFC_cutoff=1
padj=0.05
#定义两个变量：logFC_cutoff=1判断差异有显著差异的依据， padj是调整后的p值，
dge <- DGEList(counts=rt,group=group) #将数据打包为edgeR包识别的数据 group=group表示分类的因子变量是刚刚定义的分组变量
dge$samples$lib.size <- colSums(dge$counts) #colsums()函数——计算列求和 lib.size一个向量，表示每个样本的总计数。
dge <- calcNormFactors(dge) #标准化，消除量纲的影响
design <- model.matrix(~0+factor(group)) #看onenote解释
colnames(design) <- c('normal','tumor')
rownames(design)<-colnames(dge)
colnames(design)<-levels(group)
#这一段代码的目的看onenote
dge <- estimateGLMCommonDisp(dge,design)
dge <- estimateGLMTrendedDisp(dge, design)
dge <- estimateGLMTagwiseDisp(dge, design)
fit <- glmFit(dge, design)
fit1 <- glmLRT(fit, contrast=c(-1,1))
DEG=topTags(fit1, n=nrow(rt))
DEG=as.data.frame(DEG)
# 输出有显著差异的结果
significant = DEG[(
(DEG$FDR < padj) & (abs(DEG$logFC) > logFC_cutoff)
),] # 筛选显著差异的基因
write.table(significant, file="significant.xls",sep="\t",quote=F)
# 输出表达上调的结果
upregulation = DEG[(DEG$FDR < padj & (DEG$logFC>logFC_cutoff)),]
write.table(upregulation, file="upregulation.xls",sep="\t",quote=F)
# 输出表达下调的结果
downregulation = DEG[(DEG$FDR < padj & (DEG$logFC<(-logFC_cutoff))),]
write.table(downregulation, file="downregulation.xls",sep="\t",quote=F)
# 添加显著性标签
DEG$change <- factor(ifelse(DEG$PValue < padj & abs(DEG$logFC) > logFC_cutoff,
ifelse(DEG$logFC > logFC_cutoff ,'UP','DOWN'),
'NOT'))
#这个代码的目的是再DEG数据框中添加名称位change的列，用于显示基因表达是上调了下降了还是不变
# 输出结果
write.csv(DEG,file="miRNA-edgeR.csv")
# 结果可视化：火山图 and 热图
# ggplot2绘制火山图 volcano
library(ggplot2)
g <- ggplot(DEG,aes(logFC,-log10(as.numeric(FDR)),color=change))
g + geom_point()
g + geom_point(size = 1) +
labs(x = 'log2 Fold Change', y = '-log10 adjust p-value', title = 'volcano_plot') +
theme(plot.title = element_text(hjust = 0.5, size = 20,color = 'firebrick'),
panel.grid = element_blank(),
panel.background = element_rect(color = 'black', fill = 'transparent'),
legend.key = element_rect(fill = 'transparent')) +
geom_vline(xintercept = c(-logFC_cutoff, logFC_cutoff), lty = 2, color = 'black') + #添加阈值线
geom_hline(yintercept = -log10(padj), lty = 2, color = 'black')
library(ggrepel)# ggrepel包给点添加点标签的包
significant = subset(DEG,abs(logFC)>3 & -log10(as.numeric(FDR)) > 5)
#生成significant变量，里面的值决定了到底给哪些点添加标签
#geom_text_repel函数才是给点添加标签的函数
g + geom_point(size = 1) +
labs(x = 'log2 Fold Change', y = '-log10 adjust p-value', title = 'volcano_plot') +
theme(plot.title = element_text(hjust = 0.5, size = 20,color = 'firebrick'),
panel.grid = element_blank(),
panel.background = element_rect(color = 'black', fill = 'transparent'),
legend.key = element_rect(fill = 'transparent')) +
geom_vline(xintercept = c(-logFC_cutoff, logFC_cutoff), lty = 2, color = 'black') + #添加阈值线
geom_hline(yintercept = -log10(padj), lty = 2, color = 'black') +
geom_text_repel(data = significant,aes(label = rownames(significant)),
max.overlaps = 10000, # 最大覆盖率，当点很多时，有些标记会被覆盖，调大该值则不被覆盖
size=3,col = 'black')
#视频解释：【TCGA-09-miRNA差异表达】 【精准空降到 11:07】 https://www.bilibili.com/video/BV1Zv4y1J7D7/?share_source=copy_web&;vd_source=3d598f1156bea5867cef7412e80ded78&t=667
# 绘制热图heatmap
# df <- tibble::column_to_rownames(df,var = 'id')
pheatmap::pheatmap(rt,
#annotation_row=dfGene, # （可选）指定行分组文件
#annotation_col=dfSample, # （可选）指定列分组文件
show_colnames = TRUE, # 是否显示列名
show_rownames=TRUE, # 是否显示行名
fontsize=1, # 字体大小
color = colorRampPalette(c('#0000ff','#ffffff','#ff0000'))(50), # 指定热图的颜色
#annotation_legend=TRUE, # 是否显示图例
border_color=NA, # 边框颜色 NA表示没有
scale="row", # 指定归一化的方式。"row"按行归一化，"column"按列归一化，"none"不处理
cluster_rows = TRUE, # 是否对行聚类
cluster_cols = TRUE # 是否对列聚类
)
rm(list = ls())
