How to produce the fraction of reads in peaks (FRiP) using DiffBind
1
0
Entering edit mode
Gary ▴ 20
@gary-7967
Last seen 2.6 years ago

Hi,

We have three mouse ATAC-Seq data. However, there is an error message when I would like to estimate their FRiP values using DiffBind below (i.e., Error in if (colnames(res)[i] == "Peak caller") { : argument is of length zero). Could you help me? Many thanks.

Best,

Gary

> sessionInfo()
R version 3.5.0 (2018-04-23)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Sierra 10.12.6

Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] parallel  stats4    stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] DiffBind_2.10.0             SummarizedExperiment_1.12.0 DelayedArray_0.8.0          BiocParallel_1.16.6         matrixStats_0.54.0         
 [6] Biobase_2.42.0              GenomicRanges_1.34.0        GenomeInfoDb_1.18.2         IRanges_2.16.0              S4Vectors_0.20.1           
[11] BiocGenerics_0.28.0        

loaded via a namespace (and not attached):
 [1] Category_2.48.1          bitops_1.0-6             bit64_0.9-7              RColorBrewer_1.1-2       progress_1.2.0          
 [6] httr_1.4.0               Rgraphviz_2.26.0         backports_1.1.3          tools_3.5.0              R6_2.4.0                
[11] KernSmooth_2.23-15       DBI_1.0.0                lazyeval_0.2.1           colorspace_1.4-0         tidyselect_0.2.5        
[16] prettyunits_1.0.2        bit_1.1-14               compiler_3.5.0           sendmailR_1.2-1          graph_1.60.0            
[21] rtracklayer_1.42.2       checkmate_1.9.1          caTools_1.17.1.2         scales_1.0.0             BatchJobs_1.7           
[26] genefilter_1.64.0        RBGL_1.58.1              stringr_1.4.0            digest_0.6.18            Rsamtools_1.34.1        
[31] AnnotationForge_1.24.0   XVector_0.22.0           base64enc_0.1-3          pkgconfig_2.0.2          limma_3.38.3            
[36] rlang_0.3.1              rstudioapi_0.9.0         RSQLite_2.1.1            BBmisc_1.11              GOstats_2.48.0          
[41] hwriter_1.3.2            gtools_3.8.1             dplyr_0.8.0.1            RCurl_1.95-4.12          magrittr_1.5            
[46] GO.db_3.7.0              GenomeInfoDbData_1.2.0   Matrix_1.2-16            Rcpp_1.0.0               munsell_0.5.0           
[51] stringi_1.3.1            edgeR_3.24.3             zlibbioc_1.28.0          gplots_3.0.1.1           plyr_1.8.4              
[56] grid_3.5.0               blob_1.1.1               ggrepel_0.8.0            gdata_2.18.0             crayon_1.3.4            
[61] lattice_0.20-38          Biostrings_2.50.2        splines_3.5.0            GenomicFeatures_1.34.4   annotate_1.60.1         
[66] hms_0.4.2                locfit_1.5-9.1           pillar_1.3.1             rjson_0.2.20             systemPipeR_1.16.1      
[71] biomaRt_2.38.0           XML_3.98-1.19            glue_1.3.0               ShortRead_1.40.0         latticeExtra_0.6-28     
[76] data.table_1.12.0        gtable_0.2.0             purrr_0.3.1              amap_0.8-16              assertthat_0.2.0        
[81] ggplot2_3.1.0            xtable_1.8-3             survival_2.43-3          pheatmap_1.0.12          tibble_2.0.1            
[86] GenomicAlignments_1.18.1 AnnotationDbi_1.44.0     memoise_1.1.0            brew_1.0-6               GSEABase_1.44.0         
> atac <- dba(sampleSheet = "atac.csv", bRemoveM = TRUE, bRemoveRandom = FALSE, bCorPlot = TRUE, peakFormat = "narrow")
G12_1 Dermis Central Replicate WIHN 1 narrow
G12_2 Dermis Central Replicate WIHN 2 narrow
G12_Combine Dermis Central Combine WIHN 3 narrow
> atac
3 Samples, 999 sites in matrix (2690 total):
           ID Tissue  Factor Condition Treatment Replicate Caller Intervals
1       G12_1 Dermis Central Replicate      WIHN         1 narrow       416
2       G12_2 Dermis Central Replicate      WIHN         2 narrow       881
3 G12_Combine Dermis Central   Combine      WIHN         3 narrow      2597
> atac <- dba.count(atac, minOverlap = 1)
> atac
3 Samples, 2690 sites in matrix:
           ID Tissue  Factor Condition Treatment Replicate Caller Intervals
1       G12_1 Dermis Central Replicate      WIHN         1 counts      2690
2       G12_2 Dermis Central Replicate      WIHN         2 counts      2690
3 G12_Combine Dermis Central   Combine      WIHN         3 counts      2690
> dba.show(atac, atac$masks$Dermis, bContrasts=FALSE, DBA_FRIP)
Error in if (colnames(res)[i] == "Peak caller") { : 
  argument is of length zero
DiffBind FRiP ATAC-Seq dba.show dba.count • 609 views
ADD COMMENT
0
Entering edit mode
Rory Stark ★ 4.1k
@rory-stark-5741
Last seen 19 hours ago
CRUK, Cambridge, UK

Hmmn, for some reason the FRiP values are not being computed by dba.count(), or else they would have shown up then you listed the DBA object atac after counting. If you can send me copies of the DBA object atac before and after the call to dba,count() (or a link where I can download them), I may be able to see what is going on.

-Rory

ADD COMMENT
0
Entering edit mode
 > atac <- dba(sampleSheet = "atac.csv", bRemoveM = TRUE, bRemoveRandom = FALSE, 
               bCorPlot = TRUE, peakFormat = "narrow")
 > dba.save(atac,"atac_peaks")
 > atac <- dba.count(atac, minOverlap = 1)
 > dba.save(atac,"atac_counts")

which will create files called dba_ata_peaks.RData and dba_ata_counts.RData

ADD REPLY
0
Entering edit mode

Hi Rory,

Many thanks. Here you are.

http://68.181.92.180/~Gary/temporal/20190325DiffBindATACseqFRiP/dbaatac_peaks.RData

http://68.181.92.180/~Gary/temporal/20190325DiffBindATACseqFRiP/dbaatac_counts.RData

By the way, I find that DiffBind can estimate FRiP values based on bam files directly produced by Bowtie2 on Partek Flow. However, after trimming duplicate reads (Picard), mitochondrial reads, low mapping quality reads, and improper paired-end reads (Bamtools) on public Galaxy, DiffBind cannot estimate FRiP values based on these trimmed bam files.

Best,

Gary

ADD REPLY

Login before adding your answer.

Traffic: 297 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6