Question: How to produce the fraction of reads in peaks (FRiP) using DiffBind
0
gravatar for Gary
29 days ago by
Gary10
Gary10 wrote:

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
ADD COMMENTlink modified 29 days ago by Rory Stark2.8k • written 29 days ago by Gary10
Answer: How to produce the fraction of reads in peaks (FRiP) using DiffBind
0
gravatar for Rory Stark
29 days ago by
Rory Stark2.8k
CRUK, Cambridge, UK
Rory Stark2.8k wrote:

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 COMMENTlink written 29 days ago by Rory Stark2.8k

Hi Rory, Thank you very much. Below are their bam, bai, narrowPeak, and CSV files. Is it OK? Could you teach me the command lines how to save the DBA object (atac) before and after the call to dba, count()? After that, I can provide the files you need. Thank you very much.

http://68.181.92.180/~Gary/temporal/20190325DiffBindATACseqFRiP/G121WithoutPartekTrimX2000GalaxyTrim.bam http://68.181.92.180/~Gary/temporal/20190325DiffBindATACseqFRiP/G122WithoutPartekTrimX2000GalaxyTrim.bam http://68.181.92.180/~Gary/temporal/20190325DiffBindATACseqFRiP/G12CombineWithoutPartekTrimX2000GalaxyTrim.bam http://68.181.92.180/~Gary/temporal/20190325DiffBindATACseqFRiP/G121WithoutPartekTrimX2000GalaxyTrim.bam.bai http://68.181.92.180/~Gary/temporal/20190325DiffBindATACseqFRiP/G122WithoutPartekTrimX2000GalaxyTrim.bai http://68.181.92.180/~Gary/temporal/20190325DiffBindATACseqFRiP/G12CombineWithoutPartekTrimX2000GalaxyTrim.bam.bai http://68.181.92.180/~Gary/temporal/20190325DiffBindATACseqFRiP/G121peaks.narrowPeak http://68.181.92.180/~Gary/temporal/20190325DiffBindATACseqFRiP/G122peaks.narrowPeak http://68.181.92.180/~Gary/temporal/20190325DiffBindATACseqFRiP/G12Combinepeaks.narrowPeak http://68.181.92.180/~Gary/temporal/20190325DiffBindATACseqFRiP/atac.csv

Best,

Gary

ADD REPLYlink modified 28 days ago • written 28 days ago by Gary10
 > 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 REPLYlink written 28 days ago by Rory Stark2.8k

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 REPLYlink written 27 days ago by Gary10
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 16.09
Traffic: 243 users visited in the last hour