内容

1简介

选择性聚腺苷酸化(Alternative polyadenylation, APA)是真核生物中最重要的转录后调控机制之一。像选择性剪接一样，APA可以增加转录组多样性。此外，它还定义了3 ' UTR，并导致基因表达的改变。这是一个严格控制的过程，APA的错误调节可以影响许多生物过程，如不受控制的细胞周期和生长。虽然已经开发了几种高通量测序方法，但用于识别APA事件的数据仍然有限。

然而，最初为量化全基因组基因表达而创建的大规模RNA-seq数据集可在GEO和TCGA等公共数据库中获得。这些RNA-seq数据集还包含全基因组APA的信息。因此，我们开发了InPAS包，用于从传统的RNA-seq数据中识别APA。

InPAS工作流的主要流程如下:

此外，InPAS包还提供了对RNA-seq数据覆盖进行质量控制的功能，可视化感兴趣基因的近端和远端CP位点的差异使用，并为基因集富集分析(GSEA)准备必要的文件，以揭示具有替代CP位点的基因的生物学见解。

2如何运行InPAS

首先，加载所需软件包，包括InPAS，以及种特异性基因组和基因组注释数据库:BSgenome、TxDb和EnsDb。

logger <- file(tempfile()， open =" wt") sink(logger, type="message") suppressPackageStartupMessages({library(InPAS) library(BSgenome.Mmusculus.UCSC.mm10) library(BSgenome.Hsapiens.UCSC.hg19) library(TxDb.Mmusculus.UCSC.mm10.knownGene) library(EnsDb.Hsapiens.v86) library(EnsDb.Mmusculus.v79) library(cleanUpdTSeq) library(RSQLite) library(future.apply)})

2.1第一步:建立SQLite数据库

在数据库中创建了7个表。表“元数据”存储元数据，包括标签(样本名称)、条件(实验处理组)、bedgraph_file (BEDGraph文件路径)和深度(全基因组覆盖深度)，这些信息最初设置为零，随后在分析过程中更新。表“sample_coverage”，“chromosome_coverage”，“total_coverage”，“utr3_total_coverage”，“CPsites”和“utr3cds_coverage”存储中间文件的名称以及与文件相关的染色体和标记(示例名称)。

计划(顺序)data_dir <- system。file("extdata"， package = "InPAS") bedgraphs <- c(file. txt)path(data_dir， "Baf3.extract.bedgraph")，文件。path(data_dir， "UM15.extract.bedgraph")) hugeData <- FALSE genome <- BSgenome.Mmusculus.UCSC. path(data_dir， "UM15.extract.bedgraph")mm10 tags <- c("Baf3"， "UM15") metadata <- data.frame(tag = tags, condition = c("Baf3"， "UM15")， bedgraph_file = bedgraphs) ##在现实中，不要使用临时目录进行分析。相反，使用##持久目录保存分析输出。Outdir = tempdir()写入。表(元数据，文件=文件。path(outdir =outdir， "metadata.txt")， sep = "\t"， quote = FALSE, row.names = FALSE) sqlite_db <- setup_sqlitedb(metadata = file. txt)##查看数据库db_conn <- dbConnect(drv = RSQLite::SQLite()， dbname = sqlite_db) dbListTables(db_conn)
##[1]“CPsites”“chromosome_coverage”“genome_anno”##[4]“metadata”“sample_coverage”“total_coverage”##[7]“utr3_total_coverage”“utr3cds_coverage”
dbReadTable (db_conn,“元数据”)
##标签条件## 1 Baf3 Baf3 ## 2 UM15 UM15 ## bedgraph_file depth ## 1 /tmp/RtmpQEX3nM/Rinst1c93e55c584678/InPAS/extdata/Baf3.extract。床图0 ## 2 /tmp/RtmpQEX3nM/Rinst1c93e55c584678/InPAS/extdata/UM15.extract. txt . txt . txt . txt . txt。bedgraph 0
dbDisconnect (db_conn)

