RNAseqQC：一个用于RNA-seq数据质控评估的R包

RNAseqQC通过一个提供一系列数据可视化功能来帮助RNAseq数据质量控制的R包，它允许识别具有不需要的生物或技术影响的样品，并探索差异测试结果。

源码：https://github.com/frederikziebell/RNAseqQC

安装指南

安装过程简单明了，只需在R环境中执行以下命令：

install.packages("RNAseqQC")

使用示例

首先，您需要准备一个counts矩阵和相应的表型数据：

library("RNAseqQC")
library("DESeq2")
library("ensembldb")
library("dplyr")
library("ggplot2")
library("purrr")
library("tidyr")
library("tibble")
library("magrittr")

count_mat <- counts(T47D)
meta <- data.frame(colData(T47D))
# 选取部分数据展示
count_mat[head(which(rowSums(count_mat) > 0)), 1:10]
#>                 veh_WT_rep1 veh_WT_rep2 veh_WT_rep3 veh_WT_rep4 veh_D538G_rep1
#> ENSG00000223972           5           3           3           6              4
#> ENSG00000278267           2           3           2           4              2
#> ENSG00000227232          80          91          95          97             65
#> ENSG00000243485           1           2           1           2              0
#> ENSG00000237613           0           0           0           1              0
#> ENSG00000239945           7           5           7           4              5
#>                 veh_D538G_rep2 veh_D538G_rep4 veh_Y537S_rep1 veh_Y537S_rep2
#> ENSG00000223972              1              6              3              5
#> ENSG00000278267              2              2              2              2
#> ENSG00000227232             70             73             74             65
#> ENSG00000243485              2              3              1              1
#> ENSG00000237613              0              0              0              0
#> ENSG00000239945              5              5              1              5
#>                 veh_Y537S_rep3
#> ENSG00000223972              3
#> ENSG00000278267              2
#> ENSG00000227232             89
#> ENSG00000243485              1
#> ENSG00000237613              0
#> ENSG00000239945              7
# 展示样本的元数据信息：
meta
#>                treatment mutation replicate
#> veh_WT_rep1          veh       WT      rep1
#> veh_WT_rep2          veh       WT      rep2
#> veh_WT_rep3          veh       WT      rep3
#> veh_WT_rep4          veh       WT      rep4
#> veh_D538G_rep1       veh    D538G      rep1
#> veh_D538G_rep2       veh    D538G      rep2
#> veh_D538G_rep4       veh    D538G      rep4
#> veh_Y537S_rep1       veh    Y537S      rep1
#> veh_Y537S_rep2       veh    Y537S      rep2
#> veh_Y537S_rep3       veh    Y537S      rep3
#> veh_Y537S_rep4       veh    Y537S      rep4
#> veh_D538G_rep3       veh    D538G      rep3
#> E2_WT_rep1            E2       WT      rep1
#> E2_WT_rep2            E2       WT      rep2
#> E2_WT_rep3            E2       WT      rep3
#> E2_WT_rep4            E2       WT      rep4
#> E2_Y537S_rep1         E2    Y537S      rep1
#> E2_Y537S_rep2         E2    Y537S      rep2
#> E2_Y537S_rep3         E2    Y537S      rep3
#> E2_Y537S_rep4         E2    Y537S      rep4
#> E2_D538G_rep1         E2    D538G      rep1
#> E2_D538G_rep2         E2    D538G      rep2
#> E2_D538G_rep3         E2    D538G      rep3
#> E2_D538G_rep4         E2    D538G      rep4

然后，利用这些数据创建一个DESeqDataSet对象，即可进行质量图表的绘制。

dds <- make_dds(counts = count_mat, metadata = meta, ah_record = "AH89426")

mcols(AnnotationHub::AnnotationHub()) %>%
  as_tibble(rownames = "record_id") %>%
  dplyr::filter(rdataclass == "EnsDb")

质量评估图表

• 总counts查看：```
  plot_total_counts(dds)

• 文库复杂性评估：```
  plot_library_complexity(dds)

• 基因检测与counts关系：```
  plot_gene_detection(dds)

• 基因类型分布：```
  plot_biotypes(dds)

• vst标准化及均值-标准差图：```
  vsd <- vst(dds)
  mean_sd_plot(vsd)

• 基因在染色体上的表达分布：```
  map(c("1", "5", "14"), ~plot_chromosome(vsd, .x))

• 生物学重复样本相关性：```
  colData(vsd)$trt_mut <- paste0(colData(vsd)$treatment, "_", colData(vsd)$mutation)
  ma_plots <- plot_sample_MAs(vsd, group = "trt_mut")
  cowplot::plot_grid(plotlist = ma_plots[17:24], ncol = 2)

• 热图聚类：```
  set.seed(1)
  plot_sample_clustering(vsd, anno_vars = c("treatment", "mutation", "replicate"))

• PCA分析：```
  plot_pca(vsd, PC_x = 1, PC_y = 2, color_by = "treatment")

• 差异表达分析：```
  dds <- estimateSizeFactors(dds)
  plot_gene("CLEC2B", dds, x_var = "mutation", color_by = "treatment")

• 多因子差异分析：```
  dds$mutation <- as.factor(dds$mutation)
  dds$treatment <- as.factor(dds$treatment)
  design(dds) <- ~ mutation + treatment
  dds <- DESeq(dds, parallel = T)
  plotDispEsts(dds)

• MA  plot:

# get testing results
de_res <- lfcShrink(dds, coef = "mutation_WT_vs_D538G", lfcThreshold = log2(1.5), type = "normal", parallel = TRUE)

# MA plot
plot_ma(de_res, dds, annotate_top_n = 5)

plot_ma(de_res, dds, highlight_genes = c("CLEC2B", "PAGE5", "GAPDH"))

参考文档：

https://cran.r-project.org/web/packages/RNAseqQC/vignettes/introduction.html

作者：生物信息与育种，请关注同名微信公众号：生物信息与育种。