CountOverlap
1
0
Entering edit mode
michela ▴ 20
@216f0d39
Last seen 15 months ago
Netherlands
# Loading a library
library(rtracklayer)

# assigning bed file location to variable
summits_file <- 'C:/Users/utente/Desktop/Tesi/Dati/peak_C1_H3K4me3_peaks.narrowPeak'

# import the BED data
twist.summits <- import(summits_file)

# what kind of object is this?
class(twist.summits)
twist.summits

> GRanges object with 226 ranges and 2 metadata columns:
        seqnames    ranges strand |                   name     score
           <Rle> <IRanges>  <Rle> |            <character> <numeric>
    [1]        I    801380      * | peak_C1_H3K4me3_peak_1   4.01076
    [2]        I    915512      * | peak_C1_H3K4me3_peak_2   3.18610
    [3]        I   2646543      * | peak_C1_H3K4me3_peak_3   3.94499
    [4]        I   2647293      * | peak_C1_H3K4me3_peak_4   4.32684
    [5]        I   2648461      * | peak_C1_H3K4me3_peak_5   8.31980

library(ChIPpeakAnno)
txdb <- TxDb.Celegans.UCSC.ce11.ensGene
annoData <- toGRanges(TxDb.Celegans.UCSC.ce11.ensGene, feature="gene")
annoData

### Get all the transcripts for each gene
allTx <- transcriptsBy(txdb, by='gene')

# examine the result
allTx
> GRangesList object of length 46904:
$WBGene00000001.1
GRanges object with 2 ranges and 2 metadata columns:
      seqnames          ranges strand |     tx_id      tx_name
         <Rle>       <IRanges>  <Rle> | <integer>  <character>
  [1]     chrI 5107833-5110183      + |      1083 Y110A7A.10.1
  [2]     chrI 5107843-5110164      + |      1084 Y110A7A.10.2
  -------
  seqinfo: 7 sequences (1 circular) from ce11 genome


BiocManager::install('biomaRt')
library(biomaRt)
ensamble<- biomaRt::useMart('ensembl', dataset = 'celegans_gene_ensembl') 
wormGenes <- biomaRt::getBM(attributes = c('external_gene_name','chromosome_name','start_position', 'end_position','gene_biotype','description'), mart= ensamble)


##################################################
## CODE BLOCK #3
##################################################

### Getting the foot print of all transcripts

# collapse all transcripts into a single range
genic.tx <- reduce(allTx)

# examine the result
genic.tx

# Then, let's get rid of the list format
genes <- unlist(genic.tx)

# examine the result
genes
> GRanges object with 46905 ranges and 0 metadata columns:
                   seqnames            ranges strand
                      <Rle>         <IRanges>  <Rle>
  WBGene00000001.1     chrI   5107833-5110183      +
  WBGene00000002.1    chrIV   9598986-9601695      -
  WBGene00000003.1     chrV   9244402-9246360      -


### Count how many peaks overlap each gene
overlap.genes <- suppressWarnings(countOverlaps(genes,twist.summits), classes = "warning")


head(overlap.genes)
> head(overlap.genes)
WBGene00000001.1 WBGene00000002.1 WBGene00000003.1 WBGene00000004.1 WBGene00000005.1 
               0                0                0                0                0 
WBGene00000006.1 
               0

This are my code but I don't understand why the overlap genes return zero

What I have to change?

GeneExpression ChIPSeqData • 906 views
ADD COMMENT
0
Entering edit mode
@james-w-macdonald-5106
Last seen 10 hours ago
United States

When you suppressWarnings you won't get warnings that tell you that there are problems! The main warning you will be suppressing will look something like this:

## These are not identical things!
> gr1 <- GRanges(rep("I", 5), IRanges(c(2,8,23,45,78), width = 5))
> gr2 <- GRanges(rep("chrI", 5), IRanges(c(2,8,23,45,78), width = 5))
> countOverlaps(gr1, gr2)
[1] 0 0 0 0 0
Warning message:
In .Seqinfo.mergexy(x, y) :
  The 2 combined objects have no sequence levels in common. (Use
  suppressWarnings() to suppress this warning.)

I mean, it does say you can suppress the warnings, but it's good to know why there isn't any overlap, no?

ADD COMMENT
0
Entering edit mode

overlap.genes <- countOverlaps(genes,twist.summits) Warning message: In .Seqinfo.mergexy(x, y) : The 2 combined objects have no sequence levels in common. (Use suppressWarnings() to suppress this warning.)

Yes but this is the only warning that I have, it don't say nothing more

ADD REPLY
0
Entering edit mode

That's true, but it tells you everything you need to know. Do you see the difference between the two GRanges objects I created?

ADD REPLY

Login before adding your answer.

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