2．2步骤2:提取3 ' UTR注释

使用该函数提取3 ' UTR注释，包括起始和结束坐标，3 ' UTR的链信息，最后一个CDS以及3 ' UTR的3 '末端与直接下游外显子之间的间隙extract_UTR3Anno来自基因组注释数据库:一个TxDb数据库和一个ensemble bldb数据库用于感兴趣的物种。为了演示，下面的R脚本片段展示了如何从人类参考基因组(hg19)的简化TxDb中提取3 ' UTR注释。实际上，用户应该使用TxDb对用于RNA-seq读取比对的PRIMARY参考基因组组装(不包括替代补丁)进行最可靠的基因组注释。如果对感兴趣的物种没有TxDb，用户可以使用基因组特征包中的makeTxDbFromUCSC、makeTxDbFromBiomart、makeTxDbFromEnsembl或makeTxDbFromGFF函数构建TxDb，具体取决于基因组注释文件的来源。

Samplefile <- system。文件(“extdata”、“hg19_knownGene_sample。sqlite"， package="GenomicFeatures") TxDb <- loadDb(samplefile) seqnames <- seqnames(BSgenome.Hsapiens.UCSC.hg19) #排除线粒体基因组和替代单倍型chr2exclude <- c("chrM"， "chrMT"， seqnames[grepl("_(hap\\d+|fix|alt)$"， seqnames, perl = TRUE)]) #设置InPAS分析的全局变量set_globals(genome = BSgenome.Hsapiens.UCSC。hg19, TxDb = TxDb, EnsDb = EnsDb. hsapiens。V86, outdir = tempdir()， chr2exclude = chr2exclude, lockfile = tempfile())
##设置默认基因组为hg19。
##设置默认EnsDb为EnsDb。
##设置默认TxDb为TxDb。
utr3_anno <- extract_UTR3Anno(sqlite_db = sqlite_db, TxDb = getInPASTxDb()， edb = getInPASEnsDb()， genome = getInPASGenome()， outdir = getInPASOutputDirectory()， chr2exclude = getChr2Exclude()))
##警告:无法映射42个请求id中的3个。
头(utr3_anno chr1美元)
##具有6个范围和8个元数据列的GRanges对象:## seqnames ranges ##   ## chr1: lastut3:uc001bum2_5|IQCC|utr3 chr1 32673684-32674288 ## chr1:lastutr3:uc001fbq。3._3|S100A9|utr3 chr1 153333315-153333503 ## chr1:lastutr3:uc031pqa.1_13|100129405|utr3 chr1 155719929-155720673 ## chr1:lastutr3:uc001gde.2_2|LRRC52|utr3 chr1 165533062-165533185 ## chr1:lastutr3:uc001hfg.3_15|PFKFB2|utr3 chr1 207245717-207251162 ## chr1:lastutr3:uc001hfh.3_15|PFKFB2|utr3 chr1 207252365-207254368 ## strand | exon_rank transcript ##  |   ## chr1:lastutr3:uc001bum.2_5|IQCC|utr3 + | 5 uc001bum.2 ## chr1:lastutr3:uc001fbq.3_3|S100A9|utr3 + | 3 uc001fbq.3 ## chr1:lastutr3:uc031pqa.1_13|100129405|utr3 + | 13 uc031pqa.1 ## chr1:lastutr3:uc001gde.2_2|LRRC52|utr3 + | 2 uc001gde.2 ## chr1:lastutr3:uc001hfg.3_15|PFKFB2|utr3 + | 15 uc001hfg.3 ## chr1:lastutr3:uc001hfh.3_15|PFKFB2|utr3 + | 15 uc001hfh.3 ## feature gene ##   ## chr1:lastutr3:uc001bum.2_5|IQCC|utr3 utr3 55721 ## chr1:lastutr3:uc001fbq.3_3|S100A9|utr3 utr3 6280 ## chr1:lastutr3:uc031pqa.1_13|100129405|utr3 utr3 100129405 ## chr1:lastutr3:uc001gde.2_2|LRRC52|utr3 utr3 440699 ## chr1:lastutr3:uc001hfg.3_15|PFKFB2|utr3 utr3 5208 ## chr1:lastutr3:uc001hfh.3_15|PFKFB2|utr3 utr3 5208 ## exon symbol ##   ## chr1:lastutr3:uc001bum.2_5|IQCC|utr3 chr1:lastutr3:uc001b.. IQCC ## chr1:lastutr3:uc001fbq.3_3|S100A9|utr3 chr1:lastutr3:uc001f.. S100A9 ## chr1:lastutr3:uc031pqa.1_13|100129405|utr3 chr1:lastutr3:uc031p.. 100129405 ## chr1:lastutr3:uc001gde.2_2|LRRC52|utr3 chr1:lastutr3:uc001g.. LRRC52 ## chr1:lastutr3:uc001hfg.3_15|PFKFB2|utr3 chr1:lastutr3:uc001h.. PFKFB2 ## chr1:lastutr3:uc001hfh.3_15|PFKFB2|utr3 chr1:lastutr3:uc001h.. PFKFB2 ## annotatedProximalCP truncated ##   ## chr1:lastutr3:uc001bum.2_5|IQCC|utr3 unknown FALSE ## chr1:lastutr3:uc001fbq.3_3|S100A9|utr3 unknown FALSE ## chr1:lastutr3:uc031pqa.1_13|100129405|utr3 proximalCP_155720479 FALSE ## chr1:lastutr3:uc001gde.2_2|LRRC52|utr3 unknown FALSE ## chr1:lastutr3:uc001hfg.3_15|PFKFB2|utr3 unknown FALSE ## chr1:lastutr3:uc001hfh.3_15|PFKFB2|utr3 unknown FALSE ## ------- ## seqinfo: 83 sequences from an unspecified genome

