ribo-seq數(shù)據(jù)可以提供不同狀態(tài)下的翻譯效率，本帖子主要介紹2016年發(fā)表在NC上的簡便的工具：Xtail合搅。詳情請參照github網(wǎng)址：https://github.com/xryanglab/xtail
- Z. Xiao, Q. Zou, Y. Liu, X. Yang, Genome-wide assessment of differential translations with ribosome profiling data. Nat Commun 7, 11194 (2016).
背景及簡介
mRNA翻譯的嚴格調(diào)控對于精確控制蛋白質(zhì)的豐度和質(zhì)量至關(guān)重要蒋歌。核糖體分析是一種結(jié)合核糖體足跡測序和RNA深度測序的技術(shù)，已在多種研究中用于定量全基因組范圍的mRNA翻譯汰蓉。在此绷蹲，研究者開發(fā)了Xtail，一種專為核糖體分析數(shù)據(jù)設(shè)計的分析流程顾孽，能夠全面且準確地識別成對比較中差異翻譯的基因祝钢。通過對模擬數(shù)據(jù)和實際數(shù)據(jù)集的應(yīng)用，Xtail展現(xiàn)出高靈敏度和極低的假陽性率若厚，在差異翻譯量化的準確性方面優(yōu)于現(xiàn)有方法太颤。利用已發(fā)布的核糖體分析數(shù)據(jù)集，Xtail不僅揭示了具有生物學(xué)意義的差異翻譯基因盹沈，還發(fā)現(xiàn)了在人類癌細胞mTOR信號擾動和人類初級巨噬細胞干擾素γ（IFN-γ）處理中的新差異翻譯事件龄章。這表明吃谣，Xtail在揭示涉及翻譯失調(diào)的分子機制方面具有提供新見解的價值。
原理+流程
對于特定的mRNA轉(zhuǎn)錄本做裙，RPF（核糖體保護片段）的豐度表明所有該mRNA物種副本上的核糖體占有率岗憋。映射到特定基因編碼區(qū)的RPF counts取決于mRNA豐度及其翻譯速率。因此锚贱，在比較兩種條件時仔戈，差異翻譯可以通過兩種條件下mRNA和RPF表達變化的不同來表征。
Xtail通過三個主要步驟（統(tǒng)計建模拧廊、概率分布建立和差異翻譯評估）定量分析兩種條件下基因的差異翻譯监徘。最終通過兩個并行管道比較RPF和mRNA的log2折疊變化或RPF對mRNA的log2比率，得出差異翻譯的評估結(jié)果吧碾。
安裝及使用
1.安裝
library(devtools)
install_github("xryanglab/xtail")
2.使用
setwd('~/project/HKU/data/20240325/riboseq/xtail/N0/')
library(xtail)
library(tidyverse)
library(ggplot2)
# 讀取riboseq和mrna的reads counts matrix凰盔。
dir_ribo <- '~/project/HKU/data/20240325/riboseq/HTSeq/'
dir_mrna <- '~/project/HKU/data/20240325/RNA/METTL3_KO/HTSeq'
ribo0<- read.table(file.path(dir_ribo,'HTSeq_merged.txt'),header=T,quote='',check.names=F, sep='\t',row.names=1)
mrna0<- read.table(file.path(dir_mrna,'HTSeq_merged.txt'),header=T,quote='',check.names=F, sep='\t',row.names=1)
head(ribo0)
# 去除任何樣本表達都是0的
ribo_f <- ribo0[rowSums(ribo0 == 0) == 0, ]
mrna_f <- mrna0[rowSums(mrna0 == 0) == 0, ]
GENE <- intersect(rownames(ribo_f), rownames(mrna_f)) # 剩余交集基因：8722
# 提取目標基因行 ；PS：如果不去除在樣本中表達為0的基因就會報錯
ribo <- ribo_f[GENE,]
mrna <- mrna_f[GENE,]
# 樣本信息
condition <- c("control","control","treat","treat") # 每個樣本對應(yīng)的分組
# 最關(guān)鍵的一步倦春，xtail按照q-value排序户敬，選擇輸出log2FCs和log2R
results <- xtail(mrna,ribo,condition,baseLevel = "control",threads = 55,minMeanCount=1,bins=10000)
# log2FC: 2366
# log2R: 12500
# 輸出xtail分析結(jié)果
saveRDS(results,"./files/ribo_EPSC.rds")
results_tab <-resultsTable(results,sort.by="pvalue.adjust",log2FCs=TRUE,log2Rs=TRUE)
write.table(results_tab,"./files/ribo_EPSC.csv",quote=F,sep=",")
得到的結(jié)果如下所示，行是基因睁本，列是結(jié)果尿庐，有具體的mRNA和RPF(對應(yīng)ribo-seq)的差異結(jié)果，最后三列給出了最終確定的差異翻譯基因呢堰，可以后續(xù)進行分析抄瑟、繪圖。
3.可視化
results <- readRDS("./files/ribo_EPSC.rds")
results_tab <- as.data.frame(resultsTable(results,sort.by="pvalue.adjust",log2FCs=TRUE,log2Rs=TRUE))
# 轉(zhuǎn)錄組/核糖體組枉疼、翻譯組總體情況
pdf("./plots/ribo_RNA_EPSC.pdf", width = 5, height = 3)
plotFCs(results) # mRNA差異散點圖
plotRs(results) # 差異翻譯散點圖
volcanoPlot(results) # 差異翻譯火山圖
dev.off()
4. 功能富集
#### 定義上下調(diào)基因 ---
library(tidyverse)
library(clusterProfiler)
library(org.Hs.eg.db)
colnames(results_tab) #"log2FC_TE_final" "pvalue_final" "pvalue.adjust"
results_tab <- results_tab %>% mutate(change = ifelse(
pvalue_final <= 0.05 & abs(log2FC_TE_final) >= 1,
ifelse(log2FC_TE_final > 1, "Up", "Down"),
"Stable" ))
results_tab %>% dplyr::count(change)
# change n
# 1 Down 462
# 2 Stable 7766
# 3 Up 494
# 名稱轉(zhuǎn)換為易懂的SYMBOL格式
results_tab$ENSEMBL <- rownames(results_tab)
results_tab$ENSEMBL <- sub("\\..*", "", results_tab$ENSEMBL)
ref <- bitr(results_tab$ENSEMBL, fromType="ENSEMBL", toType=c("SYMBOL"), OrgDb="org.Hs.eg.db")
results_tab <- left_join(results_tab, ref, by = c("ENSEMBL" = "ENSEMBL"))
write.table(results_tab,"./files/ribo_EPSC_2.csv",quote=F,sep=",")
#### 功能富集 ----
outdir <- './files/FUN'
outdir2 <- './plots/FUN'
for(i in c("Up", "Down")){
s_gene <- results_tab |> filter(change %in% i) |> pull('SYMBOL')
bp <-
enrichGO(
s_gene,
OrgDb = org.Hs.eg.db,
keyType = 'SYMBOL',
ont = "BP",
pAdjustMethod = "BH",
pvalueCutoff = 0.05,
qvalueCutoff = 0.2)
term <- bp@result
write.table(term,
file = file.path(outdir , paste0(i,'_EPSC_terms.csv')),
quote = F,sep = ",", row.names = F)
#繪制TOP30
df <- term[1:30,]
df$labelx=rep(0,nrow(df))
df$labely=seq(nrow(df),1)
ggplot(data = df,
aes(x = -log10(pvalue),
y = reorder(Description,-log10(pvalue)))) +
geom_bar(stat="identity",
alpha=1,
fill= "#BEBADA",
width = 0.8) +
geom_text(aes(x=labelx,
y=labely,
label = df$Description),
size=3.5,
hjust =0)+
theme_classic()+
theme(axis.text.y = element_blank(),
axis.line.y = element_blank(),
axis.title.y = element_blank(),
axis.ticks.y = element_blank(),
axis.line.x = element_line(colour = 'black', linewidth = 1),
axis.text.x = element_text(colour = 'black', size = 10),
axis.ticks.x = element_line(colour = 'black', linewidth = 1),
axis.title.x = element_text(colour = 'black', size = 12))+
xlab("-log10(pvalue)")+
ggtitle(i)+
scale_x_continuous(expand = c(0,0))
ggsave(file.path(outdir2, paste0(i, '_EPSC_Barplot_TOP30.pdf')),
ggplot2::last_plot(),#最后一張圖
height=6,
width=6)
}
根據(jù)上調(diào)的差異翻譯基因可能發(fā)現(xiàn)血管運輸皮假、物質(zhì)轉(zhuǎn)運、細胞黏附等過程在實驗組中上調(diào)往衷；而胞質(zhì)翻譯钞翔、物質(zhì)合成严卖、RNA剪切等生物學(xué)過程在實驗組中下調(diào)席舍。
歡迎大家評論交流！
（每帖分享：人生一直游蕩于“自信巔峰”到“絕望之谷”的動態(tài)過程中）
