其他
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
微信扫一扫付费阅读本文
可试读18%
微信扫一扫付费阅读本文