本插图将使用为小鼠参考基因组mm10准备的3 ' UTR注释进行后续演示

##设置小鼠InPAS分析的全局变量set_globals(genome = BSgenome.Mmusculus.UCSC.基因组= BSgenome.Mmusculus.UCSC.基因组。mm10, TxDb = TxDb. mmusculus . ucsc .mm10。knownGene, EnsDb = EnsDb. mmusculus。v79, outdir = tempdir()， chr2exclude = "chrM"， lockfile = tempfile()))
##设置默认基因组为mm10。
##设置默认EnsDb为EnsDb。
##设置默认TxDb为TxDb。
tx <- parse_TxDb(sqlite_db = sqlite_db, TxDb = getInPASTxDb()， edb = getInPASEnsDb()， genome = getInPASGenome()， outdir = getInPASOutputDirectory()， chr2exclude = getChr2Exclude()))
##警告:无法映射24594个请求id中的6032个
#加载R对象:ut3。mm10 data(utr3.mm10) ##将3' UTR注释utr3的grange转换为GRangesList。Mm10 <- split(utr3。mm10, seqnames(utr3.mm10)， drop = TRUE)

2．3步骤3:重新格式化覆盖率数据

在此步骤之前，应从BAM文件中准备以BEDGraph格式的基因组覆盖，这些BAM文件由RNA-seq数据对齐使用genomecov命令在BEDTools套件。可以过滤BAM文件以删除多映射对齐、映射质量低的对齐等。参考命令如下:

对于与star# # -q 255对齐的单端RNA-seq数据，唯一映射samtools view -bu -h -q 255 /path/to/XXX.SE。-bu -h -q 255 /path/to/XXX.PE. bga -split > XXX.SE.uniq.bedgraph ##| \ bedtools genome ecov -ibam - -bga -split > XXX.PE.uniq.bedgraph

