Entering edit mode
Hello there,
I have some problem when I was doing alternative splicing analysis with my RNAseq data via edgeR package.
I have six RNAseq samples(3 control, 3 treated) and I want to analyze the alternative splicing events by R. Here is the Syntax I used.
library(edgeR)
fc <- featureCounts(files=c("file control1.bam",
"file control2.bam",
"file control3.bam",
"file treated1.bam",
"file treated2.bam",
"file treated3.bam"),
annot.ext="file.gtf",
isGTFAnnotationFile=TRUE,
GTF.featureType = "exon",
GTF.attrType = "gene_id",
useMetaFeatures = FALSE,
allowMultiOverlap = TRUE,
isPairedEnd=TRUE)
fccounts <- fc$counts
metadata <- data.frame(sample_id = colnames(fccounts))
sample <- rep(c("con","treat"),each=3)
metadata$sample <- relevel(factor(sample),"con")
y.all <- DGEList(counts = fc$counts, genes = fc$annotation)
dim(y.all)
head(y.all$genes)
y <- y.all
y$samples
head(y$genes)
TREAT <- factor(metadata$sample, levels = c("con","treat"))
keep <- filterByExpr(y, group = TREAT)
table(keep)
y <- y[keep, keep.lib.sizes = FALSE]
y <- calcNormFactors(y)
y$samples
plotMDS(y)
Batch <- factor(c(1,2,3,1,2,3))
TREAT <- factor(metadata$sample, levels = c("con","treat"))
design <- model.matrix(~ Batch + TREAT)
design
y <- estimateDisp(y, design, robust = TRUE)
plotBCV(y)
fit <- glmQLFit(y, design, robust = TRUE)
plotQLDisp(fit)
qlf <- glmQLFTest(fit, coef = 4)
topTags(qlf)
is.de <- decideTests(qlf, p.value = 0.05)
summary(is.de)
sp <- diffSpliceDGE(fit, coef = 4, geneid = "GeneID", exonid = "Start")
But at this step the R returned like this:
Error in quantile.default(zresid, probs = prob) :
missing values and NaN's not allowed if 'na.rm' is FALSE
I have checked the NA values by sum(is.na(fit)) but no NA value was found. So I would like to ask how can I solve this "missing values and NaN's" problems?
Thank you very much for your kindly help!
Qiao
sessionInfo()
R version 4.1.0 (2021-05-18)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19041)
Matrix products: default
locale:
[1] LC_COLLATE=English_Netherlands.936 LC_CTYPE=English_Netherlands.936
[3] LC_MONETARY=English_Netherlands.936 LC_NUMERIC=C
[5] LC_TIME=English_Netherlands.936
attached base packages:
[1] parallel stats4 stats graphics grDevices utils datasets methods base
other attached packages:
[1] edgeR_3.32.1 DESeq2_1.30.1 SummarizedExperiment_1.20.0
[4] Biobase_2.50.0 MatrixGenerics_1.2.1 matrixStats_0.59.0
[7] GenomicRanges_1.42.0 GenomeInfoDb_1.26.7 IRanges_2.24.1
[10] S4Vectors_0.28.1 BiocGenerics_0.36.1 limma_3.46.0
[13] Rsubread_2.4.3
loaded via a namespace (and not attached):
[1] genefilter_1.72.1 statmod_1.4.36 locfit_1.5-9.4
[4] tidyselect_1.1.1 purrr_0.3.4 splines_4.1.0
[7] lattice_0.20-44 generics_0.1.0 colorspace_2.0-1
[10] vctrs_0.3.8 utf8_1.2.1 blob_1.2.1
[13] XML_3.99-0.6 survival_3.2-11 rlang_0.4.11
[16] pillar_1.6.1 glue_1.4.2 DBI_1.1.1
[19] BiocParallel_1.24.1 bit64_4.0.5 RColorBrewer_1.1-2
[22] GenomeInfoDbData_1.2.4 lifecycle_1.0.0 zlibbioc_1.36.0
[25] munsell_0.5.0 gtable_0.3.0 memoise_2.0.0
[28] geneplotter_1.68.0 fastmap_1.1.0 AnnotationDbi_1.52.0
[31] fansi_0.5.0 Rcpp_1.0.6 xtable_1.8-4
[34] scales_1.1.1 cachem_1.0.5 DelayedArray_0.16.3
[37] annotate_1.68.0 XVector_0.30.0 bit_4.0.4
[40] ggplot2_3.3.3 dplyr_1.0.6 grid_4.1.0
[43] tools_4.1.0 bitops_1.0-7 magrittr_2.0.1
[46] RCurl_1.98-1.3 RSQLite_2.2.7 tibble_3.1.2
[49] pkgconfig_2.0.3 crayon_1.4.1 ellipsis_0.3.2
[52] Matrix_1.3-4 httr_1.4.2 R6_2.5.0
[55] compiler_4.1.0
