assembly and annotation bed file

[Lei Jiang]

Hi, I am using methykit to analyze the DNA methylation difference in a marine invertebrate with a reference genome, however, I do not know how to set up the assembly and bed file in codes. Also the bed file is not in the ucsc website and how I can get this bed file. Very new in this field and my question may sound stupid.

[Alexander Blume]

Hi,

The assembly variable for methylKit is just there for bookkeeping,  we are not using it internally, so you do not have to do anything special here. 

Concerning the gene annotation bed file, you can perform all the steps of the normal comparative analysis without the need for the annotation file. 

If UCSC supports your organism, the BED file can be obtained from the table browser (https://genome-euro.ucsc.edu/cgi-bin/hgTables).

If you have a gtf/gff instead of a bed file there is some untested code to use these instead of the BED file (https://github.com/BIMSBbioinfo/genomation/issues/192).

Best,

Alex

[Georges Kanaan]

I have a similar issue. I have a GFF file produced by prodigal for my bacterial genome, and using the code in issue 192 fails for me.

```
gene.obj = gffToGRanges("contigs.emapper.genepred.gff")                      
head(gene.obj)                                                                                                                                        
split = as(split(gene.obj, gene.obj$type), "GRangesList")                                                                                                
head(split)                                                                                                                                              
                                                                                                                                                           
diffAnn=annotateWithGeneParts(as(myDiffFiltered,"GRanges"), split)
```

I get:
```
GRanges object with 6 ranges and 17 metadata columns:
          seqnames    ranges strand |          source     type     score
             <Rle> <IRanges>  <Rle> |        <factor> <factor> <numeric>
  [1] contig_17563    1-1596      + | Prodigal_v2.6.3      CDS     330.7
  [2] contig_17563 1598-1891      + | Prodigal_v2.6.3      CDS      30.2
  [3] contig_17563 1888-2376      + | Prodigal_v2.6.3      CDS      60.8
  [4] contig_17563 2426-3826      + | Prodigal_v2.6.3      CDS     288.1
  [5] contig_17563 3898-4434      + | Prodigal_v2.6.3      CDS      68.8
  [6] contig_17563 4457-5341      - | Prodigal_v2.6.3      CDS      68.3
          phase          ID     partial  start_type   rbs_motif  rbs_spacer
      <integer> <character> <character> <character> <character> <character>
  [1]         0         1_1          10        Edge        None        None
  [2]         0         1_2          00         ATG GGA/GAG/AGG      5-10bp
  [3]         0         1_3          00         TTG        None        None
  [4]         0         1_4          00         ATG        None        None
  [5]         0         1_5          00         ATG        None        None
  [6]         0         1_6          00         ATG        None        None
          gc_cont        conf     score.1      cscore      sscore      rscore
      <character> <character> <character> <character> <character> <character>
  [1]       0.336       99.99      330.69      327.47        3.22        0.00
  [2]       0.320       99.90       30.23       27.48        2.75        0.05
  [3]       0.342      100.00       60.77       67.19       -6.42       -0.32
  [4]       0.297       99.99      288.13      282.39        5.74       -0.32
  [5]       0.302      100.00       68.78       66.14        2.64       -0.32
  [6]       0.288      100.00       68.35       65.75        2.59       -0.32
           uscore      tscore
      <character> <character>
  [1]        0.00        3.22
  [2]       -1.08        2.53
  [3]       -3.25       -7.72
  [4]        2.79        2.53
  [5]        0.43        2.53
  [6]        1.04        2.53
  -------
  seqinfo: 30 sequences from an unspecified genome; no seqlengths
GRangesList object of length 1:
$CDS
GRanges object with 1640 ranges and 17 metadata columns:
             seqnames    ranges strand |          source     type     score
                <Rle> <IRanges>  <Rle> |        <factor> <factor> <numeric>
     [1] contig_17563    1-1596      + | Prodigal_v2.6.3      CDS     330.7
     [2] contig_17563 1598-1891      + | Prodigal_v2.6.3      CDS      30.2
     [3] contig_17563 1888-2376      + | Prodigal_v2.6.3      CDS      60.8
     [4] contig_17563 2426-3826      + | Prodigal_v2.6.3      CDS     288.1
     [5] contig_17563 3898-4434      + | Prodigal_v2.6.3      CDS      68.8
     ...          ...       ...    ... .             ...      ...       ...
  [1636] contig_12307  129-1148      - | Prodigal_v2.6.3      CDS     274.0
  [1637] contig_12307 1167-1610      - | Prodigal_v2.6.3      CDS      81.5
  [1638] contig_12307 1718-2371      - | Prodigal_v2.6.3      CDS     120.9
  [1639] contig_12307 2364-2537      - | Prodigal_v2.6.3      CDS       4.0
  [1640] contig_12307 2515-3009      - | Prodigal_v2.6.3      CDS     104.2
             phase          ID     partial  start_type   rbs_motif  rbs_spacer
         <integer> <character> <character> <character> <character> <character>
     [1]         0         1_1          10        Edge        None        None
     [2]         0         1_2          00         ATG GGA/GAG/AGG      5-10bp
     [3]         0         1_3          00         TTG        None        None
     [4]         0         1_4          00         ATG        None        None
     [5]         0         1_5          00         ATG        None        None
     ...       ...         ...         ...         ...         ...         ...
  [1636]         0        30_1          00         ATG       AGGAG      5-10bp
  [1637]         0        30_2          00         GTG        None        None
  [1638]         0        30_3          00         GTG        None        None
  [1639]         0        30_4          00         TTG        None        None
  [1640]         0        30_5          01        Edge        None        None
             gc_cont        conf     score.1      cscore      sscore
         <character> <character> <character> <character> <character>
     [1]       0.336       99.99      330.69      327.47        3.22
     [2]       0.320       99.90       30.23       27.48        2.75
     [3]       0.342      100.00       60.77       67.19       -6.42
     [4]       0.297       99.99      288.13      282.39        5.74
     [5]       0.302      100.00       68.78       66.14        2.64
     ...         ...         ...         ...         ...         ...
  [1636]       0.315       99.99      274.00      259.67       14.34
  [1637]       0.279      100.00       81.49       88.77       -7.28
  [1638]       0.317      100.00      120.92      124.82       -3.90
  [1639]       0.230       71.31        3.96       16.36      -12.40
  [1640]       0.277      100.00      104.21      100.99        3.22
              rscore      uscore      tscore
         <character> <character> <character>
     [1]        0.00        0.00        3.22
     [2]        0.05       -1.08        2.53
     [3]       -0.32       -3.25       -7.72
     [4]       -0.32        2.79        2.53
     [5]       -0.32        0.43        2.53
     ...         ...         ...         ...
  [1636]        7.97        2.76        2.53
  [1637]       -0.32       -0.55       -6.41
  [1638]       -0.32        1.57       -6.41
  [1639]       -0.46       -0.65      -11.28
  [1640]        0.00        0.00        3.22
  -------
  seqinfo: 30 sequences from an unspecified genome; no seqlengths

Error in (function (classes, fdef, mtable)  :
  impossible to find a méthode inheriting for the fonction ‘countOverlaps’ for the signature ‘"GRanges", "NULL"’
Call : annotateWithGeneParts ... [<- -> [<-.data.frame -> countOverlaps -> <Anonymous>
```
I would appreciate any tips that made this work for others.

[Alexander Blume]

Hi,

for any gff-based annotation prepared with gffToGRanges() please use `annotateWithFeatures()` instead of `annotateWithGeneParts()` for the final step.

annotateWithGeneParts() expects as features (second arg) a "GRangesList object containing GRanges object for promoter, exons, introns and transcription start sites, or simply output of readTranscriptFeatures function",
while annotateWithFeatures() just needs "a named GRangesList object containing GRanges objects different set of features".

To fix your code: 

```
gene.obj = gffToGRanges("contigs.emapper.genepred.gff")                      
head(gene.obj)                                                                                                                                        
split = as(split(gene.obj, gene.obj$type), "GRangesList")                                                                                                
head(split)                                                                                                                                              
                                                                                                                                                           

diffAnn=annotateWithFeatures(as(myDiffFiltered,"GRanges"), split)
```

Best,
Alex