2024-06-17

tidybulk—摈弃繁琐代码，在巨人肩膀学习转录组（二）

原创小果生信果 2023-12-14 19:00:11

收录于话题

今天小果来和一起学习tidybulk包的第二节。不知道大家有没有看完小果第一节的讲解呢？如果还没有的话，可以进入公众号一学习一下前面的内容哦~ 小果在第一节中为大家讲解了转录组分析的行名整合，数据的标准化，筛选差异基因以及降维部分的代码。这些在转录组学的分析中都是十分重要的内容哦~因此小果在第二章中也会省略一部分重复的内容哦，当然小果也会进行相应的注释，实在有不懂的小伙伴翻看一下第一节呢。好啦~那现在就让小果陪着大家一起开始今天第二节的学习吧！

加载R包

如果是两节一起学习的小伙伴这一步就可以省略啦~

library(tidyverse)library(tidybulk)library(SummarizedExperiment)

差异分析

正如小果上节所讲的，我们做转录组分析想要研究的最主要的一个方向就是不同样本之间的差异，这一节就是在研究这一个方向。在我们现有的研究中，我们可以使用edgeR或是limma等软件进行计算。通过test_differential_abundance函数，我们就可以得出一个包含对数倍数变化，p-value的tibble文件。接下随着小果的步伐一起来看看如何使用tidybulk包进行差异分析吧！

# 这里相当于其中的一个样本和其余的所有样本进行对比            counts_SE.de =              counts_SE %>%              test_differential_abundance( ~ condition, action="get")    counts_SE.de %>% tidySummarizedExperiment::as_tibble() -> res1            head(res1)## # A tibble: 6 × 13            ##   .feature .sample    count Cell.type time  condition batch factor_of_interest            ##                                                                                                           ## 1 A1BG     SRR1740034   153 b_cell    0 d   TRUE      0     TRUE                          ## 2 A1BG-AS1 SRR1740034    83 b_cell    0 d   TRUE      0     TRUE                          ## 3 AAAS     SRR1740034   868 b_cell    0 d   TRUE      0     TRUE                          ## 4 AACS     SRR1740034   222 b_cell    0 d   TRUE      0     TRUE                          ## 5 AAGAB    SRR1740034   590 b_cell    0 d   TRUE      0     TRUE                          ## 6 AAMDC    SRR1740034    48 b_cell    0 d   TRUE      0     TRUE                          # 添加.contrasts参数过后，小果在这指定了两个样本进行差异分析。counts_SE.de.contr =              counts_SE %>%              identify_abundant(factor_of_interest = condition) %>%              test_differential_abundance(~0+ condition, .contrasts =c( "conditionFALSE - conditionTRUE"),                action="get"              )counts_SE.de.contr %>%              tidySummarizedExperiment::as_tibble() -> res2            

            head(res2)## # A tibble: 6 × 14            ##   .feature .sample    count Cell.type time  condition batch factor_of_interest            ##                                                                                                           ## 1 A1BG     SRR1740034   153 b_cell    0 d   TRUE      0     TRUE                          ## 2 A1BG-AS1 SRR1740034    83 b_cell    0 d   TRUE      0     TRUE                          ## 3 AAAS     SRR1740034   868 b_cell    0 d   TRUE      0     TRUE                          ## 4 AACS     SRR1740034   222 b_cell    0 d   TRUE      0     TRUE                          ## 5 AAGAB    SRR1740034   590 b_cell    0 d   TRUE      0     TRUE                          ## 6 AAMDC    SRR1740034    48 b_cell    0 d   TRUE      0     TRUE

去批次

当进行转录组的数据整合的时候，测序的批次会对数据可能会有一定的影响，这时候各位小伙伴就需要进行去批次的处理！

# 数据的准备，这里在第一节中已经运行过啦，不过为了分开两节学习的小伙伴，小果在这就重复一下代码，各位小伙伴直接跑就可以啦            rowData(counts_SE)$gene_name =rownames(counts_SE)            

            counts_SE.aggr = counts_SE %>%aggregate_duplicates(.transcript = gene_name)## tidybulk says: your object does not have duplicates along the gene_name column. The input dataset is returned.counts_SE.norm = counts_SE.aggr %>%identify_abundant(factor_of_interest = condition) %>%scale_abundance()## tidybulk says: the sample with largest library size SRR1740080 was chosen as reference for scaling# 在这一步中可能需要下载R包，那因为环境不一样，各位小伙伴只需要根据提示进行R包的安装就可以啦，如果遇到不好解决的问题，也可以在评论留言，小果看到后也会为大家的问题进行一一解答哦!            counts_SE.norm.adj =                counts_SE.norm %>%adjust_abundance(    ~ factor_of_interest + batch)## Found2batches## Adjusting for1covariate(s) or covariate level(s)## Standardizing Data across genes## Fitting L/S model and finding priors## Finding parametric adjustments## Adjusting the Data

