ribotish 质控结果复现及重新绘制

测试开头

一个努力中的公众号

1引言

之前有提到过 ribotish 这个软件,用于质控比对到基因组上的 Ribo-seq 的 bam 文件。结果如下:

2安装

github 地址:

https://github.com/zhpn1024/ribotish

$ pip install ribotish

使用:

$ ribotish quality -b chx.bam -g gene.gtf

然后目录下会生成一个 pdf 和一个绘图数据和一个参数文件结果:

pdf 里的图如果放到文章里可能还不太美观,绘图的数据都在 txt 文件里,以字典格式保存,因为软件本身也是基于 python 分析的:

我打红勾的五个分别为左边 5 end matched 的数据,分别为reads 长度分布,距离 start codon,距离 stop codon 条形图,frame 分布及cds profile五个类型的图。

我们尝试拿绘图数据提取出来在 R 里进行自定义绘制这五种图!还可以美化。

下面是复现及拓展结果:

3reads 长度分布

读取数据进来:

library(ggplot2)
library(ggsci)
library(reshape2)
library(tidyverse)
library(cowplot)
library(aplot)
library(stringr)

# load qc-data
qc <- read.delim('ribo-EB.sorted_qual.txt',header = F)[1:5,] %>% as.data.frame()

处理数据转化成数据框:

################################################### reads length

# process raw data
rawlen <- qc[1,1] %>% str_replace(.,pattern = '\{',replacement = '') %>%
  str_replace(.,pattern = '\}',replacement = '') %>%
  str_split(.,pattern = ',') %>% unlist()

# get length and reads number
readlength = sapply(strsplit(rawlen,split = ':'), '[',1) %>% as.numeric()
readsnumber = sapply(strsplit(rawlen,split = ':'), '[',2) %>% as.numeric()

# to dataframe
df_len <- data.frame(readlength,readsnumber)

# check
head(df_len,3)
#   readlength readsnumber
# 1         25      168849
# 2         26      236643
# 3         27      487806

画图:

# plot
p1 <- ggplot(df_len,aes(y = factor(readlength),x = readsnumber/10^6)) +
  geom_col(fill = 'grey40') + theme_classic(base_size = 16) +
  theme(axis.ticks.length = unit(0.3,'cm'),
        # aspect.ratio = 3,
        axis.text.x = element_text(angle = 0),
        axis.title = element_text(size = 16),
        plot.margin = margin(1,1,1,0)) +
  scale_y_discrete(position = 'left') +
  # scale_x_reverse() +
  # ylab('Footprint length (nt)')
  ylab('') +
  xlab('Reads (x 10^6)')

p1