##############################################################################
##############################################################################
###
### Running command:
###
### C:\Users\biocbuild\bbs-3.9-bioc\R\bin\R.exe CMD check --force-multiarch --install=check:exomePeak.install-out.txt --library=C:\Users\biocbuild\bbs-3.9-bioc\R\library --no-vignettes --timings exomePeak_2.18.0.tar.gz
###
##############################################################################
##############################################################################
* using log directory 'C:/Users/biocbuild/bbs-3.9-bioc/meat/exomePeak.Rcheck'
* using R version 3.6.1 (2019-07-05)
* using platform: x86_64-w64-mingw32 (64-bit)
* using session charset: ISO8859-1
* using option '--no-vignettes'
* checking for file 'exomePeak/DESCRIPTION' ... OK
* checking extension type ... Package
* this is package 'exomePeak' version '2.18.0'
* checking package namespace information ... OK
* checking package dependencies ... OK
* checking if this is a source package ... OK
* checking if there is a namespace ... OK
* checking for hidden files and directories ... OK
* checking for portable file names ... OK
* checking whether package 'exomePeak' can be installed ... OK
* checking installed package size ... OK
* checking package directory ... OK
* checking 'build' directory ... OK
* checking DESCRIPTION meta-information ... OK
* checking top-level files ... OK
* checking for left-over files ... OK
* checking index information ... OK
* checking package subdirectories ... OK
* checking R files for non-ASCII characters ... OK
* checking R files for syntax errors ... OK
* loading checks for arch 'i386'
** checking whether the package can be loaded ... OK
** checking whether the package can be loaded with stated dependencies ... OK
** checking whether the package can be unloaded cleanly ... OK
** checking whether the namespace can be loaded with stated dependencies ... OK
** checking whether the namespace can be unloaded cleanly ... OK
* loading checks for arch 'x64'
** checking whether the package can be loaded ... OK
** checking whether the package can be loaded with stated dependencies ... OK
** checking whether the package can be unloaded cleanly ... OK
** checking whether the namespace can be loaded with stated dependencies ... OK
** checking whether the namespace can be unloaded cleanly ... OK
* checking dependencies in R code ... NOTE
Packages in Depends field not imported from:
'GenomicAlignments' 'GenomicFeatures' 'Rsamtools' 'rtracklayer'
These packages need to be imported from (in the NAMESPACE file)
for when this namespace is loaded but not attached.
* checking S3 generic/method consistency ... OK
* checking replacement functions ... OK
* checking foreign function calls ... OK
* checking R code for possible problems ... NOTE
.bed2grangeslist: no visible global function definition for 'head'
.bed2grangeslist: no visible global function definition for 'makeTxDb'
.bed2grangeslist: no visible global function definition for 'exonsBy'
.exomepeak_notification: no visible global function definition for
'mcols<-'
.get.bam.chrs: no visible global function definition for 'indexBam'
.get.bam.chrs: no visible global function definition for 'ScanBamParam'
.get.bam.chrs: no visible global function definition for 'scanBam'
.get.bam.read.length: no visible global function definition for
'indexBam'
.get.bam.read.length: no visible global function definition for
'ScanBamParam'
.get.bam.read.length: no visible global function definition for
'scanBam'
.get.bam.read.length: no visible global function definition for
'median'
.get.check.points.reads.count: no visible global function definition
for 'IRangesList'
.get.check.points.reads.count: no visible global function definition
for 'IRanges'
.get.check.points.reads.count: no visible global function definition
for 'ScanBamParam'
.get.check.points.reads.count: no visible global function definition
for 'scanBam'
.read.gtf: no visible global function definition for 'makeTxDbFromUCSC'
.read.gtf: no visible global function definition for 'makeTxDbFromGFF'
.read.gtf: no visible global function definition for 'columns'
.read.gtf: no visible global function definition for 'keys'
.read.gtf: no visible global function definition for 'select'
.readTxDb2: no visible global function definition for
'makeTxDbFromUCSC'
.readTxDb2: no visible global function definition for 'makeTxDbFromGFF'
.report.diff.peak.based.on.result: no visible global function
definition for 'write.table'
.report.peak.based.on.result: no visible global function definition for
'write.table'
.xls2Grangeslist: no visible global function definition for
'read.table'
.xls2Grangeslist: no visible global function definition for 'head'
.xls2Grangeslist: no visible global function definition for 'makeTxDb'
.xls2Grangeslist: no visible global function definition for 'exonsBy'
.xls2Grangeslist: no visible global function definition for 'mcols<-'
RMT: no visible global function definition for 'read.table'
RMT: no visible global function definition for 'write.table'
RMT: no visible global function definition for 'findOverlaps'
RMT: no visible global function definition for 'queryHits'
RMT: no visible global function definition for 'subjectHits'
RMT: no visible global function definition for 'ranges'
RMT: no visible global function definition for 'exonsBy'
RMT: no visible global function definition for 'readGAlignments'
RMT: no visible global function definition for 'ScanBamParam'
RMT: no visible global function definition for 'scanBam'
RMT: no visible global function definition for 'countOverlaps'
RMT: no visible global function definition for 'width'
RMT: no visible global function definition for 'var'
bltest: no visible global function definition for 'pchisq'
bltest: no visible global function definition for 'p.adjust'
ctest: no visible global function definition for 'pbinom'
ctest: no visible global function definition for 'p.adjust'
exomepeak: no visible global function definition for 'indexBam'
rhtest: no visible global function definition for 'phyper'
rhtest: no visible global function definition for 'p.adjust'
Undefined global functions or variables:
IRanges IRangesList ScanBamParam columns countOverlaps exonsBy
findOverlaps head indexBam keys makeTxDb makeTxDbFromGFF
makeTxDbFromUCSC mcols<- median p.adjust pbinom pchisq phyper
queryHits ranges read.table readGAlignments scanBam select
subjectHits var width write.table
Consider adding
importFrom("stats", "median", "p.adjust", "pbinom", "pchisq", "phyper",
"var")
importFrom("utils", "head", "read.table", "write.table")
to your NAMESPACE file.
* checking Rd files ... NOTE
prepare_Rd: RMT.Rd:60-62: Dropping empty section \details
prepare_Rd: RMT.Rd:74-76: Dropping empty section \note
prepare_Rd: RMT.Rd:80-82: Dropping empty section \seealso
* checking Rd metadata ... OK
* checking Rd cross-references ... OK
* checking for missing documentation entries ... OK
* checking for code/documentation mismatches ... OK
* checking Rd \usage sections ... OK
* checking Rd contents ... OK
* checking for unstated dependencies in examples ... OK
* checking files in 'vignettes' ... OK
* checking examples ...
** running examples for arch 'i386' ... ERROR
Running examples in 'exomePeak-Ex.R' failed
The error most likely occurred in:
> base::assign(".ptime", proc.time(), pos = "CheckExEnv")
> ### Name: RMT
> ### Title: Extract a combinatorial RNA methylome from multiple biological
> ### conditions or replicates.
> ### Aliases: RMT
> ### Keywords: MeRIP-Seq
>
> ### ** Examples
>
> # the RESPECT R-package has two main functions:
> # 1. exome-based peak
> # 2. compute reads counts, rpkm and fold change
> # please feel free to contact [email protected] for any questions
>
> # load library and specify the parameters
> gene_anno_gtf <- system.file("extdata", "example.gtf", package="exomePeak")
> f1 <- system.file("extdata", "IP1.bam", package="exomePeak")
> f2 <- system.file("extdata", "IP2.bam", package="exomePeak")
> f3 <- system.file("extdata", "IP3.bam", package="exomePeak")
> f4 <- system.file("extdata", "IP4.bam", package="exomePeak")
> ip_bam <- c(f1,f2,f3,f4)
> f1 <- system.file("extdata", "Input1.bam", package="exomePeak")
> f2 <- system.file("extdata", "Input2.bam", package="exomePeak")
> f3 <- system.file("extdata", "Input3.bam", package="exomePeak")
> input_bam <- c(f1,f2,f3)
> input2bam <- list(4)
> input2bam[[1]] <- c(1,2) # the 1st ip sample uses the 1st & 2nd input samples as control
> input2bam[[2]] <- c(1,2) # the 2nd ip sample uses the 1st & 2nd input samples as control
> input2bam[[3]] <- c(1,2) # the 3rd ip sample uses the 1st & 2nd input samples as control
> input2bam[[4]] <- c(3) # the 4th ip sample uses the 3rd input sample as control
>
> # Extract the combinatorial RNA methylome
> # unfortunately, this function has not been optimized for parallel processing.
> # This part will take a really long time on real data
> RMT(INPUT_BAM=input_bam, IP_BAM=ip_bam, INPUT2IP=input2bam, GENE_ANNO_GTF=gene_anno_gtf)
[1] "Processing IP sample 1"
Import genomic features from the file as a GRanges object ... OK
Prepare the 'metadata' data frame ... OK
Make the TxDb object ... OK
'select()' returned 1:many mapping between keys and columns
[1] "Divide transcriptome into chr-gene-batch sections ..."
[1] "Get Reads Count ..."
[1] "This step may take a few hours ..."
[1] "100 %"
[1] "Get all the peaks ..."
[1] "Get the consistent peaks ..."
[1] "---------------------------------"
[1] "The bam files used:"
[1] "1 IP replicate(s)"
[1] "2 Input replicate(s)"
[1] "---------------------------------"
[1] "Peak calling result: "
[1] "15 peaks detected on merged data."
[1] "Please check 'peak.bed/xls' under C:/Users/biocbuild/bbs-3.9-bioc/meat/exomePeak.Rcheck/examples_i386/RMT_Result/temp/IP_rep_1"
[1] "15 consistent peaks detected on every replicates. (Recommended list)"
[1] "Please check 'con_peak.bed/xls' under C:/Users/biocbuild/bbs-3.9-bioc/meat/exomePeak.Rcheck/examples_i386/RMT_Result/temp/IP_rep_1"
[1] "Processing IP sample 2"
Import genomic features from the file as a GRanges object ... OK
Prepare the 'metadata' data frame ... OK
Make the TxDb object ... OK
'select()' returned 1:many mapping between keys and columns
[1] "Divide transcriptome into chr-gene-batch sections ..."
[1] "Get Reads Count ..."
[1] "This step may take a few hours ..."
[1] "100 %"
[1] "Get all the peaks ..."
[1] "Get the consistent peaks ..."
[1] "---------------------------------"
[1] "The bam files used:"
[1] "1 IP replicate(s)"
[1] "2 Input replicate(s)"
[1] "---------------------------------"
[1] "Peak calling result: "
[1] "19 peaks detected on merged data."
[1] "Please check 'peak.bed/xls' under C:/Users/biocbuild/bbs-3.9-bioc/meat/exomePeak.Rcheck/examples_i386/RMT_Result/temp/IP_rep_2"
[1] "19 consistent peaks detected on every replicates. (Recommended list)"
[1] "Please check 'con_peak.bed/xls' under C:/Users/biocbuild/bbs-3.9-bioc/meat/exomePeak.Rcheck/examples_i386/RMT_Result/temp/IP_rep_2"
[1] "Processing IP sample 3"
Import genomic features from the file as a GRanges object ... OK
Prepare the 'metadata' data frame ... OK
Make the TxDb object ... OK
'select()' returned 1:many mapping between keys and columns
[1] "Divide transcriptome into chr-gene-batch sections ..."
[1] "Get Reads Count ..."
[1] "This step may take a few hours ..."
[1] "100 %"
[1] "Get all the peaks ..."
[1] "Get the consistent peaks ..."
[1] "---------------------------------"
[1] "The bam files used:"
[1] "1 IP replicate(s)"
[1] "2 Input replicate(s)"
[1] "---------------------------------"
[1] "Peak calling result: "
[1] "14 peaks detected on merged data."
[1] "Please check 'peak.bed/xls' under C:/Users/biocbuild/bbs-3.9-bioc/meat/exomePeak.Rcheck/examples_i386/RMT_Result/temp/IP_rep_3"
[1] "14 consistent peaks detected on every replicates. (Recommended list)"
[1] "Please check 'con_peak.bed/xls' under C:/Users/biocbuild/bbs-3.9-bioc/meat/exomePeak.Rcheck/examples_i386/RMT_Result/temp/IP_rep_3"
[1] "Processing IP sample 4"
Import genomic features from the file as a GRanges object ... OK
Prepare the 'metadata' data frame ... OK
Make the TxDb object ... OK
'select()' returned 1:many mapping between keys and columns
[1] "Divide transcriptome into chr-gene-batch sections ..."
[1] "Get Reads Count ..."
[1] "This step may take a few hours ..."
[1] "100 %"
[1] "Get all the peaks ..."
[1] "Get the consistent peaks ..."
[1] "---------------------------------"
[1] "The bam files used:"
[1] "1 IP replicate(s)"
[1] "1 Input replicate(s)"
[1] "---------------------------------"
[1] "Peak calling result: "
[1] "5 peaks detected on merged data."
[1] "Please check 'peak.bed/xls' under C:/Users/biocbuild/bbs-3.9-bioc/meat/exomePeak.Rcheck/examples_i386/RMT_Result/temp/IP_rep_4"
[1] "5 consistent peaks detected on every replicates. (Recommended list)"
[1] "Please check 'con_peak.bed/xls' under C:/Users/biocbuild/bbs-3.9-bioc/meat/exomePeak.Rcheck/examples_i386/RMT_Result/temp/IP_rep_4"
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
exomePeak
--- call from context ---
RMT(INPUT_BAM = input_bam, IP_BAM = ip_bam, INPUT2IP = input2bam,
GENE_ANNO_GTF = gene_anno_gtf)
--- call from argument ---
if (x < y) {
num[matrix_overlap[i, 1][1][[1]]] <- matrix_overlap[i, 1][1][[1]]
} else if (x > y) {
num[matrix_overlap[i, 2][1][[1]]] <- matrix_overlap[i, 2][1][[1]]
} else {
num[matrix_overlap[i, 1][1][[1]]] <- matrix_overlap[i, 1][1][[1]]
num[matrix_overlap[i, 2][1][[1]]] <- matrix_overlap[i, 2][1][[1]]
}
--- R stacktrace ---
where 1: RMT(INPUT_BAM = input_bam, IP_BAM = ip_bam, INPUT2IP = input2bam,
GENE_ANNO_GTF = gene_anno_gtf)
--- value of length: 4 type: logical ---
[1] FALSE TRUE FALSE TRUE
--- function from context ---
function (INPUT_BAM, IP_BAM, INPUT2IP = NA, GENE_ANNO_GTF = NA,
GENOME = NA, UCSC_TABLE_NAME = "knownGene", TXDB = NA, EXOME_OUTPUT_DIR = NA,
RMT_OUTPUT_DIR = NA, RMT_EXPERIMENT_NAME = "RMT_Result")
{
PARAMETERS = list()
PARAMETERS$INPUT_BAM = INPUT_BAM
PARAMETERS$IP_BAM = IP_BAM
PARAMETERS$INPUT2IP = INPUT2IP
PARAMETERS$EXOME_OUTPUT_DIR = EXOME_OUTPUT_DIR
PARAMETERS$GO_OUTPUT_DIR = RMT_OUTPUT_DIR
PARAMETERS$GO_EXPERIMENT_NAME = RMT_EXPERIMENT_NAME
PARAMETERS$GENE_ANNO_GTF = GENE_ANNO_GTF
PARAMETERS$GENOME = GENOME
PARAMETERS$UCSC_TABLE_NAME = UCSC_TABLE_NAME
PARAMETERS$TXDB = TXDB
if (is.na(PARAMETERS$EXOME_OUTPUT_DIR)) {
PARAMETERS$EXOME_OUTPUT_DIR <- getwd()
}
if (is.na(PARAMETERS$GO_OUTPUT_DIR)) {
PARAMETERS$GO_OUTPUT_DIR <- getwd()
}
input_length <- length(PARAMETERS$INPUT_BAM)
ip_length <- length(PARAMETERS$IP_BAM)
all_filepath <- vector()
temp1 <- is.na(PARAMETERS$INPUT2IP)
flag <- temp1[1]
if (flag) {
for (i in 1:ip_length) {
experiment_name = paste(PARAMETERS$GO_EXPERIMENT_NAME,
"/temp/IP_rep_", i, sep = "")
all_filepath[i] = PARAMETERS$IP_BAM[i]
print(paste("Processing IP sample", i))
temp <- exomepeak(GENE_ANNO_GTF = GENE_ANNO_GTF,
GENOME = GENOME, UCSC_TABLE_NAME = UCSC_TABLE_NAME,
TXDB = TXDB, IP_BAM = c(PARAMETERS$IP_BAM[i]),
INPUT_BAM = c(PARAMETERS$INPUT_BAM[i]), OUTPUT_DIR = PARAMETERS$EXOME_OUTPUT_DIR,
POISSON_MEAN_RATIO = 1, EXPERIMENT_NAME = experiment_name)
}
}
else {
for (i in 1:ip_length) {
experiment_name = paste(PARAMETERS$GO_EXPERIMENT_NAME,
"/temp/IP_rep_", i, sep = "")
input_bam = PARAMETERS$INPUT_BAM[PARAMETERS$INPUT2IP[[i]][1]]
if (length(PARAMETERS$INPUT2IP[[i]]) > 1) {
for (j in 2:length(PARAMETERS$INPUT2IP[[i]])) {
input_bam = c(input_bam, c(PARAMETERS$INPUT_BAM[PARAMETERS$INPUT2IP[[i]][j]]))
}
}
all_filepath[i] <- paste(PARAMETERS$EXOME_OUTPUT_DIR,
experiment_name, sep = "/")
print(paste("Processing IP sample", i))
temp <- exomepeak(GENE_ANNO_GTF = GENE_ANNO_GTF,
GENOME = GENOME, UCSC_TABLE_NAME = UCSC_TABLE_NAME,
TXDB = TXDB, IP_BAM = c(PARAMETERS$IP_BAM[i]),
INPUT_BAM = input_bam, OUTPUT_DIR = PARAMETERS$EXOME_OUTPUT_DIR,
POISSON_MEAN_RATIO = 1, EXPERIMENT_NAME = experiment_name)
}
}
for (i in 1:length(all_filepath)) {
path <- paste(all_filepath[i], "peak.xls", sep = "/")
bam_file <- read.table(path, sep = "\t", header = FALSE)
if (i == 1) {
all_peak <- bam_file
}
else {
all_peak = rbind(all_peak, bam_file[2:nrow(bam_file),
])
}
}
dir.create(paste(PARAMETERS$GO_OUTPUT_DIR, PARAMETERS$GO_EXPERIMENT_NAME,
sep = "/"), recursive = TRUE, showWarnings = FALSE)
dir = paste(PARAMETERS$GO_OUTPUT_DIR, PARAMETERS$GO_EXPERIMENT_NAME,
sep = "/")
write.table(all_peak, paste(dir, "all_peak.xls", sep = "/"),
sep = "\t", quote = FALSE, col.names = FALSE, row.names = FALSE)
filepath = paste(dir, "all_peak.xls", sep = "/")
matrix_peak_read = read.table(filepath, header = TRUE, stringsAsFactors = FALSE)
ori_peak <- .xls2Grangeslist(filepath)
overlaps <- findOverlaps(ori_peak, ori_peak)
matrix_overlap <- matrix(0, length(overlaps), 2)
matrix_overlap[, 1] <- queryHits(overlaps)
matrix_overlap[, 2] <- subjectHits(overlaps)
num <- c(rep(0, length(ori_peak)))
for (i in 1:nrow(matrix_overlap)) {
if (matrix_overlap[i, 1][1][[1]] != matrix_overlap[i,
2][1][[1]]) {
x <- ranges(ori_peak[matrix_overlap[i, 1]])[[1]]@width
y <- ranges(ori_peak[matrix_overlap[i, 2]])[[1]]@width
if (x < y) {
num[matrix_overlap[i, 1][1][[1]]] <- matrix_overlap[i,
1][1][[1]]
}
else if (x > y) {
num[matrix_overlap[i, 2][1][[1]]] <- matrix_overlap[i,
2][1][[1]]
}
else {
num[matrix_overlap[i, 1][1][[1]]] <- matrix_overlap[i,
1][1][[1]]
num[matrix_overlap[i, 2][1][[1]]] <- matrix_overlap[i,
2][1][[1]]
}
}
}
peak <- ori_peak[num == 0]
tab <- which(num == 0)
write.table(matrix_peak_read[tab, ], paste(dir, "merged_peak.xls",
sep = "/"), sep = "\t", quote = FALSE, col.names = TRUE,
row.names = FALSE)
write.table(matrix_peak_read[tab, 1:12], paste(dir, "merged_peak.bed",
sep = "/"), sep = "\t", quote = FALSE, col.names = FALSE,
row.names = FALSE)
txdb <- .readTxDb2(PARAMETERS)
exonRanges <- exonsBy(txdb, "gene")
rpkm_input <- matrix(0, length(peak), length(PARAMETERS$INPUT_BAM))
rpkm_ip <- matrix(0, length(peak), length(PARAMETERS$IP_BAM))
count_input <- matrix(0, length(peak), length(PARAMETERS$INPUT_BAM))
count_ip <- matrix(0, length(peak), length(PARAMETERS$IP_BAM))
for (i in 1:length(PARAMETERS$INPUT_BAM)) {
filepath <- PARAMETERS$INPUT_BAM[i]
aligns <- readGAlignments(filepath)
para <- ScanBamParam(what = "mapq")
mapq <- scanBam(filepath, param = para)[[1]][[1]]
mapq[is.na(mapq)] <- 255
ID_keep <- (mapq > 30)
filtered <- aligns[ID_keep]
id <- countOverlaps(filtered, exonRanges)
transcriptome_filtered_aligns <- filtered[id > 0]
counts <- countOverlaps(peak, transcriptome_filtered_aligns)
count_input[, i] <- counts
numBases <- sum(width(peak))
geneLengthsInKB <- numBases/1000
millionsMapped <- length(transcriptome_filtered_aligns)/1e+06
rpm <- counts/millionsMapped
rpkm_input[, i] <- rpm/geneLengthsInKB
}
write.table(count_input, paste(dir, "count_input.xls", sep = "/"),
sep = "\t", quote = FALSE, col.names = TRUE, row.names = FALSE)
write.table(rpkm_input, paste(dir, "rpkm_input.xls", sep = "/"),
sep = "\t", quote = FALSE, col.names = TRUE, row.names = FALSE)
for (i in 1:length(PARAMETERS$IP_BAM)) {
filepath <- PARAMETERS$IP_BAM[i]
aligns <- readGAlignments(filepath)
para <- ScanBamParam(what = "mapq")
mapq <- scanBam(filepath, param = para)[[1]][[1]]
mapq[is.na(mapq)] <- 255
ID_keep <- (mapq > 30)
filtered <- aligns[ID_keep]
id <- countOverlaps(filtered, exonRanges)
transcriptome_filtered_aligns <- filtered[id > 0]
counts <- countOverlaps(peak, transcriptome_filtered_aligns)
count_ip[, i] <- counts
numBases <- sum(width(peak))
geneLengthsInKB <- numBases/1000
millionsMapped <- length(transcriptome_filtered_aligns)/1e+06
rpm <- counts/millionsMapped
rpkm_ip[, i] <- rpm/geneLengthsInKB
}
write.table(count_ip, paste(dir, "count_ip.xls", sep = "/"),
sep = "\t", quote = FALSE, col.names = TRUE, row.names = FALSE)
write.table(rpkm_ip, paste(dir, "rpkm_ip.xls", sep = "/"),
sep = "\t", quote = FALSE, col.names = TRUE, row.names = FALSE)
matrix_log2_fc = matrix(0, length(peak), ip_length)
if (is.na(PARAMETERS$INPUT2IP)) {
for (i in 1:ip_length) {
matrix_log2_fc[, i] <- log2(rpkm_ip[, i] + 0.01) -
log2(rpkm_input[, i] + 0.01)
}
}
else {
for (i in 1:ip_length) {
if (length(PARAMETERS$INPUT2IP[[i]]) == 1) {
matrix_log2_fc[, i] <- log2(rpkm_ip[, i] + 0.01) -
log2(rpkm_input[, PARAMETERS$INPUT2IP[[i]]] +
0.01)
}
else {
sum = matrix(0, length(peak), 1)
for (j in 1:length(PARAMETERS$INPUT2IP[i])) {
sum = sum + rpkm_input[, PARAMETERS$INPUT2IP[[i]][j]]
}
matrix_log2_fc[, i] <- log2(rpkm_ip[, i] + 0.01) -
log2(sum + 0.01)
}
}
}
v <- matrix_log2_fc[, 1]
for (i in 1:nrow(matrix_log2_fc)) {
v[i] <- var(matrix_log2_fc[i, ])
}
matrix_log2_fc <- cbind(tab, matrix_log2_fc[1:nrow(matrix_log2_fc),
], matrix_peak_read[tab, ])
write.table(matrix_log2_fc, paste(dir, "all_information.xls",
sep = "/"), sep = "\t", quote = FALSE, col.names = TRUE,
row.names = FALSE)
file.remove(paste(dir, "all_peak.xls", sep = "/"))
file.remove(paste(dir, "all_information.xls", sep = "/"))
print("***************************")
print("***************************")
print(paste("The combinatorial RNA methylome is generated under:",
dir))
print("Including:")
print("1. Merged Peaks")
print("2. Merged Peaks in BED format")
print("3. Reads count for every peak in every bam file")
print("4. RPKM for every peak in every bam file")
}
--- function search by body ---
Function RMT in namespace exomePeak has this body.
----------- END OF FAILURE REPORT --------------
Fatal error: the condition has length > 1
** running examples for arch 'x64' ... ERROR
Running examples in 'exomePeak-Ex.R' failed
The error most likely occurred in:
> base::assign(".ptime", proc.time(), pos = "CheckExEnv")
> ### Name: RMT
> ### Title: Extract a combinatorial RNA methylome from multiple biological
> ### conditions or replicates.
> ### Aliases: RMT
> ### Keywords: MeRIP-Seq
>
> ### ** Examples
>
> # the RESPECT R-package has two main functions:
> # 1. exome-based peak
> # 2. compute reads counts, rpkm and fold change
> # please feel free to contact [email protected] for any questions
>
> # load library and specify the parameters
> gene_anno_gtf <- system.file("extdata", "example.gtf", package="exomePeak")
> f1 <- system.file("extdata", "IP1.bam", package="exomePeak")
> f2 <- system.file("extdata", "IP2.bam", package="exomePeak")
> f3 <- system.file("extdata", "IP3.bam", package="exomePeak")
> f4 <- system.file("extdata", "IP4.bam", package="exomePeak")
> ip_bam <- c(f1,f2,f3,f4)
> f1 <- system.file("extdata", "Input1.bam", package="exomePeak")
> f2 <- system.file("extdata", "Input2.bam", package="exomePeak")
> f3 <- system.file("extdata", "Input3.bam", package="exomePeak")
> input_bam <- c(f1,f2,f3)
> input2bam <- list(4)
> input2bam[[1]] <- c(1,2) # the 1st ip sample uses the 1st & 2nd input samples as control
> input2bam[[2]] <- c(1,2) # the 2nd ip sample uses the 1st & 2nd input samples as control
> input2bam[[3]] <- c(1,2) # the 3rd ip sample uses the 1st & 2nd input samples as control
> input2bam[[4]] <- c(3) # the 4th ip sample uses the 3rd input sample as control
>
> # Extract the combinatorial RNA methylome
> # unfortunately, this function has not been optimized for parallel processing.
> # This part will take a really long time on real data
> RMT(INPUT_BAM=input_bam, IP_BAM=ip_bam, INPUT2IP=input2bam, GENE_ANNO_GTF=gene_anno_gtf)
[1] "Processing IP sample 1"
Import genomic features from the file as a GRanges object ... OK
Prepare the 'metadata' data frame ... OK
Make the TxDb object ... OK
'select()' returned 1:many mapping between keys and columns
[1] "Divide transcriptome into chr-gene-batch sections ..."
[1] "Get Reads Count ..."
[1] "This step may take a few hours ..."
[1] "100 %"
[1] "Get all the peaks ..."
[1] "Get the consistent peaks ..."
[1] "---------------------------------"
[1] "The bam files used:"
[1] "1 IP replicate(s)"
[1] "2 Input replicate(s)"
[1] "---------------------------------"
[1] "Peak calling result: "
[1] "15 peaks detected on merged data."
[1] "Please check 'peak.bed/xls' under C:/Users/biocbuild/bbs-3.9-bioc/meat/exomePeak.Rcheck/examples_x64/RMT_Result/temp/IP_rep_1"
[1] "15 consistent peaks detected on every replicates. (Recommended list)"
[1] "Please check 'con_peak.bed/xls' under C:/Users/biocbuild/bbs-3.9-bioc/meat/exomePeak.Rcheck/examples_x64/RMT_Result/temp/IP_rep_1"
[1] "Processing IP sample 2"
Import genomic features from the file as a GRanges object ... OK
Prepare the 'metadata' data frame ... OK
Make the TxDb object ... OK
'select()' returned 1:many mapping between keys and columns
[1] "Divide transcriptome into chr-gene-batch sections ..."
[1] "Get Reads Count ..."
[1] "This step may take a few hours ..."
[1] "100 %"
[1] "Get all the peaks ..."
[1] "Get the consistent peaks ..."
[1] "---------------------------------"
[1] "The bam files used:"
[1] "1 IP replicate(s)"
[1] "2 Input replicate(s)"
[1] "---------------------------------"
[1] "Peak calling result: "
[1] "19 peaks detected on merged data."
[1] "Please check 'peak.bed/xls' under C:/Users/biocbuild/bbs-3.9-bioc/meat/exomePeak.Rcheck/examples_x64/RMT_Result/temp/IP_rep_2"
[1] "19 consistent peaks detected on every replicates. (Recommended list)"
[1] "Please check 'con_peak.bed/xls' under C:/Users/biocbuild/bbs-3.9-bioc/meat/exomePeak.Rcheck/examples_x64/RMT_Result/temp/IP_rep_2"
[1] "Processing IP sample 3"
Import genomic features from the file as a GRanges object ... OK
Prepare the 'metadata' data frame ... OK
Make the TxDb object ... OK
'select()' returned 1:many mapping between keys and columns
[1] "Divide transcriptome into chr-gene-batch sections ..."
[1] "Get Reads Count ..."
[1] "This step may take a few hours ..."
[1] "100 %"
[1] "Get all the peaks ..."
[1] "Get the consistent peaks ..."
[1] "---------------------------------"
[1] "The bam files used:"
[1] "1 IP replicate(s)"
[1] "2 Input replicate(s)"
[1] "---------------------------------"
[1] "Peak calling result: "
[1] "14 peaks detected on merged data."
[1] "Please check 'peak.bed/xls' under C:/Users/biocbuild/bbs-3.9-bioc/meat/exomePeak.Rcheck/examples_x64/RMT_Result/temp/IP_rep_3"
[1] "14 consistent peaks detected on every replicates. (Recommended list)"
[1] "Please check 'con_peak.bed/xls' under C:/Users/biocbuild/bbs-3.9-bioc/meat/exomePeak.Rcheck/examples_x64/RMT_Result/temp/IP_rep_3"
[1] "Processing IP sample 4"
Import genomic features from the file as a GRanges object ... OK
Prepare the 'metadata' data frame ... OK
Make the TxDb object ... OK
'select()' returned 1:many mapping between keys and columns
[1] "Divide transcriptome into chr-gene-batch sections ..."
[1] "Get Reads Count ..."
[1] "This step may take a few hours ..."
[1] "100 %"
[1] "Get all the peaks ..."
[1] "Get the consistent peaks ..."
[1] "---------------------------------"
[1] "The bam files used:"
[1] "1 IP replicate(s)"
[1] "1 Input replicate(s)"
[1] "---------------------------------"
[1] "Peak calling result: "
[1] "5 peaks detected on merged data."
[1] "Please check 'peak.bed/xls' under C:/Users/biocbuild/bbs-3.9-bioc/meat/exomePeak.Rcheck/examples_x64/RMT_Result/temp/IP_rep_4"
[1] "5 consistent peaks detected on every replicates. (Recommended list)"
[1] "Please check 'con_peak.bed/xls' under C:/Users/biocbuild/bbs-3.9-bioc/meat/exomePeak.Rcheck/examples_x64/RMT_Result/temp/IP_rep_4"
----------- FAILURE REPORT --------------
--- failure: the condition has length > 1 ---
--- srcref ---
:
--- package (from environment) ---
exomePeak
--- call from context ---
RMT(INPUT_BAM = input_bam, IP_BAM = ip_bam, INPUT2IP = input2bam,
GENE_ANNO_GTF = gene_anno_gtf)
--- call from argument ---
if (x < y) {
num[matrix_overlap[i, 1][1][[1]]] <- matrix_overlap[i, 1][1][[1]]
} else if (x > y) {
num[matrix_overlap[i, 2][1][[1]]] <- matrix_overlap[i, 2][1][[1]]
} else {
num[matrix_overlap[i, 1][1][[1]]] <- matrix_overlap[i, 1][1][[1]]
num[matrix_overlap[i, 2][1][[1]]] <- matrix_overlap[i, 2][1][[1]]
}
--- R stacktrace ---
where 1: RMT(INPUT_BAM = input_bam, IP_BAM = ip_bam, INPUT2IP = input2bam,
GENE_ANNO_GTF = gene_anno_gtf)
--- value of length: 4 type: logical ---
[1] FALSE TRUE FALSE TRUE
--- function from context ---
function (INPUT_BAM, IP_BAM, INPUT2IP = NA, GENE_ANNO_GTF = NA,
GENOME = NA, UCSC_TABLE_NAME = "knownGene", TXDB = NA, EXOME_OUTPUT_DIR = NA,
RMT_OUTPUT_DIR = NA, RMT_EXPERIMENT_NAME = "RMT_Result")
{
PARAMETERS = list()
PARAMETERS$INPUT_BAM = INPUT_BAM
PARAMETERS$IP_BAM = IP_BAM
PARAMETERS$INPUT2IP = INPUT2IP
PARAMETERS$EXOME_OUTPUT_DIR = EXOME_OUTPUT_DIR
PARAMETERS$GO_OUTPUT_DIR = RMT_OUTPUT_DIR
PARAMETERS$GO_EXPERIMENT_NAME = RMT_EXPERIMENT_NAME
PARAMETERS$GENE_ANNO_GTF = GENE_ANNO_GTF
PARAMETERS$GENOME = GENOME
PARAMETERS$UCSC_TABLE_NAME = UCSC_TABLE_NAME
PARAMETERS$TXDB = TXDB
if (is.na(PARAMETERS$EXOME_OUTPUT_DIR)) {
PARAMETERS$EXOME_OUTPUT_DIR <- getwd()
}
if (is.na(PARAMETERS$GO_OUTPUT_DIR)) {
PARAMETERS$GO_OUTPUT_DIR <- getwd()
}
input_length <- length(PARAMETERS$INPUT_BAM)
ip_length <- length(PARAMETERS$IP_BAM)
all_filepath <- vector()
temp1 <- is.na(PARAMETERS$INPUT2IP)
flag <- temp1[1]
if (flag) {
for (i in 1:ip_length) {
experiment_name = paste(PARAMETERS$GO_EXPERIMENT_NAME,
"/temp/IP_rep_", i, sep = "")
all_filepath[i] = PARAMETERS$IP_BAM[i]
print(paste("Processing IP sample", i))
temp <- exomepeak(GENE_ANNO_GTF = GENE_ANNO_GTF,
GENOME = GENOME, UCSC_TABLE_NAME = UCSC_TABLE_NAME,
TXDB = TXDB, IP_BAM = c(PARAMETERS$IP_BAM[i]),
INPUT_BAM = c(PARAMETERS$INPUT_BAM[i]), OUTPUT_DIR = PARAMETERS$EXOME_OUTPUT_DIR,
POISSON_MEAN_RATIO = 1, EXPERIMENT_NAME = experiment_name)
}
}
else {
for (i in 1:ip_length) {
experiment_name = paste(PARAMETERS$GO_EXPERIMENT_NAME,
"/temp/IP_rep_", i, sep = "")
input_bam = PARAMETERS$INPUT_BAM[PARAMETERS$INPUT2IP[[i]][1]]
if (length(PARAMETERS$INPUT2IP[[i]]) > 1) {
for (j in 2:length(PARAMETERS$INPUT2IP[[i]])) {
input_bam = c(input_bam, c(PARAMETERS$INPUT_BAM[PARAMETERS$INPUT2IP[[i]][j]]))
}
}
all_filepath[i] <- paste(PARAMETERS$EXOME_OUTPUT_DIR,
experiment_name, sep = "/")
print(paste("Processing IP sample", i))
temp <- exomepeak(GENE_ANNO_GTF = GENE_ANNO_GTF,
GENOME = GENOME, UCSC_TABLE_NAME = UCSC_TABLE_NAME,
TXDB = TXDB, IP_BAM = c(PARAMETERS$IP_BAM[i]),
INPUT_BAM = input_bam, OUTPUT_DIR = PARAMETERS$EXOME_OUTPUT_DIR,
POISSON_MEAN_RATIO = 1, EXPERIMENT_NAME = experiment_name)
}
}
for (i in 1:length(all_filepath)) {
path <- paste(all_filepath[i], "peak.xls", sep = "/")
bam_file <- read.table(path, sep = "\t", header = FALSE)
if (i == 1) {
all_peak <- bam_file
}
else {
all_peak = rbind(all_peak, bam_file[2:nrow(bam_file),
])
}
}
dir.create(paste(PARAMETERS$GO_OUTPUT_DIR, PARAMETERS$GO_EXPERIMENT_NAME,
sep = "/"), recursive = TRUE, showWarnings = FALSE)
dir = paste(PARAMETERS$GO_OUTPUT_DIR, PARAMETERS$GO_EXPERIMENT_NAME,
sep = "/")
write.table(all_peak, paste(dir, "all_peak.xls", sep = "/"),
sep = "\t", quote = FALSE, col.names = FALSE, row.names = FALSE)
filepath = paste(dir, "all_peak.xls", sep = "/")
matrix_peak_read = read.table(filepath, header = TRUE, stringsAsFactors = FALSE)
ori_peak <- .xls2Grangeslist(filepath)
overlaps <- findOverlaps(ori_peak, ori_peak)
matrix_overlap <- matrix(0, length(overlaps), 2)
matrix_overlap[, 1] <- queryHits(overlaps)
matrix_overlap[, 2] <- subjectHits(overlaps)
num <- c(rep(0, length(ori_peak)))
for (i in 1:nrow(matrix_overlap)) {
if (matrix_overlap[i, 1][1][[1]] != matrix_overlap[i,
2][1][[1]]) {
x <- ranges(ori_peak[matrix_overlap[i, 1]])[[1]]@width
y <- ranges(ori_peak[matrix_overlap[i, 2]])[[1]]@width
if (x < y) {
num[matrix_overlap[i, 1][1][[1]]] <- matrix_overlap[i,
1][1][[1]]
}
else if (x > y) {
num[matrix_overlap[i, 2][1][[1]]] <- matrix_overlap[i,
2][1][[1]]
}
else {
num[matrix_overlap[i, 1][1][[1]]] <- matrix_overlap[i,
1][1][[1]]
num[matrix_overlap[i, 2][1][[1]]] <- matrix_overlap[i,
2][1][[1]]
}
}
}
peak <- ori_peak[num == 0]
tab <- which(num == 0)
write.table(matrix_peak_read[tab, ], paste(dir, "merged_peak.xls",
sep = "/"), sep = "\t", quote = FALSE, col.names = TRUE,
row.names = FALSE)
write.table(matrix_peak_read[tab, 1:12], paste(dir, "merged_peak.bed",
sep = "/"), sep = "\t", quote = FALSE, col.names = FALSE,
row.names = FALSE)
txdb <- .readTxDb2(PARAMETERS)
exonRanges <- exonsBy(txdb, "gene")
rpkm_input <- matrix(0, length(peak), length(PARAMETERS$INPUT_BAM))
rpkm_ip <- matrix(0, length(peak), length(PARAMETERS$IP_BAM))
count_input <- matrix(0, length(peak), length(PARAMETERS$INPUT_BAM))
count_ip <- matrix(0, length(peak), length(PARAMETERS$IP_BAM))
for (i in 1:length(PARAMETERS$INPUT_BAM)) {
filepath <- PARAMETERS$INPUT_BAM[i]
aligns <- readGAlignments(filepath)
para <- ScanBamParam(what = "mapq")
mapq <- scanBam(filepath, param = para)[[1]][[1]]
mapq[is.na(mapq)] <- 255
ID_keep <- (mapq > 30)
filtered <- aligns[ID_keep]
id <- countOverlaps(filtered, exonRanges)
transcriptome_filtered_aligns <- filtered[id > 0]
counts <- countOverlaps(peak, transcriptome_filtered_aligns)
count_input[, i] <- counts
numBases <- sum(width(peak))
geneLengthsInKB <- numBases/1000
millionsMapped <- length(transcriptome_filtered_aligns)/1e+06
rpm <- counts/millionsMapped
rpkm_input[, i] <- rpm/geneLengthsInKB
}
write.table(count_input, paste(dir, "count_input.xls", sep = "/"),
sep = "\t", quote = FALSE, col.names = TRUE, row.names = FALSE)
write.table(rpkm_input, paste(dir, "rpkm_input.xls", sep = "/"),
sep = "\t", quote = FALSE, col.names = TRUE, row.names = FALSE)
for (i in 1:length(PARAMETERS$IP_BAM)) {
filepath <- PARAMETERS$IP_BAM[i]
aligns <- readGAlignments(filepath)
para <- ScanBamParam(what = "mapq")
mapq <- scanBam(filepath, param = para)[[1]][[1]]
mapq[is.na(mapq)] <- 255
ID_keep <- (mapq > 30)
filtered <- aligns[ID_keep]
id <- countOverlaps(filtered, exonRanges)
transcriptome_filtered_aligns <- filtered[id > 0]
counts <- countOverlaps(peak, transcriptome_filtered_aligns)
count_ip[, i] <- counts
numBases <- sum(width(peak))
geneLengthsInKB <- numBases/1000
millionsMapped <- length(transcriptome_filtered_aligns)/1e+06
rpm <- counts/millionsMapped
rpkm_ip[, i] <- rpm/geneLengthsInKB
}
write.table(count_ip, paste(dir, "count_ip.xls", sep = "/"),
sep = "\t", quote = FALSE, col.names = TRUE, row.names = FALSE)
write.table(rpkm_ip, paste(dir, "rpkm_ip.xls", sep = "/"),
sep = "\t", quote = FALSE, col.names = TRUE, row.names = FALSE)
matrix_log2_fc = matrix(0, length(peak), ip_length)
if (is.na(PARAMETERS$INPUT2IP)) {
for (i in 1:ip_length) {
matrix_log2_fc[, i] <- log2(rpkm_ip[, i] + 0.01) -
log2(rpkm_input[, i] + 0.01)
}
}
else {
for (i in 1:ip_length) {
if (length(PARAMETERS$INPUT2IP[[i]]) == 1) {
matrix_log2_fc[, i] <- log2(rpkm_ip[, i] + 0.01) -
log2(rpkm_input[, PARAMETERS$INPUT2IP[[i]]] +
0.01)
}
else {
sum = matrix(0, length(peak), 1)
for (j in 1:length(PARAMETERS$INPUT2IP[i])) {
sum = sum + rpkm_input[, PARAMETERS$INPUT2IP[[i]][j]]
}
matrix_log2_fc[, i] <- log2(rpkm_ip[, i] + 0.01) -
log2(sum + 0.01)
}
}
}
v <- matrix_log2_fc[, 1]
for (i in 1:nrow(matrix_log2_fc)) {
v[i] <- var(matrix_log2_fc[i, ])
}
matrix_log2_fc <- cbind(tab, matrix_log2_fc[1:nrow(matrix_log2_fc),
], matrix_peak_read[tab, ])
write.table(matrix_log2_fc, paste(dir, "all_information.xls",
sep = "/"), sep = "\t", quote = FALSE, col.names = TRUE,
row.names = FALSE)
file.remove(paste(dir, "all_peak.xls", sep = "/"))
file.remove(paste(dir, "all_information.xls", sep = "/"))
print("***************************")
print("***************************")
print(paste("The combinatorial RNA methylome is generated under:",
dir))
print("Including:")
print("1. Merged Peaks")
print("2. Merged Peaks in BED format")
print("3. Reads count for every peak in every bam file")
print("4. RPKM for every peak in every bam file")
}
--- function search by body ---
Function RMT in namespace exomePeak has this body.
----------- END OF FAILURE REPORT --------------
Fatal error: the condition has length > 1
* checking for unstated dependencies in vignettes ... OK
* checking package vignettes in 'inst/doc' ... OK
* checking running R code from vignettes ... SKIPPED
* checking re-building of vignette outputs ... SKIPPED
* checking PDF version of manual ... OK
* DONE
Status: 2 ERRORs, 3 NOTEs
See
'C:/Users/biocbuild/bbs-3.9-bioc/meat/exomePeak.Rcheck/00check.log'
for details.