引言

❝
polarity score 可以鉴别核糖体在基因上的翻译差别。下面是计算原则：

原文：

相关图：

❝
我来解释一下，pi 代表基因每个位置的 polarity score，pi 由 di（每个位置的 ribosome density）和 wi（每个位置相对于基因中部的相对位置）的乘积再比上所有位置处的 di（ribosome density）。这样我们就得到了每个位置处的 pi，进行加和就能得到该基因整体的 polarity score。知道计算原理了我们就可以自己去算，当然位置里包含一些筛选规则。
示例代码，我们先算基因总的 ribosome density，然后再算 wi（相对位置），最后算 pi 进行加和即可：
ps_df <- dt |>
dplyr::left_join(y = total_density,by = “tid”,multiple = “all”) |>
dplyr::mutate(abs_pos = pos – utr5,
wi = (2abs_pos -(cds + 1))/(cds – 1), pi = densitywi/sum_density) |>
dplyr::group_by(gene_name) |>
dplyr::summarise(sum_pi = sum(pi)) |>
dplyr::mutate(group = group,sample = sam_file)
安装

install.packages(“devtools”)

devtools::install_github(“junjunlab/RiboProfiler”)

or

remotes::install_github(“junjunlab/RiboProfiler”)

library(RiboProfiler)
测试

这里我们利用 RiboProfiler 进行上游处理，获得每个位置处的 ribosome density，然后根据上面分析定义一个计算 polarity score 的函数 calculatePolarity，下面是相关的参数，包括筛选基因总的 count 数，基因前后 15 nt 处位置的筛选：
calculatePolarity <- function(gene_anno_file = NULL,
gtf_file = NULL,
sam_file = NULL,
density_file = NULL,
minCounts = 64,
upstreamNt = 15,
downstreamNt = 15,
group = NULL)
我们拿数据来试试：
library(ggplot2)
library(RiboProfiler)

sam_file = list.files(“../5.map-data”,pattern = “Ribo.*sam$”,full.names = TRUE)
density_file <- list.files(“2.density-data/”,pattern = “Ribo.trans.txt$”,full.names = TRUE)

df <- calculatePolarity(gene_anno_file = “longest_info.txt”,
gtf_file = “../../0.index-data/Homo_sapiens.GRCh38.110.gtf”,
sam_file = sam_file[1:2],
density_file = density_file[1:2],
group = c(“sample1″,”sample2”))

check

head(df,3)

# A tibble: 3 × 4

gene_name sum_pi group sample

1 AAMP 0.108 sample1 ../5.map-data/1Ribo.sam

2 AARS1 0.0397 sample1 ../5.map-data/1Ribo.sam

3 ABCC1 0.0512 sample1 ../5.map-data/1Ribo.sam

输出结果就是每个样本每个基因的总的 pi 值。这里提供了 polarity_score_plot 来对结果进行绘图：
polarity_score_plot(polarity_score_df = df)

分面：
polarity_score_plot(polarity_score_df = df,facet = T)

Ribo QC

❝
增加了一个热图的可视化函数 footprint_heatmap 来展示 QC 结果，文献里长下面这样：

我们准备一下 qc 的数据：

run

pre_qc_data(longest_trans_file = “longest_info.txt”,
sam_file = sam_file,
out_file = paste(sample_name,”.qc.txt”,sep = “”),
mapping_type = “genome”,
seq_type = “pairedEnd”)

load qc data

qc_df <- load_qc_data()

qc_df$length <- factor(qc_df$length)

check

head(qc_df,3)

length framest relst framesp relsp feature counts sample group

1 29 0 216 0 -444 3 1 1Ribo NA

2 32 2 62 1 -1384 3 1 1Ribo NA

3 32 0 186 0 -807 3 1 1Ribo NA

绘图：
footprint_heatmap(qc_data = qc_df)

绘制离终止密码子处的：
footprint_heatmap(qc_data = qc_df,
type = “relsp”)

结尾

❝
路漫漫其修远兮,吾将上下而求索。
欢迎加入生信交流群。加我微信我也拉你进 微信群聊 老俊俊生信交流群 (微信交流群需收取 20 元入群费用,一旦交费,拒不退还!(防止骗子和便于管理)) 。QQ 群可免费加入, 记得进群按格式修改备注哦。