內(nèi)容獲取 1、購(gòu)買打包合集（《KS科研分享與服務(wù)》付費(fèi)內(nèi)容打包集合），價(jià)格感人，可以加入微信VIP群（答疑交流群，甚至有小伙伴覺(jué)得群比代碼更好），可以獲取建號(hào)以來(lái)所有內(nèi)容，群成員專享視頻教程，提前更新，其他更多福利！ 2、《KS科研分享與服務(wù)》公眾號(hào)有QQ群，進(jìn)入門檻是20元（完全是為了防止白嫖黨，請(qǐng)理解），請(qǐng)考慮清楚。群里有免費(fèi)推文的注釋代碼和示例數(shù)據(jù)（終身?yè)碛校瑳](méi)有付費(fèi)內(nèi)容，群成員福利是購(gòu)買單個(gè)付費(fèi)內(nèi)容半價(jià)！需要者詳情請(qǐng)聯(lián)系作者（非需要者勿擾，處理太費(fèi)時(shí)間）：
我們之前介紹了三元圖用于轉(zhuǎn)錄組數(shù)據(jù)三組的展示（Ternary plot：三元圖及其在轉(zhuǎn)錄組中的部分應(yīng)用示例），也復(fù)現(xiàn)過(guò)重疊的火山圖用于差異基因展示（復(fù)現(xiàn)nature圖表：雙重組合-漸變背景火山圖-效果？無(wú)需多言杠杠的！）。小伙伴也提出有另外的包Triwise實(shí)現(xiàn)三組數(shù)據(jù)的展示。所以我們嘗試展示一下，對(duì)于bulk和scRNA都可以，這個(gè)包還可以展示通路基因集，如果你的數(shù)據(jù)是三組，或者轉(zhuǎn)錄組你要比較三種celltype結(jié)果展示，那么可以選擇這種方式！Triwise包出現(xiàn)很早，很多高分文章也使用過(guò)它的圖，官網(wǎng)有詳細(xì)的示例教程，我們這里按照實(shí)際情況，實(shí)際的數(shù)據(jù)進(jìn)行演示，也進(jìn)行了拓展。示例數(shù)據(jù)及注釋代碼已發(fā)布群里！ GitHub鏈接：https://github.com/saeyslab/triwise
paper：Yolk Sac Macrophages, Fetal Liver, and Adult Monocytes Can Colonize an Empty Niche and Develop into Functional Tissue-Resident Macrophages. van de Laar et al. Immunity, Volume 44, Issue 4, 755 - 768devtools::install_github("saeyslab/triwise")library(triwise)library(limma)library(Biobase)library(ggplot2)library(ggrepel)
setwd("D:\\KS項(xiàng)目\\公眾號(hào)文章\\Triwise轉(zhuǎn)錄組三組結(jié)果展示")
df <- read.csv('Bulk_FPKM.csv', header = T, row.names = 1)avg_exp <- data.frame(E1 = matrixStats::rowMeans2(as.matrix(df[ 1:3])), E18 = matrixStats::rowMeans2(as.matrix(df[ 4:6])), E30 = matrixStats::rowMeans2(as.matrix(df[ 7:9])))barycoords = transformBarycentric(avg_exp)#將包含三種生物狀態(tài)的表達(dá)式矩陣轉(zhuǎn)換為質(zhì)心坐標(biāo)plotDotplot(barycoords)
design = as.data.frame(colnames(df))rownames(design) <- colnames(df)design$sample <- c(rep("E1",3),rep("E18",3),rep("E30",3))design$sample = factor(as.character(design$sample))design <- model.matrix(~0+sample, design)
fit <- lmFit(df, design)fit = contrasts.fit(fit, matrix(c(1, -1, 0, 0, 1, -1, 1, 0, -1), ncol=3))fit = eBayes(fit)
top = topTable(fit, p.value=0.05, lfc=1, number=Inf, sort.by = "none")Gdiffexp = rownames(top)
#triwise作圖可以與ggplot對(duì)接，所以可以用ggplot主題做一些修飾plotDotplot(barycoords, Gdiffexp)+ scale_color_manual(values = c("#B40F20","#C1AE8D"))
plotRoseplot(barycoords, Gdiffexp, relative = F)
也可以是你需要展示的每組代表性的關(guān)鍵基因，Gdiffexp參數(shù)傳入基因向量即可。clipHexagon <- function(barypoints, rmax=5, baseangle=0) { barypoints["rclip"] = mapply(function(angle, r) {min(hexagonPolar(angle, rmax, baseangle), r)}, barypoints$angle, barypoints$r) barypoints["xclip"] = cos(barypoints["angle"]+baseangle) * barypoints["rclip"] barypoints["yclip"] = sin(barypoints["angle"]+baseangle) * barypoints["rclip"] barypoints}
#需要標(biāo)記的基因E18_gene <- c("DHFR2","TMPRSS15","BBS10","TRADD","ZNF35")E18_gene = barycoords[E18_gene,]E18_gene = clipHexagon(E18_gene)E18_gene$label = rownames(E18_gene)
E30_gene <- c("PSD3","TNPO1","RBBP9","DGCR10","HELZ2","SOX6","OR1L8","PAQR7")E30_gene = barycoords[E30_gene,]E30_gene = clipHexagon(E30_gene)E30_gene$label = rownames(E30_gene)
E1_gene <- c("OR2S2","LRRC61","SOHLH2","MPL","CUL5", "TAMM41","MIR652","STX10","RAD51C")E1_gene = barycoords[E1_gene,]E1_gene = clipHexagon(E1_gene)E1_gene$label = rownames(E1_gene)
plotDotplot(barycoords, Gdiffexp)+ scale_color_manual(values = c("#eb6841","#C1AE8D"))+ geom_point(data = E1_gene, aes(x=xclip, y=yclip), size=3, color="#053061")+ geom_text_repel(data = E1_gene, aes(x=xclip, y=yclip, label=label), nudge_y = 1.5, nudge_x = 1, color = 'black', size = 3,fontface = 'italic', min.segment.length = 0, max.overlaps = Inf, box.padding = unit(0.5, 'mm'), point.padding = unit(0.5, 'mm'))+ geom_point(data = E30_gene, aes(x=xclip, y=yclip), size=3, color='#24796C')+ geom_text_repel(data = E30_gene, aes(x=xclip, y=yclip, label=label), nudge_y = -1, nudge_x = -1, color = 'black', size = 3,fontface = 'italic', min.segment.length = 0, max.overlaps = Inf, box.padding = unit(0.5, 'mm'), point.padding = unit(0.5, 'mm'))+ geom_point(data = E18_gene, aes(x=xclip, y=yclip), size=3, color='#FA8E24')+ geom_text_repel(data = E18_gene, aes(x=xclip, y=yclip, label=label), nudge_y = 0, nudge_x = -2, color = 'black', size = 3,fontface = 'italic', min.segment.length = 0, max.overlaps = Inf, box.padding = unit(0.5, 'mm'), point.padding = unit(0.5, 'mm'))
 #上述的展示在scRNA中依然可以使用，剛好我有一個(gè)示例數(shù)據(jù)是三個(gè)分組#也可以展示三個(gè)亞群
