assembly and annotation bed file

Lei Jiang

unread,

Apr 13, 2023, 12:08:44 AM4/13/23

to methylkit_discussion

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

unread,

Apr 19, 2023, 9:06:22 AM4/19/23

to methylkit_discussion

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

unread,

Sep 26, 2023, 3:36:07 PM9/26/23

to methylkit_discussion

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

unread,

Oct 9, 2023, 8:41:27 AM10/9/23

to methylkit_discussion

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

Adeline M

unread,

Nov 4, 2025, 10:38:34 AM (5 days ago) Nov 4

to methylkit_discussion

Hello,

I'm in the same situation with a GTF file and I am having trouble understanding exactly what the annotateWithFeatures function does. When I use this function, 60% of my DMCs match the "gene" feature, but when I manually match each cytosine with the gene identifiers, only about 40% match LOC identifiers, meaning the others are in intergenic regions. Could you give me more details about this function?

Thanks and have a good day,

Adeline

Alexander Blume

unread,

8:08 AM (22 minutes ago) 8:08 AM

to methylkit_...@googlegroups.com

Hi Adeline,

This function is part of the genomation package, and the code can be found here: https://github.com/BIMSBbioinfo/genomation/blob/fa285d07802273fe3cf1a0b2a7b4de034a34544d/R/readAnnotate.R#L871
You could try to compare the output of your code with our implementation and see if there are any differences. I would be happy to hear about any issues you might find or suggestions on how we can improve it.

Best,
Alex

--
You received this message because you are subscribed to the Google Groups "methylkit_discussion" group.
To unsubscribe from this group and stop receiving emails from it, send an email to methylkit_discus...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/methylkit_discussion/81b4a724-32cd-40e4-b193-8882b92063dan%40googlegroups.com.

Reply all

Reply to author

Forward