Issue with annotateWithFeatures #214
Description
I am running into an error when running annotateWithFeatures, I have attached the code run and the error. I have verified that the gff file I use matches the file used for alignment. I believe the issue however still might be a mismatch between the methylation object and the gff, so I have attached the head() output here as well. I have also run the code using a bed file with annotateWithGeneParts however there I received a different error. Does anyone happen to know what could be causing this issue?
CODE:
Reading methylation data
`library(methylKit)
library(genomation)
library(GenomicRanges)
library(rtracklayer) # for reading GFF files
file.list <- list(
"D:/bismark/SF-1F_Liver_S67_R1_001_val_1_bismark_bt2_pe.deduplicated.bismark.cov",
"D:/bismark/SF-2F_Liver_S73_R1_001_val_1_bismark_bt2_pe.deduplicated.bismark.cov",
"D:/bismark/SF-3F_Liver_S79_R1_001_val_1_bismark_bt2_pe.deduplicated.bismark.cov",
"D:/bismark/SF-4F_Liver_S85_R1_001_val_1_bismark_bt2_pe.deduplicated.bismark.cov",
"D:/bismark/SF-5F_Liver_S90_R1_001_val_1_bismark_bt2_pe.deduplicated.bismark.cov",
"D:/bismark/SF-6F_Liver_S95_R1_001_val_1_bismark_bt2_pe.deduplicated.bismark.cov"
)
myobj <- methRead(file.list,
sample.id=list("SF_Liver_1", "SF_Liver_2", "SF_Liver_3",
"SF_Liver_4", "SF_Liver_5", "SF_Liver_6"),
pipeline="bismarkCoverage",
assembly="GCF_023375975.1", # genome assembly number for Astyanax
treatment=c(1,1,1,0,0,0),
mincov=1
)
myobj.filt <- filterByCoverage(myobj, lo.count=1, lo.perc=NULL, hi.count=NULL, hi.perc=99.9)
meth <- unite(myobj.filt, destrand=FALSE)
read gff file and prepare GRanges
na.omit(meth)
#bed = readTranscriptFeatures("C:/Users/marcf/Downloads/BedTest.bed/BedTest.bed",remove.unusual=FALSE)
#head(bed)
gff = gffToGRanges("D:/bismark/GCF_023375975.1_AstMex3_surface_genomic.gff")
head(gff)
split = as(split(gff, gff$type), "GRangesList")
methGrange <- as(meth, "GRanges")
head(split)
head(methGrange)
annotateWithFeatures and processing
annotated_meth <- annotateWithFeatures(methGrange, split)
dm_genes <- as.data.frame(annotated_meth$gene.obj)
write.csv(dm_genes, file="methylated_genes_in_liver.csv", row.names=FALSE)
plotTargetAnnotation(annotated_meth, main="Methylation Annotation to Gene Regions",
precedence=TRUE)
myDiff <- calculateDiffMeth(meth, covariates=NULL)
dmr_annotated <- annotateWithGeneParts(as(myDiff, "GRanges"), split)
dmr_genes <- as.data.frame(dmr_annotated$gene.obj)
write.csv(dmr_genes, file="dmr_genes_in_liver.csv", row.names=FALSE)`
ERROR:
Error in h(simpleError(msg, call)) : error in evaluating the argument 'x' in selecting a method for function 'as.matrix': Item 1 of j is 2 which is outside the column number range [1,ncol=1] In addition: Warning message: In .merge_two_Seqinfo_objects(x, y) : The 2 combined objects have no sequence levels in common. (Use suppressWarnings() to suppress this warning.)
HEAD GFF:
seqnames ranges strand | source type score phase ID Dbxref Name chromosome
|
[1] NC_064408.1 1-134019835 + | RefSeq region NA NC_064408.1:1..13401.. taxon:7994 1 1
[2] NC_064409.1 1-78662058 + | RefSeq region NA NC_064409.1:1..78662.. taxon:7994 2 2
[3] NC_064410.1 1-65479127 + | RefSeq region NA NC_064410.1:1..65479.. taxon:7994 3 3
[4] NC_064411.1 1-58968335 + | RefSeq region NA NC_064411.1:1..58968.. taxon:7994 4 4
[5] NC_064412.1 1-58297222 + | RefSeq region NA NC_064412.1:1..58297.. taxon:7994 5 5
... ... ... ... . ... ... ... ... ... ... ... ...
[105] NW_026040105.1 1-40726 + | RefSeq region NA NW_026040105.1:1..40.. taxon:7994 Unknown Unknown
[106] NW_026040106.1 1-40445 + | RefSeq region NA NW_026040106.1:1..40.. taxon:7994 Unknown Unknown
[107] NW_026040107.1 1-38578 + | RefSeq region NA NW_026040107.1:1..38.. taxon:7994 Unknown Unknown
[108] NW_026040108.1 1-38391 + | RefSeq region NA NW_026040108.1:1..38.. taxon:7994 Unknown Unknown
[109] NW_026040109.1 1-36882 + | RefSeq region NA NW_026040109.1:1..36.. taxon:7994 Unknown Unknown
dev-stage gbkey genome isolate mol_type sex tissue-type description gene gene_biotype Parent
HEAD METHYLATION OBJECT:
seqnames ranges strand | score name
|
[1] NC_064408.1 45299-45933 + | 1 XM_022663405.2
[2] NC_064408.1 51644-51709 + | 2 XM_022663405.2
[3] NC_064408.1 52981-53133 + | 3 XM_022663405.2
[4] NC_064408.1 54101-54235 + | 4 XM_022663405.2
[5] NC_064408.1 55099-55275 + | 5 XM_022663405.2
... ... ... ... . ... ...
[800709] NW_026040109.1 267-385 + | 1 XM_049474205.1
[800710] NW_026040109.1 1839-1990 + | 2 XM_049474205.1
[800711] NW_026040109.1 6861-6973 + | 3 XM_049474205.1
[800712] NW_026040109.1 17441-17535 + | 4 XM_049474205.1
[800713] NW_026040109.1 17655-17899 + | 5 XM_049474205.1