R包：APAlyzer从RNA-seq数据计算APA表达丰度

介绍

今天安利APAlyzer工具，它是通过RNA-seq数据获取3′UTR APA, intronic APA等表达谱的R包。

APAlyzer将bam文件比对到PolyA-DB数据库识别APA。

Most eukaryotic genes produce alternative polyadenylation (APA) isoforms. APA is dynamically regulated under different growth and differentiation conditions. Here, we present a bioinformatics package, named APAlyzer, for examining 3′UTR APA, intronic APA and gene expression changes using RNA-seq data and annotated polyadenylation sites in the PolyA_DB database. Using APAlyzer and data from the GTEx database, we present APA profiles across human tissues.

教程

library(APAlyzer)
library(TBX20BamSubset)
library(Rsamtools)# RNA-seq BAM files
flsall = getBamFileList()# Genomic reference
library(repmis)
URL="https://github.com/RJWANGbioinfo/PAS_reference_RData/blob/master/"
file="mm9_REF.RData"
source_data(paste0(URL,file,"?raw=True"))# Building 3’UTR and intronic PAS reference region at once
refUTRraw=refUTRraw[which(refUTRraw$Chrom=='chr19'),]
dfIPAraw=dfIPA[which(dfIPA$Chrom=='chr19'),]
dfLEraw=dfLE[which(dfLE$Chrom=='chr19'),]   
PASREF=REF4PAS(refUTRraw,dfIPAraw,dfLEraw)
UTRdbraw=PASREF$UTRdbraw
dfIPA=PASREF$dfIPA
dfLE=PASREF$dfLE # Building 3’UTR PAS and IPA reference using GTF files
download.file(url='ftp://ftp.ensembl.org/pub/release-99/gtf/mus_musculus/Mus_musculus.GRCm38.99.gtf.gz',destfile='Mus_musculus.GRCm38.99.gtf.gz')           
GTFfile="Mus_musculus.GRCm38.99.gtf.gz" 
PASREFraw=PAS2GEF(GTFfile)  
refUTRraw=PASREFraw$refUTRraw
dfIPAraw=PASREFraw$dfIPA
dfLEraw=PASREFraw$dfLE
PASREF=REF4PAS(refUTRraw,dfIPAraw,dfLEraw)# Building aUTR and cUTR references
refUTRraw=refUTRraw[which(refUTRraw$Chrom=='chr19'),]
UTRdbraw=REF3UTR(refUTRraw)# Calculation of relative expression
DFUTRraw=PASEXP_3UTR(UTRdbraw, flsall, Strandtype="forward")# Building intronic polyA references
URL="https://github.com/RJWANGbioinfo/PAS_reference_RData/blob/master/"
file="mm9_REF.RData"
source_data(paste0(URL,file,"?raw=True"))# Calculation of relative expression
dfIPA=dfIPA[which(dfIPA$Chrom=='chr19'),]
dfLE=dfLE[which(dfLE$Chrom=='chr19'),]
IPA_OUTraw=PASEXP_IPA(dfIPA, dfLE, flsall, Strandtype="forward", nts=1)# Significance analysis of APA events
sampleTable1 = data.frame(samplename = c(names(flsall)),condition = c(rep("NT",3),rep("KD",3)))# Significantly regulated APA in 3’UTRs
test_3UTRsing=APAdiff(sampleTable2,DFUTRraw, conKET='NT',trtKEY='KD',PAS='3UTR',CUTreads=0,p_adjust_methods="fdr")
# Visualization of analysis results
APAVolcano(test_3UTRsing, PAS='3UTR', Pcol = "pvalue", top=5, main='3UTR APA')

实战案例

数据

下列样本存成bam_file.tsv

SampleID	BamPath
SRR316184	/Library/Frameworks/R.framework/Versions/4.1/Resources/library/TBX20BamSubset/extdata/SRR316184.bam
SRR316185	/Library/Frameworks/R.framework/Versions/4.1/Resources/library/TBX20BamSubset/extdata/SRR316185.bam
SRR316186	/Library/Frameworks/R.framework/Versions/4.1/Resources/library/TBX20BamSubset/extdata/SRR316186.bam
SRR316187	/Library/Frameworks/R.framework/Versions/4.1/Resources/library/TBX20BamSubset/extdata/SRR316187.bam
SRR316188	/Library/Frameworks/R.framework/Versions/4.1/Resources/library/TBX20BamSubset/extdata/SRR316188.bam
SRR316189	/Library/Frameworks/R.framework/Versions/4.1/Resources/library/TBX20BamSubset/extdata/SRR316189.bam

脚本

下列代码存成APAlyzer_Expression.R

