gmat

[hjack123]

Hi, When creating gmat, is it possible to just make counts under a promoter region of a gene, instead of everywhere on the gene?

Thanks!!!

[ravipatel4]

@hjack123 This solution is not exactly what you are looking for, but it comes pretty close, I think. You can use coordinates of promoter regions (e.g. 1kb upstream and 0.5kb downstream of TSS, or your favorite definition of the promoter) instead of the coordinates of genes in the following code: genes = read.table("gencode.vM16.gene.bed"). Such that when you run createGmatFromMat, it will include counts from only those bins/peaks (depending on if you utilize "bmat" or "pmat") that overlap with the promoter regions.

[hjack123]

@ravipatel4 thanks I did the same at the end

[JihedC]

Hi @ravipatel4 ! Could you elaborate a bit more? I would like to count how many reads I get into repetitive regions upstream certains genes. I have imported my own bed file and created the IRange object.

x.sp = createGmatFromMat(
    obj=x.sp,
    input.mat="bmat",
    genes=DAR.select,
    do.par=TRUE,
    num.cores=2
  );

I get this error when I am trying to create the matrix for these intervals:

Error in createGmatFromMat.default(obj = x.sp, input.mat = "bmat", genes = DAR.gr, : genes does not contain gene names

My DAR.select Grange object looks like this :

GRanges object with 4 ranges and 1 metadata column:
        seqnames              ranges strand |       names
           <Rle>           <IRanges>  <Rle> | <character>
  peak2     chr1 143849553-143849709      * |       peak2
  peak3     chr1 105026965-105027121      * |       peak3
  peak4     chr1     4611985-4612141      * |       peak4
  peak5     chr1   51385207-51385363      * |       peak5
  -------
  seqinfo: 20 sequences from an unspecified genome; no seqlengths

while the one created form the read.table("gencode.vM16.gene.bed") looks like this:

GRanges object with 4 ranges and 1 metadata column:
      seqnames              ranges strand |     name
         <Rle>           <IRanges>  <Rle> | <factor>
  [1]     chr3   34650405-34652461      * |     Sox2
  [2]     chr4   55527143-55532466      * |     Klf4
  [3]     chr6 122707489-122714633      * |    Nanog
  [4]    chr17   35506018-35510772      * |   Pou5f1
  -------
  seqinfo: 22 sequences from an unspecified genome; no seqlengths

So I basically have no gene names to create the matrix but I have the range and a column called name.

If you could help me figure out what is wrong, this would be great!

Thanks in advance!

[ravipatel4]

@JihedC, Can you tell me how you generated your GRanges object?

It seems to work for me. I first read my bed intervals into peaksToGenes data frame. And then used the following code to generate a GRanges object.

# Using only first 6 bed intervals for this example.
genes.gr = GRanges(peaksToGenes[1:6,1], IRanges(peaksToGenes[1:6,2], peaksToGenes[1:6,3]), name= peaksToGenes[1:6,4] )

Here is how my GRanges look.

GRanges object with 6 ranges and 1 metadata column:
      seqnames          ranges strand |        name
         <Rle>       <IRanges>  <Rle> | <character>
  [1]     chr1 4492499-4492660      * |       peak1
  [2]     chr1 4492666-4492816      * |       peak2
  [3]     chr1 4496326-4497088      * |       peak3
  [4]     chr1 4525578-4525728      * |       peak4
  [5]     chr1 4571893-4572129      * |       peak5
  [6]     chr1 4600873-4601130      * |       peak6
  -------
  seqinfo: 22 sequences from an unspecified genome; no seqlengths

It works for me when I use the createGmatFromMat function using the GRanges created above.

More importantly though, another problem your approach has, apart from the error, is that you are using "bmat" matrix in createGmatFromMat, which, based on my understanding of SnapATAC, means that you won't get reads mapping to your bed intervals, but instead you will get reads mapping to the "bins" overlapping with your bed intervals. So, be careful how you use the "bmat".

If you are specifically interested in the number of peaks mapping to your bed intervals, you would first need to add a cell x repetitive regions count matrix to the snap object. For example,

snaptools snap-add-pmat \
        --snap-file my.snap \
        --peak-file repetitive.bed

where, repetitive.bed will have all regions that you are interested in and my.snap is the snap file you have been using so far. And then load the my.snap file into R and use "pmat" instead of "bmat" in createGmatFromMat.

DISCLAIMER: I haven't tested this approach; I cooked it up based on my understanding of how SnapATAC works. You might need to tweak it depending on the issues you may face. Perhaps that authors might be able to shed more light on this.

Good luck.