将BEDGraph格式中的基因组覆盖数据转换为rle类的R对象get_ssRleCov每个样本的每个染色体的功能。每个染色体的对象保存到outdir．文件名、标签(样本名)和染色体名保存到表" sample_coverage "中。随后，所有样本的染色体特异性Rle对象被组装成一个两级Rle对象列表，第1级是染色体名称，第2级是每个标记(样本名称)的Rle。值得注意的是，这里使用的样本BEDGraph文件仅包含小鼠参考基因组mm10的“chr6”的覆盖数据。

coverage <- list() for (i in seq_along(bedgraphs)){coverage[[tags[i]]]] <- get_ssRleCov(bedgraph = bedgraphs[i]， tag = tags[i]， genome = genome, sqlite_db = sqlite_db, outdir = outdir, chr2exclude = getChr2Exclude())} coverage <- assemble_allCov(sqlite_db, seqname = "chr6"， outdir, genome = getInPASGenome()))

此时，用户可以根据所有表达基因和3 ' utr的覆盖率来检查数据质量run_coverageQC．这个函数输出概括的覆盖率指标:率,expressed.gene.coverage。率,UTR3.coverage。和utr3 . expression .gene.sub - set.coverage.rate。对于表达基因的3 ' utr，质量数据的覆盖率应大于0.75。

如果(.Platform $ OS。type == "windows"){计划(多会话)}else{计划(多核)}run_coverageQC(sqlite_db, TxDb = getInPASTxDb()， edb = getInPASEnsDb()， genome = getInPASGenome()， chr2exclude = getChr2Exclude()， which = GRanges("chr6"， ranges = IRanges(98013000, 140678000))))
##链信息将被忽略。
UTR3.coverage. raterate ## Baf3 0.003463505 0.5778441 0.01419771 ## UM15 0.003428528 0.5719748 0.01405159 ## utr3 . expression .gene.sub - set.coverage。利率## Baf3 0.8035821 ## UM15 0.7953112
计划(顺序)

2.4步骤4:识别潜在CP位点

深度权重，Z-score截止阈值，以及在每个条件(大数据)或单个样本(非大数据)中生物重复合并的3 ' UTRs的总覆盖率setup_CPsSearch函数。潜在的新CP位点被识别为每条染色体使用search_CPs函数。通过提供的朴素贝叶斯分类器可以过滤和/或调整这些潜在的CP位点cleanUpdTseq和/或通过简单地匹配六聚体聚腺苷酸信号的位置权重矩阵(PWM)来使用聚腺苷酸评分(AAUAAA等)。

##从cleanUpdTseq包数据(分类器)prepared_data <- setup_CPsSearch(sqlite_db, genome = getInPASGenome()， chr.;Utr3 = Utr3。mm10$chr6, seqname = "chr6"， background = "10K"， TxDb = getInPASTxDb()， hugeData = TRUE, outdir = outdir, silence = TRUE, coverage_threshold = 5) CPsites <- search_CPs(seqname = "chr6"， sqlite_db = sqlite_db, genome = getInPASGenome()， MINSIZE = 10, window_size = 100, search_point_START = 50, search_point_END = NA, cutEnd = 0, filter. filter.)last = TRUE, adjust_distal_polyA_end = TRUE, long_coverage_threshold = 2, PolyA_PWM = NA, classifier = classifier, classifier_cutoff = 0.8, shift_range = 50, step = 5, outdir = outdir, silence = TRUE)
没有找到可读的配置文件
##创建注册表在'/tmp/RtmpMeoG6g/006.CPsitesout_chr6'使用集群函数'交互式'
##增加1个工作…
##使用集群函数“交互式”在1个块中提交1个作业…

2.5步骤5:估计近端和远端CP位点的使用情况

根据短3 ' utr和长3 ' utr的读覆盖估计近端和远端CP位点的使用情况

utr3_cds_cov <- get_regionCov(chr. cov)Utr3 = Utr3。mm10[["chr6"]]， sqlite_db, outdir, phmm = FALSE) eSet <- get_UTR3eSet(sqlite_db, normalize ="none"， singleSample = FALSE)