suppressPackageStartupMessages({ library(dplyr)library(tibble)library(optparse)library(data.table)library(APAlyzer)library(TBX20BamSubset)library(Rsamtools)
})option_list <- list(make_option(c("-b", "--bam"), type = "character",help = "bam csv file (1st column: sampleID; 2nd: bam path)", metavar = "character"),make_option(c("-r", "--reference"), type = "character", # RData/gtfhelp = "genomic reference type", metavar = "character"),    make_option(c("-g", "--genome"), type = "character",help = "genomic reference file", metavar = "character"), make_option(c("-c", "--chromosome"), type = "character",default = "all", # chr19help = "chromosome to be selected", metavar = "character"),  make_option(c("-e", "--expression"), type = "character", default = "all", # 3UTR/IPAhelp = "APA expression: 3UTR and intronic APA", metavar = "character"),  make_option(c("-o", "--out"), type = "character",help = "output file path", metavar = "character")
)opt_parser <- OptionParser(option_list = option_list)
opt <- parse_args(opt_parser)# input parameters
bam_path <- opt$bam
ref_type <- opt$reference
ref_path <- opt$genome
chrom <- opt$chromosome
expr_type <- opt$expression
dir <- opt$out# bam_path <- "bam_file.tsv"
# ref_type <- "RData"
# ref_path <- "mm9_REF.RData"
# chrom <- "chr19"
# expr_type <- "3UTR"
# dir <- "result"# step1: bam file
bam_vector <- read.table("bam_file.tsv", header = TRUE)
bam_file <- bam_vector$BamPath
names(bam_file) <- bam_vector$SampleID# step2: genomic reference
if (ref_type == "RData") {# data from built referencerequire(repmis)URL <- "https://github.com/RJWANGbioinfo/PAS_reference_RData/blob/master/"source_data(paste0(URL, ref_path, "?raw=True"))if (ref_path == "mm9_REF.RData") {refUTRraw_temp <- refUTRrawdfIPAraw_temp <- dfIPAdfLEraw_temp <- dfLE} else if (ref_path == "hg19_REF.RData") {refUTRraw_temp <- refUTRraw_hg19dfIPAraw_temp <- dfIPA_hg19dfLEraw_temp <- dfLE_hg19}} else if (ref_type == "gtf") {# building reference from gtf filePASREFraw <- PAS2GEF(ref_path)  refUTRraw_temp <- PASREFraw$refUTRrawdfIPAraw_temp <- PASREFraw$dfIPAdfLEraw_temp <- PASREFraw$dfLE
}# step3: whether to choose chromosome
if (chrom == "all") {UTRdbraw <- refUTRraw_tempdfIPAraw <- dfIPAraw_tempdfLEraw <- dfLEraw_temp   
} else {# multiple chromosome or notif (length(grep(":", chrom)) > 0) {chroms <- unlist(strsplit(chrom, ":"))} else {chroms <- chrom}UTRdbraw <- refUTRraw_temp[which(refUTRraw_temp$Chrom %in% chroms), ]dfIPAraw <- dfIPAraw_temp[which(dfIPAraw_temp$Chrom %in% chroms), ]dfLEraw <- dfLEraw_temp[which(dfLEraw_temp$Chrom %in% chroms), ]
}
## aUTR cUTR
PASREF_temp <- REF4PAS(UTRdbraw, dfIPAraw, dfLEraw)
UTRdb <- PASREF_temp$UTRdbraw
dfIPA <- PASREF_temp$dfIPA
dfLE <- PASREF_temp$dfLE  # step4: APA expression (3UTR and IPA)
if (expr_type == "all") {# 3UTRUTR_APA_OUT <- PASEXP_3UTR(UTRdb, bam_file, Strandtype = "forward")# IPAIPA_OUT <- PASEXP_IPA(dfIPA, dfLE, bam_file, Strandtype = "invert", nts = 4)final_OUT <- list(UTR = UTR_APA_OUT,IPA = IPA_OUT)
} else if (expr_type == "3UTR") { # 3UTRfinal_OUT <- PASEXP_3UTR(UTRdb, bam_file, Strandtype = "forward")  
} else if (expr_type == "IPA") { final_OUT <- PASEXP_IPA(dfIPA, dfLE, bam_file, Strandtype = "invert", nts = 4)
}# step5: output
if (!dir.exists(dir)) {dir.create(dir, recursive = TRUE)
}if (!is.data.frame(final_OUT)) {file_name <- paste0(dir, "/APA_Expr_", expr_type, ".RDS")saveRDS(final_OUT, file_name, compress = TRUE)
} else {file_name <- paste0(dir, "/APA_Expr_", expr_type, ".tsv")write.table(final_OUT, file_name, quote = F, row.names = F, sep = "\t")
}print("Program Ended without Problems")

运行

在命令行模式下运行该命令

Rscript APAlyzer_Expression.R \-b bam_file.tsv \-r RData \-g mm9_REF.RData \-c chr19 \-e 3UTR \-o result