Question: How to use DiffBind to get consensus peaks I am interested in?
0
gravatar for Gary
3.7 years ago by
Gary20
Gary20 wrote:

Hi,

I wound like to get 2,913 consensus peaks that in Early stage, but not in Late stage samples (fig1). They are overlapping between T8E_HMC & T9E_HMC samples, but not in T8L_HMC and T9L_HMC samples (fig1). I have tried the dba.peakset command, but doesn't work. Could you show me how to do it correctly? The best way is to get a .bed file for the downstream analysis. Below are commands I used to get this Venn diagram for your reference.  Thank you so much.

Best,

Gary

 

> t8t9=dba(sampleSheet = "t8t9.csv",peakFormat = "bed")

T8E_HMC MSC Early T8  1 MACS

T9E_HMC MSC Early T9  2 MACS

T8L_HMC MSC Late T8  1 MACS

T9L_HMC MSC Late T9  2 MACS

> t8t9

4 Samples, 15919 sites in matrix (39174 total):

       ID Tissue Factor Condition Replicate Caller Intervals

1 T8E_HMC    MSC  Early        T8         1   MACS     21905

2 T9E_HMC    MSC  Early        T9         2   MACS     23327

3 T8L_HMC    MSC   Late        T8         1   MACS     12377

4 T9L_HMC    MSC   Late        T9         2   MACS     12902

> names(t8t9$masks)

 [1] "MSC"         "Early"       "Late"        "T8"          "T9"          ""            "MACS"        "Replicate.1" "Replicate.2" "All"        

[11] "None"      

> dba.plotVenn(t8t9,t8t9$masks$MSC)

 

fig1

fig1

diffbind • 1.1k views
ADD COMMENTlink modified 3.7 years ago by Rory Stark2.9k • written 3.7 years ago by Gary20
Answer: How to use DiffBind to get consensus peaks I am interested in?
1
gravatar for Rory Stark
3.7 years ago by
Rory Stark2.9k
CRUK, Cambridge, UK
Rory Stark2.9k wrote:

Hi Gary-

You can use dba.overlap() as follows:

> t8t9.overlaps <- dba.overlap(t8t9,t8t9$masks$MSC)

t8t9.overlaps will contain the peaksets for all 15 positions in the Venn diagram:

> names(t8t9.overlaps)

The one you want will be called t8t9.overlaps$AandB:

> t8t9.overlaps$AandB

This is a GRanges object by default (you can get them back as a dataframe by setting DataType=DBA_DATA_FRAME).

You can also get the same data back from the call to dba.plotVenn() by setting bReturnPeaksets=TRUE (in the development version moving forward, these are returned "invisibly" and the bReturnPeaksets parameter has been eliminated).

Cheers-

Rory

ADD COMMENTlink written 3.7 years ago by Rory Stark2.9k

Hi Rory,

Thank you so much. I have tried what you suggested, but I don't know how to save the result as a txt file. I know it's not a DiffBind question, but could you help me again? Many thanks. Below are command lines I used for your reference.

Best,

Gary

 

> t8t9=dba(sampleSheet = "t8t9.csv",peakFormat = "bed")
T8E_HMC MSC Early T8  1 MACS
T9E_HMC MSC Early T9  2 MACS
T8L_HMC MSC Late T8  1 MACS
T9L_HMC MSC Late T9  2 MACS
> t8t9
4 Samples, 15919 sites in matrix (39174 total):
       ID Tissue Factor Condition Replicate Caller Intervals
1 T8E_HMC    MSC  Early        T8         1   MACS     21905
2 T9E_HMC    MSC  Early        T9         2   MACS     23327
3 T8L_HMC    MSC   Late        T8         1   MACS     12377
4 T9L_HMC    MSC   Late        T9         2   MACS     12902
> names(t8t9$masks)
 [1] "MSC"         "Early"       "Late"        "T8"          "T9"          ""           
 [7] "MACS"        "Replicate.1" "Replicate.2" "All"         "None"       
> dba.plotVenn(t8t9, t8t9$masks$MSC)
> t8t9.overlaps <-dba.overlap(t8t9,t8t9$masks$MSC)
> names(t8t9.overlaps)
 [1] "onlyA" "onlyB" "onlyC" "onlyD" "AandB" "AandC" "AandD" "BandC" "BandD" "CandD" "notA" 
[12] "notB"  "notC"  "notD"  "inAll"
> t8t9.overlaps$AandB
GRanges object with 2913 ranges and 2 metadata columns:
        seqnames               ranges strand   |             scoreA             scoreB
           <Rle>            <IRanges>  <Rle>   |          <numeric>          <numeric>
     24     chr1   [ 564352,  570445]      *   |  0.046839174582378 0.0234009360374415
     55     chr1   [ 900161,  901002]      *   | 0.0167048804454635 0.0227769110764431
     71     chr1   [1009781, 1010682]      *   | 0.0347199475925319 0.0277691107644306
     77     chr1   [1047306, 1048373]      *   | 0.0252210940058958 0.0162246489859594
     94     chr1   [1218174, 1218917]      *   | 0.0212905339010809 0.0287051482059282
    ...      ...                  ...    ... ...                ...                ...
  39020     chrY [ 2427812,  2428479]      *   | 0.0370127743203406 0.0271450858034321
  39091     chrY [10018797, 10020074]      *   | 0.0216180805764821 0.0168486739469579
  39125     chrY [13671130, 13672096]      *   | 0.0458565345561743 0.0215288611544462
  39127     chrY [13682255, 13683846]      *   | 0.0255486406812971 0.0290171606864275
  39137     chrY [13816319, 13823737]      *   | 0.0674746151326564 0.0464898595943838
  -------
  seqinfo: 33 sequences from an unspecified genome; no seqlengths
> write.table(t8t9.overlaps$AandB,"./t8t9.overlapsAandB.txt")
Error in as.vector(x) : no method for coercing this S4 class to a vector
> write.table(t8t9.overlaps$AandB,file="t8t9overlapsAandB.txt")
Error in as.vector(x) : no method for coercing this S4 class to a vector

ADD REPLYlink written 3.7 years ago by Gary20
2

Look here

A: diffbind: exporting unique peak sets

ADD REPLYlink written 3.7 years ago by Udi Landau30

Thanks a lot.

ADD REPLYlink written 3.7 years ago by Gary20

Thanks a lot.

ADD REPLYlink written 3.7 years ago by Gary20
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: 342 users visited in the last hour