基因组可视化细节优化

引言

omicScope 最近我自己用的时候也是挺顺手的，我需要绘制几个基因的基因组的可视化表达情况，使用 coverage_plot 绘制的时候发现了一些小的问题，顺手优化了一下。赞赞点起来。

安装

# Install from GitHub
if (!requireNamespace("devtools", quietly = TRUE)) {
    install.packages("devtools")
}
devtools::install_github("junjunlab/omicScope")

优化

加载：

library(omicScope)

gtf <- rtracklayer::import("../../index-data/Mus_musculus.GRCm38.102.gtf.gz")

show_utr 参数可以展示蛋白编码基因的 utr 和 cds 区域，如果你画多个基因的时候，其中有一个非编码基因可能就会没有基因结构，比如 Malat1 这个非编码 RNA：

coverage_plot(bam_file <-  c("../5.map-data/DMSOR1Input.bam",
                             "../5.map-data/DMSOR2Input.bam",
                             "../5.map-data/DM1input.bam",
                             "../5.map-data/DM2input.bam"),
              sample_name = c("DMSO1-rep1", "DMSO1-rep2",
                              "DMSO2-rep1", "DMSO2-rep2"),
              gtf_file = gtf,
              target_gene = c("Nanog","Malat1"),
              collapse_exon = FALSE,
              show_utr = TRUE)

可以判断一下输入的基因是不是蛋白编码基因，如果有非编码的基因，直接拿 exon 来进行替代：

if(show_utr == TRUE){
            strc <- subset(gtf, type %in% c("CDS","five_prime_utr", "three_prime_utr","5UTR","3UTR"))

            # get non-coding strc info
            gin2 <- target_gene %in% unique(strc$gene_name)
            if(all(gin)){
                if(!all(gin2)){
                    strc.nc <- exon |>
                        dplyr::filter(gene_name %in% target_gene[!gin2]) |>
                        dplyr::mutate(type = "CDS")

                    # merge
                    strc <- rbind(strc,strc.nc)
                }
            }
        }else{
            strc <- exon
        }
...

再画一个看看，此外这里给基因上面添加了基因名称的标签：

如果你想知道每个转录本名称是什么，可以设置 gene_label_aes 参数：

coverage_plot(bam_file <-  c("../5.map-data/DMSOR1Input.bam",
                             "../5.map-data/DMSOR2Input.bam",
                             "../5.map-data/DM1input.bam",
                             "../5.map-data/DM2input.bam"),
              sample_name = c("DMSO1-rep1", "DMSO1-rep2",
                              "DMSO2-rep1", "DMSO2-rep2"),
              gtf_file = gtf,
              target_gene = c("Nanog","Malat1"),
              collapse_exon = FALSE,
              show_utr = TRUE,
              gene_label_aes = "transcript_id")

指定区间绘图的时候，转录本会顶在最上面，原因是 y 轴位置没有设置好：

coverage_plot(bam_file <-  c("../5.map-data/DMSOR1Input.bam",
                             "../5.map-data/DMSOR2Input.bam",
                             "../5.map-data/DM1input.bam",
                             "../5.map-data/DM2input.bam"),
              sample_name = c("DMSO1-rep1", "DMSO1-rep2",
                              "DMSO2-rep1", "DMSO2-rep2"),
              gtf_file = gtf,
              target_region = c("5:76176615-76376575",
                                "1:91236162-91493591"),
              collapse_exon = FALSE,
              show_utr = TRUE)

把上下间距控制好：

还有一个就是在折叠转录本显示的时候，标签也会折叠到堆积到一起：

coverage_plot(bam_file <-  c("../5.map-data/DMSOR1Input.bam",
                             "../5.map-data/DMSOR2Input.bam",
                             "../5.map-data/DM1input.bam",
                             "../5.map-data/DM2input.bam"),
              sample_name = c("DMSO1-rep1", "DMSO1-rep2",
                              "DMSO2-rep1", "DMSO2-rep2"),
              gtf_file = gtf,
              target_region = c("5:76176615-76376575",
                                "1:91236162-91493591"),
              collapse_exon = TRUE,
              show_utr = TRUE,
              exon_linewidth = 8)

可以选择只显示一个就行：

你重新安装一下最新版本的 omicScope，这些优化的问题都可以解决。演示的问题是使用的老的版本。对应的单细胞版本的 scCoverage_plot 也同样进行了优化。

结尾

路漫漫其修远兮,吾将上下而求索。

欢迎加入生信交流群。加我微信我也拉你进 微信群聊老俊俊生信交流群(微信交流群需收取 20 元入群费用,一旦交费,拒不退还!(防止骗子和便于管理)) 。