2.6步骤6。识别PDUI差异事件

InPAS提供该功能test_dPDUI根据实验设计，利用不同的统计模型，确定不同情况下近端CP位点和远端CP位点的使用差异。InPAS提供了单样本差分PDUI分析和单组分析的统计方法。此外，InPAS还为两组不可复制设计提供了Fisher精确检验，并为更复杂的设计提供了利用limma包的经验贝叶斯线性模型。属性可以进一步筛选测试结果filter_testOut基于每个条件下的分数样本的函数，其中包含已识别的差分PDUI事件的覆盖数据，和/或标称p值、调整p值或log2(折叠变化)的截止值。

test_out <- test_dPDUI(eset = eset, method = "fisher. "Exact "， normalize = "none"， sqlite_db = sqlite_db)
filter_out <- filter_testOut(res = test_out, gp1 = "Baf3"， gp2 = "UM15"， background_coverage_threshold = 2, P.Value_cutoff = 0.05, adj.P。Val_cutoff = 0.05, dPDUI_cutoff = 0.3, PDUI_logFC_cutoff = 0.59)

2.7步骤7。可视化dPDUI事件并为GSEA准备文件

InPAS包还提供功能，get_usage4plot，plot_utr3Usage,setup_GSEA，以可视化感兴趣基因的近端和远端CP位点的差异使用，并为基因集富集分析(GSEA)准备必要的文件，以揭示具有替代CP位点的基因的生物学见解。

##直观dPDUI事件gr <- GRanges("chr6"， IRanges(128846245, 128850081)， strand = "-") names(gr) <- "128846245-128850081" data4plot <- get_usage4plot(gr, proximalSites = 128849130, sqlite_db, hugeData = TRUE) plot_utr3Usage(usage_data = data4plot, vline_color = "紫色"，vline_type = "虚线")

##为GSEA准备一个排名文件setup_GSEA(eset = test_out, groupList= list(Baf3 =" Baf3"， UM15 ="UM15")， outdir = outdir, preranked = TRUE, rankBy =" logFC"， rnkFilename =" InPAS.rnk")

3.会话信息

sessionInfo ()

R版本4.2.1(2022-06-23)平台:x86_64-pc-linux-gnu(64位)运行在Ubuntu 20.04.5 LTS下

矩阵产品:默认BLAS: /home/biocbuild/bbs-3.16-bioc/R/lib/libRblas。so LAPACK: /home/biocbuild/bbs-3.16-bioc/R/lib/libRlapack.so

locale: [1] LC_CTYPE=en_US。utf - 8 LC_NUMERIC = C
[3] LC_TIME=en_GB LC_COLLATE=C
[5] LC_MONETARY = en_US。utf - 8LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER = en_US。utf - 8 LC_NAME = C
[9] lc_address = c lc_phone = c
[11] LC_MEASUREMENT = en_US。utf - 8 LC_IDENTIFICATION = C

附加的基本包:[1]stats4 stats graphics grDevices utils数据集方法
[8]的基础

其他附加包:[1]future.apply_1.9.1
[2] future_1.28.0
[3] RSQLite_2.2.18
[4] cleanUpdTSeq_1.36.0
[5] BSgenome.Drerio.UCSC.danRer7_1.4.0
[6] EnsDb.Mmusculus.v79_2.99.0
[7] EnsDb.Hsapiens.v86_2.99.0
[8] ensembldb_2.22.0
[9] AnnotationFilter_1.22.0
[11]基因组特征_1.50.0
[12] AnnotationDbi_1.60.0
[13] Biobase_2.58.0
[14] BSgenome.Hsapiens.UCSC.hg19_1.4.3
[15] BSgenome.Mmusculus.UCSC.mm10_1.4.3
[16] BSgenome_1.66.0
[17] rtracklayer_1.58.0
[18] Biostrings_2.66.0
[19] XVector_0.38.0
[20] GenomicRanges_1.50.0
[21] GenomeInfoDb_1.34.0
[22] IRanges_2.32.0
[23] S4Vectors_0.36.0
[24] BiocGenerics_0.44.0
[25] InPAS_2.6.0
[26] BiocStyle_2.26.0

通过命名空间加载(且未附加):[1]backports_1.4.1 BiocFileCache_2.6.0
[3] plyr_1.8.7 lazyeval_0.2.2
[5] BiocParallel_1.32.0 listenv_0.8.0
[7] ggplot2_3.3.6 digest_0.6.30
[9] htmltools_0.5.3 magick_2.7.3
[11] fansi_1.0.3 magrittr_2.0.3
[13] Rsolnp_1.16 checkmate_2.1.0
[15] memoise_2.0.1 base64url_1.4
[17] tzdb_0.3.0 limma_3.54.0
[19] globals_0.16.1 readr_2.1.3
[21] matrixStats_0.62.0 vroom_1.6.0
[23] prettyunits_1.1.1 colorspace_2.0-3
[25] blob_1.2.3 rappdirs_0.3.3
[27] xfun_0.34 dplyr_1.0.10
[29] crayon_1.5.2 RCurl_1.98-1.9
[31] jsonlite_1.8.3 brew_1.0-8
[33] flock_0.7 glue_1.6.2
[35] gtable_0.3.1 zlibbioc_1.44.0
[37] seqinr_2 .2-16 DelayedArray_0.24.0
[39] plyranges_1.18.0 depmixS4_1.5-0
[41] scales_1.2.1 DBI_1.1.3 .
[43] Rcpp_1.0.9 progress_1.2.2 .
[45] archive_1.1.5 bit_4.0.4
[47] proxy_0.4-27 preprocessCore_1.60.0
[49] truncnorm_1.0-8 httr_1.4.4
[51] ellipsis_0.3.2 farver_2.1.1
XML_3.99-0.12 . [53] pkgconfig_2.0.3
[55] nnet_7.3-18 sass_0.4.2 .
[57] dbplyr_2.2.1 utf8_1.2.2
[59] labeling_0.4.2 tidyselect_1.2.0
[61] rlang_1.0.6 reshape2_1.4.4
[63] [au:] munsell_0.5.0 tools_4.2.1 .
[65] cachem_1.0.6 cli_3.4.1
[67] generics_0.1.3 ade4_1.7-19
[69] evaluate_0.17 string_1 .4.1
[71] fastmap_1.1.0 yaml_2.3.6
[73] fs_1.5.2 knitr_1.40
[75] [au:] bit64_4.0.5 KEGGREST_1.38.0 .
[77] nlme_1 .1-160 xml2_1.3.3
[79] bioart_2 .54.0 debugme_1.1.0
[81] compiler_4.2.1 filelock_1.0.2
[83] curl_4.3.3 png_0.1-7
[85] e1071_1.7-12 tibble_3.1.8
[87] [font =宋体]
[89] [endnoteref: 1
[91] [font =宋体]
[93] vctrs_0.5.0 pillar_1.8.1
[95]生命周期_1.0.3 BiocManager_1.30.19
[97] jquerylib_0.1.4 data.table_1.14.4
[99] [cn
[101] BiocIO_1.8.0 bookdown_0.29
[103] parallely_1.32.1 codetools_0.2-18
[105] MASS_7.3-58.1 assertthat_0.2.1
[107] j . jjson_0.2.21
[109]与thr_2.5.0基因组比对s_1.34.0
[111] batchtools_0.9.15 Rsamtools_2.14.0
[113] GenomeInfoDbData_1.2.9 parallel_4.2.1
[115] hms_1.1.2 grid_4.2.1
[117]王晓明
[119]刘志军

水槽(type = "消息")关闭(日志)

1.Sheppard S.， Lawson N. D.， Zhu L. J.利用朴素贝叶斯分类器精确识别3 '端深测序中的多聚腺苷酸位点。生物信息学29日,2564 - 2571(2013)。