卷积分辨细胞类型（机器学习）

在转录组学的研究中，细胞类型的确定是很关键的一项，在这个R包中，为各位小伙伴提供了一个大家可能平时不太会接触到的确定细胞类型的方法，那就是卷积，但如果要讲解卷积具体是什么东西就会涉及到很多内容啦，小果今天就不详细介绍啦，有兴趣的小伙伴可以通过查看公众号的其他推文进行学习哦~ 那么下面就让小果可以实操吧！

counts_SE.cibersort =                counts_SE %>%                deconvolve_cellularity(action="get", cores=1, prefix ="cibersort__")             

            counts_SE.cibersort %>%                pivot_longer(                    names_to="Cell_type_inferred",                     values_to ="proportion",                     names_prefix ="cibersort__",                     cols=contains("cibersort__")                ) %>%              ggplot(aes(x=`Cell_type_inferred`, y=proportion, fill=`Cell.type`)) +              geom_boxplot() +              facet_wrap(~`Cell.type`) +              theme_bw() +              theme(axis.text.x =element_text(angle =90, hjust =1, vjust =0.5), aspect.ratio=1/5)## tidySummarizedExperiment says: A data frame is returned for independent data analysis.

# 差异细胞类型的丰度计算 tidybulk包可以对不同条件下的差异细胞类型丰度进行统计测试。

counts_SE %>%              test_differential_cellularity(. ~ condition ) %>%head()## # A tibble: 6 × 7            ##   .cell_type cell_type_proportions `estimate_(Intercept)` estimate_conditionTRUE            ##                                                                                                         ## 1 cibersort…                      -2.94                  2.25           ## 2 cibersort…                      -4.86                  1.48           ## 3 cibersort…                      -5.33                 -0.487           ## 4 cibersort…                      -2.33                  0.924           ## 5 cibersort…                      -2.83                 -0.620           ## 6 cibersort…                      -2.46                  0.190           ## # ℹ 3 more variables: std.error_conditionTRUE,           ## #   statistic_conditionTRUE, p.value_conditionTRUE

紧跟着，还可以使用得到的数据进行回归分析。

counts_SE_survival =                counts_SE %>%                nest(data =-sample) %>%                    mutate(                        days =sample(1:1000, size =n()),                        dead =sample(c(0,1), size =n(), replace =TRUE)                    ) %>%                unnest(data)     counts_SE_survival %>%                test_differential_cellularity(survival::Surv(days, dead) ~ .) %>%head()## # A tibble: 6 × 6            ##   .cell_type          cell_type_proportions estimate std.error statistic p.value            ##                                                                                                                   ## 1 cibersort.B.cells.…       0.379      0.685     0.554 0.580                   ## 2 cibersort.B.cells.…       0.201      0.346     0.581 0.561                   ## 3 cibersort.Plasma.c…      -3.34       1.02     -3.28  0.00102           ## 4 cibersort.T.cells.…      -2.58       1.69     -1.53  0.127                   ## 5 cibersort.T.cells.…       0.0698     0.513     0.136 0.892                   ## 6 cibersort.T.cells.…       0.383      0.493     0.778 0.437

接着小果会带着大家进行Kaplan-Meier 曲线的检验并进行可视化。

counts_stratified =                counts_SE_survival %>%            

                # Test                test_stratification_cellularity(                    survival::Surv(days, dead) ~ .,                    sample, transcript, count                )            

            head(counts_stratified)## # A tibble: 6 × 6            ##   .cell_type cell_type_proportions .low_cellularity_exp…¹ .high_cellularity_ex…²            ##                                                                                                         ## 1 cibersort…                       11.3                  11.7           ## 2 cibersort…                       17.7                   5.32           ## 3 cibersort…                       13.3                   9.72           ## 4 cibersort…                       13.8                   9.25           ## 5 cibersort…                       10.7                  12.3           ## 6 cibersort…                       13.6                   9.45           ## # ℹ abbreviated names: ¹.low_cellularity_expected, ².high_cellularity_expected            ## # ℹ 2 more variables: pvalue, plotcounts_stratified$plot[[1]]

聚类

在各位小伙伴的转录组学研究中，相信大家一定对聚类问题不陌生，成功的聚类，可以让研究人员看出个细胞间的异同点，tidybulk包中为各位准备了两种聚类的方法，分别是K-means和SNN（但SNN目前R包因Matrix的更新发生报错，暂不可使用），下面就让小果来为大家分别实操一下。

K-means