uterus <- readRDS("D:/KS項(xiàng)目/公眾號(hào)文章/uterus.rds")#一般情況肯定是比較某一個(gè)celltype中三組基因的差異，這里用全部細(xì)胞展示gene_cell_exp <- AverageExpression(uterus, features = rownames(uterus), group.by = 'orig.ident', slot = 'data')
gene_cell_exp <- as.data.frame(gene_cell_exp$RNA)
sce_barycoords = transformBarycentric(gene_cell_exp)# plotDotplot(sce_barycoords)
#用FindAllMarkers計(jì)算每組差異基因，至于要怎么label就是自己的事了Idents(uterus) <- 'orig.ident'gene_marker <- FindAllMarkers(uterus, logfc.threshold = 0.5, min.pct = 0.25, only.pos = T)
plotDotplot(sce_barycoords, gene_marker$gene)+ scale_color_manual(values = c('#B497E7',"darkgrey"))
除了標(biāo)記感興趣的部分基因，也可以標(biāo)記基因集。#我這里示例兩個(gè)通路,直接讀入整理好的通路基因集gene_set <- read.csv('gene_set.csv', header = T, check.names = F)gsets <- list(gene_set$`Glycolysis/Gluconeogenesis`[!grepl("^$", gene_set$`Glycolysis/Gluconeogenesis`)], gene_set$`Oxidative phosphorylation`)
names(gsets) <- colnames(gene_set)gsets = sapply(gsets, function(gset) intersect(rownames(gene_cell_exp), unique(as.character(gset))))
scores = testUnidirectionality(sce_barycoords, gsets, gene_marker$gene, statistic = "rank", mc.cores = 1, nsamples=1e+6)
library(cowplot)cowplot::plot_grid( plotDotplot(sce_barycoords, gene_marker$gene, Goi=gsets[["Oxidative phosphorylation"]], showlabels = F) + ggplot2::theme(legend.position = "none")+ ggplot2::ggtitle("Oxidative phosphorylation"), plotRoseplot(sce_barycoords, gene_marker$gene, Goi=gsets[["Oxidative phosphorylation"]], relative = T, showlabels = F))
覺(jué)得我們分享有些用的，點(diǎn)個(gè)贊再走唄！
|
