写在前面：本文为微信公众号:生信星球的数据挖掘线上班的随堂笔记，感谢小洁老师的付出！
作业
今天先写作业，不然往后拖又不写作业了。
1. idmap annoGene geoChina
- 对基因名进行注释/转换等的函数，出自jmzeng的annoprobe包。
如果前期没有进行对表达矩阵进行log处理，则图如下。
3.1 添加Entrezid列。
- ENTREID适合做富集分析，因此需要转换↓。
- 记得改种属！！
#4.【种属！】加ENTREZID列，用于富集分析（symbol转entrezid，然后inner_join）
library(ggplot2)
library(clusterProfiler)
library(org.Hs.eg.db)
s2e <- bitr(deg$symbol,
fromType = "SYMBOL",
toType = "ENTREZID",
OrgDb = org.Hs.eg.db)#人类!!如果不是人一定要改
#其他物种http://bioconductor.org/packages/release/BiocViews.html#___OrgDb
deg <- inner_join(deg,s2e,by=c("symbol"="SYMBOL"))
save(group_list,deg,logFC_t,P.Value_t,file = "step4output.Rdata")
4.绘制差异基因火山图/热图
- 复杂的热图：ComplexHeatmap包
rm(list = ls())
load(file = "step1output.Rdata")
load(file = "step4output.Rdata")
#1.火山图----
library(dplyr)
library(ggplot2)
dat = deg
for_label <- dat%>%
filter(abs(logFC)>logFC_t & P.Value < P.Value_t) %>%
head(10)#取差异基因的前十个标上名字
p <- ggplot(data = dat,
aes(x = logFC,
y = -log10(P.Value))) +
geom_point(alpha=0.4, size=3.5,
aes(color=change)) +
ylab("-log10(Pvalue)")+
scale_color_manual(values=c("blue", "grey","red"))+
geom_vline(xintercept=c(-logFC_t,logFC_t),lty=4,col="black",lwd=0.8) +
geom_hline(yintercept = -log10(P.Value_t),lty=4,col="black",lwd=0.8) +
theme_bw()
p
#添加图片上的文字信息
volcano_plot <- p +
geom_point(size = 3, shape = 1, data = for_label) +
ggrepel::geom_label_repel(
aes(label = symbol),
data = for_label,
color="black"
)
volcano_plot
ggsave(plot = volcano_plot,filename = paste0(gse,"volcano.png"))
#2.差异基因热图----
load(file = 'step2output.Rdata')
x=deg$logFC
names(x)=deg$probe_id
#以下那行是根据logfc选取上下调30名的差异基因热图
#cg=c(names(head(sort(x),30)),names(tail(sort(x),30)))
#以下为对所有差异基因做热图
cg = deg$probe_id[deg$change !="stable"]#先根据列值$选出逻辑值，然后用列$来取
n=exp[cg,]
dim(n)
#作热图
library(pheatmap)
annotation_col=data.frame(group=group_list)
rownames(annotation_col)=colnames(n)
library(ggplotify)
heatmap_plot <- as.ggplot(pheatmap(n,show_colnames =F,
show_rownames = F,
scale = "row",
#cluster_cols = F,
annotation_col=annotation_col))
heatmap_plot
ggsave(heatmap_plot,filename = paste0(gse,"heatmap.png"))
load("pca_plot.Rdata")
library(patchwork)
(pca_plot + volcano_plot +heatmap_plot)+ plot_annotation(tag_levels = "A")
- 挑选几个指定基因进行标注_更改for_label参数：
filter（symbol %in% c（“基因1”，“基因2”，“基因3”。。。））
#代码----可以折叠
-
browseVignettes（“clusterProfiler”）找到Y叔写的关于GSEA的包包说明。
5. 绘制富集分析图
复杂数据分析
1. 配对数据：
-
如GSE5109
差异分析代码：
-
绘制热图：有序因子
热图 -
使用ggpubr绘制差异基因的箱线图（需定义group）
2. 多分组数据
- 策略1 ：选出一个对照，其他分组分别与对照进行差异分析（需要向量/矩阵取子集）
- 策略2： 两两对比
- 策略3：wgca分析更多分组
3. 多个series联合
考虑批次效应：
1）选择来自同一芯片平台的series
2）需要处理批次效应Batch effect
蛋白互作
- 生成输入string的数据。
#制作string的输入数据
load("step4output.Rdata")
gene_up= deg[deg$change == 'up','symbol']
gene_down=deg[deg$change == 'down','symbol']
write.table(gene_up,
file="upgene.txt",
row.names = F,
col.names = F,
quote = F)
write.table(gene_down,
file="downgene.txt",
row.names = F,
col.names = F,
quote = F)
write.table(deg$symbol[1:200],
file="diffgene.txt",
row.names = F,
col.names = F,
quote = F)
- STRING网站：使用gene symbol。
- 在string生成的文件里下载tsv格式的文件。
3.导入R - 在R中导出绘制网络图的两个数据。
- 使用cytoscape重新绘制网络图。
- 寻找hub核心基因：插件cytohHubba。
- 寻找子网络：插件MCODE
#准备cytoscape的输入文件,把差异基因拿出来做蛋白质互作网络图
tsv = read.table("string_interactionsdiff.tsv",comment.char = "!",header = T)
tsv2 = tsv[,c(1,2,ncol(tsv))]
head(tsv2)
write.table(tsv2,
file = "cyto.txt",
sep = "\t",
quote = F,
row.names = F)
p = deg[deg$change != "stable",c("symbol","logFC","P.Value")]
head(p)
write.table(p,
file = "deg.txt",
sep = "\t",
quote = F,
row.names = F)