# 这里小果会用到上一节中的降维结果              counts_SE.norm.MDS =                counts_SE.norm %>%                reduce_dimensions(method="MDS", .dims =6)## Getting the 500 most variable genes## tidybulk says: to access the raw results do `attr(..., "internals")$MDS`counts_SE.norm.cluster = counts_SE.norm.MDS %>%                cluster_elements(method="kmeans", centers =2, action="get" )              

              # 可以将聚类注释添加到 MDS 降维数据集和绘图中。               counts_SE.norm.cluster %>%                  ggplot(aes(x=`Dim1`, y=`Dim2`, color=`cluster_kmeans`)) +                geom_point() +                theme_bw()

去除冗余数据

在研究的过程中各位小伙伴可能会希望从原始数据集中删除冗余的一些数据，作者也想到了这点，为大家准备了对应的函数。有两种方法：方法一：通过使用method=“correlation”参数删除高度相关的簇（保留具有代表性的簇）方法二：在降维空间中移除最近的元素对。## 方法1

counts_SE.norm.non_redundant =                counts_SE.norm.MDS %>%              remove_redundancy(    method ="correlation" )## Getting the 8513 most variable genescounts_SE.norm.non_redundant %>%                pivot_sample() %>%                ggplot(aes(x=`Dim1`, y=`Dim2`, color=`Cell.type`)) +              geom_point() +              theme_bw()

方法2

counts_SE.norm.non_redundant =                  counts_SE.norm.MDS %>%                remove_redundancy(                  method ="reduced_dimensions",                  Dim_a_column =`Dim1`,                  Dim_b_column =`Dim2`                )              

              counts_SE.norm.non_redundant %>%                  pivot_sample() %>%                  ggplot(aes(x=`Dim1`, y=`Dim2`, color=`Cell.type`)) +                geom_point() +                theme_bw()

其他包装函数

从BAM/SAM到基因计数的转换

counts =tidybulk_SAM_BAM(                  file_names,                  genome ="hg38",                  isPairedEnd =TRUE,                  requireBothEndsMapped =TRUE,                  checkFragLength =FALSE,                  useMetaFeatures =TRUE              )

ensembleID和symbolID之间的转换

counts_ensembl %>%ensembl_to_symbol(ens) %>%head()## # A tibble: 6 × 8              ##   ens    iso   `read count` sample cases_0_project_dise…¹ cases_0_samples_0_sa…²              ##                                                                                                                               ## 1 ENSG0… 13             144 TARGE… Acute Myeloid Leukemia Primary Blood Derived…              ## 2 ENSG0… 13              72 TARGE… Acute Myeloid Leukemia Primary Blood Derived…              ## 3 ENSG0… 13               0 TARGE… Acute Myeloid Leukemia Primary Blood Derived…              ## 4 ENSG0… 13            1099 TARGE… Acute Myeloid Leukemia Primary Blood Derived…              ## 5 ENSG0… 13              11 TARGE… Acute Myeloid Leukemia Primary Blood Derived…              ## 6 ENSG0… 13               2 TARGE… Acute Myeloid Leukemia Primary Blood Derived…              ## # ℹ abbreviated names: ¹cases_0_project_disease_type,              ## #   ²cases_0_samples_0_sample_type              ## # ℹ 2 more variables: transcript, ref_genome

symbol和基因全名之间的转换

counts_SE %>%                  describe_transcript() %>%                  select(feature, description, everything()) %>%print(n =5)## # A SummarizedExperiment-tibble abstraction: 408,624 × 48              ## # [90mFeatures=8513 | Samples=48 | Assays=count[0m              ##    feature  sample     count Cell.type time  condition batch factor_of_interest              ##                                                                                                                    ##  1 A1BG     SRR1740034   153 b_cell    0 d   TRUE      0     TRUE                            ##  2 A1BG-AS1 SRR1740034    83 b_cell    0 d   TRUE      0     TRUE                            ##  3 AAAS     SRR1740034   868 b_cell    0 d   TRUE      0     TRUE                            ##  4 AACS     SRR1740034   222 b_cell    0 d   TRUE      0     TRUE                            ##  5 AAGAB    SRR1740034   590 b_cell    0 d   TRUE      0     TRUE                            ##  6 AAMDC    SRR1740034    48 b_cell    0 d   TRUE      0     TRUE                            ##  7 AAMP     SRR1740034  1257 b_cell    0 d   TRUE      0     TRUE                            ##  8 AANAT    SRR1740034   284 b_cell    0 d   TRUE      0     TRUE                            ##  9 AAR2     SRR1740034   379 b_cell    0 d   TRUE      0     TRUE                            ## 10 AARS2    SRR1740034   685 b_cell    0 d   TRUE      0     TRUE                            ## # ℹ 40 more rows              ## # ℹ 2 more variables: description, gene_name