我想在我的数据集中创建一个新的列,其中的值由另一个数据集中的值确定,但这并不像一个列中的值是另一个列中值的函数那么简单。下面是一个例子:
>df1
chromosome position
1 1 1
2 1 2
3 1 4
4 1 5
5 1 7
6 1 12
7 1 13
8 1 15
9 1 21
10 1 23
11 1 24
12 2 1
13 2 5
14 2 7
15 2 8
16 2 12
17 2 15
18 2 18
19 2 21
20 2 22
和
>df2
chromosome segment_start segment_end segment.number
1 1 1 5 1.1
2 1 6 20 1.2
3 1 21 25 1.3
4 2 1 7 2.1
5 2 8 16 2.2
6 2 18 22 2.3
我想在df1中创建一个名为'segment'的新列,segment中的值将由'position'中的值属于哪个段(由df2中的'segment_start', 'segment_end'和'chromosome'确定)来确定。例如,在df1中,第7行,位置=13,染色体=1。因为13介于6和20之间,我假设的"segment"列中的条目将是1.2,来自df2的第2行,因为13位于该行的segment_start和segment_end之间(分别为6和20),并且df1第7行的"染色体"值为1,就像df2第2行的"染色体"值为1一样。
df1中的每一行都属于df2中描述的一个段;也就是说,它与其中一个片段位于同一染色体上,其"位置"为>=segment_start和<=segment_end。我想把这些信息输入到df1中,它表示每个位置属于哪个段。
我正在考虑使用if函数,并开始使用:
if(df1$position>=df2$segment_start & df1$position<=df2$segment_end & df1$chromosome==df2$chromosome) df1$segment<-df2$segment.number
但我不确定这种方法是否可行。如果没有别的,也许代码可以帮助说明我要做的是什么。基本上,我想将每一行的位置和染色体匹配到df2中的一个片段。谢谢。
这似乎是一个滚动连接。对于
,您可以使用data.table
require(data.table)
DT1 <- data.table(df1, key = c('chromosome','position'))
DT2 <- data.table(df2, key = c('chromosome','section_start'))
# this will perform the join you want (but retain all the
# columns with names names of DT2)
# DT2[DT1, roll=TRUE]
# which is why I have renamed and subset here)
DT2[DT1, roll=TRUE][ ,list(chromosome,position = segment_start,segment.number)]
# chromosome position segment.number
# 1: 1 1 1.1
# 2: 1 2 1.1
# 3: 1 4 1.1
# 4: 1 5 1.1
# 5: 1 7 1.2
# 6: 1 12 1.2
# 7: 1 13 1.2
# 8: 1 15 1.2
# 9: 1 21 1.3
# 10: 1 23 1.3
# 11: 1 24 1.3
# 12: 2 1 2.1
# 13: 2 5 2.1
# 14: 2 7 2.1
# 15: 2 8 2.2
# 16: 2 12 2.2
# 17: 2 15 2.2
# 18: 2 18 2.3
# 19: 2 21 2.3
# 20: 2 22 2.3
您确实需要查看Bioconductor的GenomicRanges包。它提供了适合您的用例的数据结构。
首先,我们创建GRanges
对象:
gr1 <- with(df1, GRanges(chromosome, IRanges(position, width=1L)))
gr2 <- with(df2, GRanges(chromosome, IRanges(segment_start, segment_end),
segment.number=segment.number))
然后我们找到重叠并进行合并:
hits <- findOverlaps(gr1, gr2)
gr1$segment[queryHits(hits)] <- gr2$segment.number[subjectHits(hits)]
我将假设df2
中的区域是不重叠的,连续的和完整的(不缺少df1
中的任何位置)。我似乎每次尝试的方法都不一样,所以这是我的最新想法。
首先,确保染色体在两个数据集中都是一个因素
df1$chromosome<-factor(df1$chromosome)
df2$chromosome<-factor(df2$chromosome)
现在我想将chr/pos展开为一个通用位置,我将使用
ends<-with(df2, tapply(segment_end, chromosome, max))
offset<-head(c(0,cumsum(ends)),-1)
names(offset)<-names(ends)
这将为所有染色体上的所有位置分配唯一的位置值,并在这个新系统中跟踪每个染色体开始的偏移量。现在我们将根据df2
seglookup <- approxfun(with(df2, offset[chromosome]+segment_start), 1:nrow(df2),
method="constant", rule=2)
我们使用approxfun
在每个片段的遗传位置空间中找到合适的区间。现在我们在df1
segid <- with(df1, seglookup(offset[chromosome]+position))
现在我们有了每个位置的正确ID。我们可以通过合并数据并查看结果
来验证这一点cbind(df1, df2[segid,-1])
chromosome position segment_start segment_end segment.number
1 1 1 1 5 1.1
2 1 2 1 5 1.1
3 1 4 1 5 1.1
4 1 5 1 5 1.1
5 1 7 6 20 1.2
6 1 12 6 20 1.2
7 1 13 6 20 1.2
8 1 15 6 20 1.2
9 1 21 21 25 1.3
10 1 23 21 25 1.3
11 1 24 21 25 1.3
12 2 1 1 7 2.1
13 2 5 1 7 2.1
14 2 7 1 7 2.1
15 2 8 8 16 2.2
16 2 12 8 16 2.2
17 2 15 8 16 2.2
18 2 18 18 22 2.3
19 2 21 18 22 2.3
20 2 22 18 22 2.3
看起来还行