使用rGADEM包GADEM函数时遇序列检索错误求助

问题详情 #

我有一份Bed格式的大峰列表，已转换为GenomicRange对象，使用rGADEM包寻找de novo基序时，运行GADEM函数出现如下错误：

[ Retrieving sequences... Error in.Call2("C_solve_user_SEW", refwidths, start, end, width, translate.negative.coord: solving row 136: 'allow.nonnarrowing' is FALSE and the supplied start (55134751) is > refwidth + 1 ]

补充说明：使用少于136行的示例输入文件时，程序可正常运行并得到基序。

运行代码：

library(GenomicRanges)
library(rGADEM)
data = makeGRangesFromDataFrame(data, keep.extra.columns = TRUE)
data = reduce(data)
data = resize(data, width = 50, fix='center')
gadem<-GADEM(data,verbose=1,genome=Hsapiens)
plot(gadem) 

20行示例数据：

1      chr6  29723590  29723790
2     chr14 103334312 103334512
3      chr1 150579030 150579230
4      chr7  76358527  76358727
5      chr6  11537891  11538091
6     chr14  49893256  49893456
7      chr5 179623200 179623400
8      chr1 228082831 228083031
9     chr12  93441644  93441844
10    chr10   3784776   3784976
11     chr3 183635833 183636033
12     chr7    975301    975501
13    chr12 123364510 123364710
14     chr1   1615578   1615778
15     chr1  36156320  36156520
16    chr14  55051781  55051981
17     chr8  11867697  11867897
18    chr22  38706135  38706335
19     chr6  44265256  44265456
20     chr1 185316658 185316858 

问题原因 #

错误指向第136行的区域：执行resize(data, width = 50, fix='center')时，该区域的中心位置过于靠近染色体末端，导致扩展后的起始位置超过了对应染色体的总长度，无法从参考基因组中提取有效序列。

解决方案 #

1. 过滤无效区域 #

在resize操作后，添加代码过滤掉起始位置小于1或终止位置超过染色体长度的区域：

library(GenomicRanges)
library(rGADEM)

data = makeGRangesFromDataFrame(data, keep.extra.columns = TRUE)
data = reduce(data)
data = resize(data, width = 50, fix='center')

# 获取参考染色体长度并过滤无效区域
chr_lengths <- seqlengths(Hsapiens)
valid_data <- data[start(data) >= 1 & end(data) <= chr_lengths[seqnames(data)]]

# 运行GADEM
gadem <- GADEM(valid_data, verbose=1, genome=Hsapiens)
plot(gadem)

2. 调整resize的锚点策略 #

如果不想丢弃区域，可以改用fix='start'或fix='end'进行resize，避免中心扩展时超出染色体边界：

# 以起始位置为锚点扩展
data = resize(data, width = 50, fix='start')
# 或以终止位置为锚点扩展
# data = resize(data, width = 50, fix='end')

3. 检查并修正原始Bed数据 #

直接查看第136行的原始Bed区域，确认其是否接近染色体末端。必要时手动调整该区域的宽度，或修改resize的宽度参数，避免超出染色体范围。

内容的提问来源于stack exchange，提问作者mehdi heidari