基于重叠合并 2 个基因组文件
merge 2 genome files based on overlap
我有两个 "BED" 文件。每个文件都指定了基因组中的一组区域（开始和结束列），并且每个文件都指定了给定基因组区域的特征（例如 NRL 和其他 returns 这些区域的 'mappability'）
它们的组织如下：
head(file1)
chr start end mappability
chr1 3000066 3000100 1.0000
chr1 3000100 3000130 0.5000
chr1 3000130 3000199 0.0625
chr1 3000199 3000277 0.0500
head(file2)
chr start end NRL
chr1 3000000 3000067 250
chr1 3000067 3000079 300
chr1 3000079 3000084 200
chr1 3000084 3000099 130
问题是这些文件是不同实验的结果，并不是两个文件之间记录的所有区域都重叠...因此我需要找出哪些区域重叠...
到目前为止，我的尝试是：
file1-read.table("file1.txt", sep='\t', header = F)
file2=read.table("file2.txt", sep='\t', header = F)
overlapping_regions<-function(file1, file2){
for(i in file1[,2]){
x<-seq(file1[i,2], file1[i,3])
for(j in file1[,2]){
y<-seq(file2[j,2], file2[j,3])
if(setequal(union(x, y), c(setdiff(x, y), intersect(x, y), setdiff(y, x)))){
####GET OVERLAP
}
}
}
}
上述策略的第一个问题是出现上述错误：
Error in seq.default(file1[i, 2], file1[i, 3]) :
'from' 不能为 NA、NaN 或无穷大
其次，我不确定该策略是否正确，因为我希望将每个文件的每一行与另一行进行比较，以找到重叠的 ANY 区域...
所以我想知道是否有人可以用 R 脚本帮助我解析这些文件，这样我就可以制作一个新文件，其中包含每个开始和结束指定列之间的重叠区域，并保留与每个原始文件...
所以我希望我的输出是这样的：
head(files_merged)
chr overlap mappability NRL GC_content more_features......
chr1 start-end 1.0000 250
chr1 start-end 0.5000 300
chr1 start-end 0.0625 200
我问这个是为了尝试应用机器学习算法来预测基因组特征。
我可以（很明显地）看出我的计划有多么有缺陷，因为一个文件中指定的区域可能比另一个文件中指定的区域小得多。因此，我也愿意接受有关更好的方法的建议？
这可能有点长，但您可以尝试一下。
我创建了类似的数据框，但不完全相同：
df1 <- data.frame(chr=rep("chr1",4),
start=c(100,200,300,400),
end=c(200,300,400,500),
mappability=c(1,0.5,0.0625,0.05))
df2 <- data.frame(chr=rep("chr1",4),
start=c(90,190,290,380),
end=c(120,220,320,390),
NRL=c(250,300,200,130))
加载使用映射和嵌套函数所需的库：
library(purrr)
library(tidyr)
一个带有开始和结束的小标题的函数，在 df1 中查找重叠的索引和 return 行号。
您可以根据您的边界、约束或重叠定义来编辑条件：
xx <- function(x){
y <- (x$start<df1$start & x$end<df1$end & x$end>df1$start) | (x$start>df1$start & x$start<df1$start & x$end>df1$end)
z <- which(y==TRUE)
ifelse((length(z)>0),z,0) %>%
as.integer()
}
嵌套 df2 并将 start-end 放在一个小标题中：
df2 <- df2 %>%
nest(start,end,.key=data.df2)
# A tibble: 4 x 3
chr NRL data.df2
<fctr> <dbl> <list>
1 chr1 250 <tibble [1 x 2]>
2 chr1 300 <tibble [1 x 2]>
3 chr1 200 <tibble [1 x 2]>
4 chr1 130 <tibble [1 x 2]>
将每一行中的小标题传递给函数 xx，它将 return 重叠的行 （如果存在多个条目的情况，则函数可能需要更改并且我们将使用 map 而不是 map_int)
df2 <- df2 %>%
mutate(idx=map_int(data.df2,xx)) %>%
unnest %>%
filter(idx!=0)
取消嵌套并删除没有交集的行后，我们将在 df2 中获得与 df1 中的条目重叠的条目。
# A tibble: 3 x 5
chr NRL idx start end
<fctr> <dbl> <int> <dbl> <dbl>
1 chr1 250 1 90 120
2 chr1 300 2 190 220
3 chr1 200 3 290 320
我们将向 df1 添加一个 idx 列以便能够合并：
df1 <- df1 %>%
变异（idx=seq_along（df1））
chr start end mappability idx
1 chr1 100 200 1.0000 1
2 chr1 200 300 0.5000 2
3 chr1 300 400 0.0625 3
4 chr1 400 500 0.0500 4
现在合并 df1 和 df2，基于索引：
df_all <- merge(df1,df2,by=c("idx"),
all.x = FALSE,
all.y = TRUE
)
TOu 会有类似的东西，您可以在其中清理和计算每一行的重叠：
idx chr.x start.x end.x mappability chr.y NRL start.y end.y
1 1 chr1 100 200 1.0000 chr1 250 90 120
2 2 chr1 200 300 0.5000 chr1 300 190 220
3 3 chr1 300 400 0.0625 chr1 200 290 320
Bioconductor support site 上也有人提出了这个问题，我在其中提供了类似的长答案。 @OmaymaS提供的数据结果为
> olaps
GRanges object with 6 ranges and 2 metadata columns:
seqnames ranges strand | mappability NRL
<Rle> <IRanges> <Rle> | <numeric> <numeric>
[1] chr1 [101, 120] * | 1 250
[2] chr1 [191, 200] * | 1 300
[3] chr1 [201, 220] * | 0.5 300
[4] chr1 [291, 300] * | 0.5 200
[5] chr1 [301, 320] * | 0.0625 200
[6] chr1 [381, 390] * | 0.0625 130
-------
seqinfo: 1 sequence from an unspecified genome; no seqlengths
将 BED 文件的基于 0 的半开区间转换为更友好的/Bioconductor 标准基于 1 的闭区间的基于 1 的偏移量。
我有两个 "BED" 文件。每个文件都指定了基因组中的一组区域（开始和结束列），并且每个文件都指定了给定基因组区域的特征（例如 NRL 和其他 returns 这些区域的 'mappability'）
它们的组织如下：
head(file1)
chr start end mappability
chr1 3000066 3000100 1.0000
chr1 3000100 3000130 0.5000
chr1 3000130 3000199 0.0625
chr1 3000199 3000277 0.0500
head(file2)
chr start end NRL
chr1 3000000 3000067 250
chr1 3000067 3000079 300
chr1 3000079 3000084 200
chr1 3000084 3000099 130
问题是这些文件是不同实验的结果，并不是两个文件之间记录的所有区域都重叠...因此我需要找出哪些区域重叠...
到目前为止，我的尝试是：
file1-read.table("file1.txt", sep='\t', header = F)
file2=read.table("file2.txt", sep='\t', header = F)
overlapping_regions<-function(file1, file2){
for(i in file1[,2]){
x<-seq(file1[i,2], file1[i,3])
for(j in file1[,2]){
y<-seq(file2[j,2], file2[j,3])
if(setequal(union(x, y), c(setdiff(x, y), intersect(x, y), setdiff(y, x)))){
####GET OVERLAP
}
}
}
}
上述策略的第一个问题是出现上述错误：
Error in seq.default(file1[i, 2], file1[i, 3]) :
'from' 不能为 NA、NaN 或无穷大
其次，我不确定该策略是否正确，因为我希望将每个文件的每一行与另一行进行比较，以找到重叠的 ANY 区域...
所以我想知道是否有人可以用 R 脚本帮助我解析这些文件，这样我就可以制作一个新文件，其中包含每个开始和结束指定列之间的重叠区域，并保留与每个原始文件...
所以我希望我的输出是这样的：
head(files_merged)
chr overlap mappability NRL GC_content more_features......
chr1 start-end 1.0000 250
chr1 start-end 0.5000 300
chr1 start-end 0.0625 200
我问这个是为了尝试应用机器学习算法来预测基因组特征。
我可以（很明显地）看出我的计划有多么有缺陷，因为一个文件中指定的区域可能比另一个文件中指定的区域小得多。因此，我也愿意接受有关更好的方法的建议？
这可能有点长，但您可以尝试一下。
我创建了类似的数据框，但不完全相同：
df1 <- data.frame(chr=rep("chr1",4),
start=c(100,200,300,400),
end=c(200,300,400,500),
mappability=c(1,0.5,0.0625,0.05))
df2 <- data.frame(chr=rep("chr1",4),
start=c(90,190,290,380),
end=c(120,220,320,390),
NRL=c(250,300,200,130))
加载使用映射和嵌套函数所需的库：
library(purrr)
library(tidyr)
一个带有开始和结束的小标题的函数，在 df1 中查找重叠的索引和 return 行号。 您可以根据您的边界、约束或重叠定义来编辑条件：
xx <- function(x){
y <- (x$start<df1$start & x$end<df1$end & x$end>df1$start) | (x$start>df1$start & x$start<df1$start & x$end>df1$end)
z <- which(y==TRUE)
ifelse((length(z)>0),z,0) %>%
as.integer()
}
嵌套 df2 并将 start-end 放在一个小标题中：
df2 <- df2 %>%
nest(start,end,.key=data.df2)
# A tibble: 4 x 3
chr NRL data.df2
<fctr> <dbl> <list>
1 chr1 250 <tibble [1 x 2]>
2 chr1 300 <tibble [1 x 2]>
3 chr1 200 <tibble [1 x 2]>
4 chr1 130 <tibble [1 x 2]>
将每一行中的小标题传递给函数 xx，它将 return 重叠的行 （如果存在多个条目的情况，则函数可能需要更改并且我们将使用 map 而不是 map_int)
df2 <- df2 %>%
mutate(idx=map_int(data.df2,xx)) %>%
unnest %>%
filter(idx!=0)
取消嵌套并删除没有交集的行后，我们将在 df2 中获得与 df1 中的条目重叠的条目。
# A tibble: 3 x 5
chr NRL idx start end
<fctr> <dbl> <int> <dbl> <dbl>
1 chr1 250 1 90 120
2 chr1 300 2 190 220
3 chr1 200 3 290 320
我们将向 df1 添加一个 idx 列以便能够合并：
df1 <- df1 %>% 变异（idx=seq_along（df1））
chr start end mappability idx
1 chr1 100 200 1.0000 1
2 chr1 200 300 0.5000 2
3 chr1 300 400 0.0625 3
4 chr1 400 500 0.0500 4
现在合并 df1 和 df2，基于索引：
df_all <- merge(df1,df2,by=c("idx"),
all.x = FALSE,
all.y = TRUE
)
TOu 会有类似的东西，您可以在其中清理和计算每一行的重叠：
idx chr.x start.x end.x mappability chr.y NRL start.y end.y
1 1 chr1 100 200 1.0000 chr1 250 90 120
2 2 chr1 200 300 0.5000 chr1 300 190 220
3 3 chr1 300 400 0.0625 chr1 200 290 320
Bioconductor support site 上也有人提出了这个问题，我在其中提供了类似的长答案。 @OmaymaS提供的数据结果为
> olaps
GRanges object with 6 ranges and 2 metadata columns:
seqnames ranges strand | mappability NRL
<Rle> <IRanges> <Rle> | <numeric> <numeric>
[1] chr1 [101, 120] * | 1 250
[2] chr1 [191, 200] * | 1 300
[3] chr1 [201, 220] * | 0.5 300
[4] chr1 [291, 300] * | 0.5 200
[5] chr1 [301, 320] * | 0.0625 200
[6] chr1 [381, 390] * | 0.0625 130
-------
seqinfo: 1 sequence from an unspecified genome; no seqlengths
将 BED 文件的基于 0 的半开区间转换为更友好的/Bioconductor 标准基于 1 的闭区间的基于 1 的偏移量。