Error in plotDI in pacakge diffHic
2
0
Entering edit mode
Andrew • 0
@andrew-11520
Last seen 5.0 years ago

Hi Aaron,

I tried to use plotDI, but it occurs error as following:

> plotDI(iset, test.fc,  first.region=GRanges("chr1", IRanges(22240001, 22280000)), second.region=GRanges("chr1", IRanges(3000001, 3040000)), diag=FALSE)
Error in if (first.min > first.max || second.min > second.max) { :
  missing value where TRUE/FALSE needed
 

Would you help me how to figure it out? Thanks.
 

diffHic plotDI • 566 views
ADD COMMENT
0
Entering edit mode
Aaron Lun ★ 27k
@alun
Last seen 1 hour ago
The city by the bay

I would guess that you haven't constructed iset properly, i.e., it doesn't contain sequence length information for chromosome 1. If you tried something like seqlengths(regions(iset))[["chr1"]], do you get a sensible value? It's probably NA instead, which leads to the observed error above.

ADD COMMENT
0
Entering edit mode
Andrew • 0
@andrew-11520
Last seen 5.0 years ago

I used method you posted at Loading a precomputed interaction matrix into diffHiC for iset file generation

------------------------------------------

> library(rtracklayer)
> bin.coords <- import("bin_4k.bed")
> bin.coords <- as(bin.coords, "GRanges")
> head(bin.coords)
GRanges object with 6 ranges and 0 metadata columns:
      seqnames           ranges strand
         <Rle>        <IRanges>  <Rle>
  [1]     chr1 [     1,  40000]      *
  [2]     chr1 [ 40001,  80000]      *
  [3]     chr1 [ 80001, 120000]      *
  [4]     chr1 [120001, 160000]      *
  [5]     chr1 [160001, 200000]      *
  [6]     chr1 [200001, 240000]      *
  -------
  seqinfo: 22 sequences from an unspecified genome; no seqlengths
> library(InteractionSet)
> a_40000.matrix <-read.table("a_40000.matrix")
> gi <- GInteractions(a_40000.matrix[,1], a_40000.matrix[,2], bin.coords, mode="reverse")
> head(gi)
ReverseStrictGInteractions object with 6 interactions and 0 metadata columns:
      seqnames1            ranges1     seqnames2            ranges2
          <Rle>          <IRanges>         <Rle>          <IRanges>
  [1]      chr1 [3000001, 3040000] ---      chr1 [3000001, 3040000]
  [2]      chr1 [3040001, 3080000] ---      chr1 [3000001, 3040000]
  [3]      chr1 [3080001, 3120000] ---      chr1 [3000001, 3040000]
  [4]      chr1 [3120001, 3160000] ---      chr1 [3000001, 3040000]
  [5]      chr1 [3160001, 3200000] ---      chr1 [3000001, 3040000]
  [6]      chr1 [3200001, 3240000] ---      chr1 [3000001, 3040000]
  -------
  regions: 68149 ranges and 0 metadata columns
  seqinfo: 22 sequences from an unspecified genome; no seqlengths

> iset<-InteractionSet(as.matrix(a_40000.matrix[,3]), gi)
> seqlengths(regions(iset))[["chr1"]]
[1] NA
> seqlengths(regions(iset))
 chr1  chr2  chr3  chr4  chr5  chr6  chr7  chr8  chr9 chr10 chr11 chr12 chr13
   NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
chr14 chr15 chr16 chr17 chr18 chr19  chrX  chrY  chrM
   NA    NA    NA    NA    NA    NA    NA    NA    NA
----------------------------------------------------------------

Like what you guessed, seqlengths(regions(iset))[["chr1"]]. how can I fix it? Thanks.

 

ADD COMMENT
0
Entering edit mode
  1. Post replies to answers as comments, rather than as new answers, unless you're answering your original post.
  2. Try using the seqlengths assignment function, e.g., seqlengths(bin.coords) <- SOME_SEQLEN_HERE.
ADD REPLY

Login before adding your answer.

Traffic: 311 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