Ошибка выбора данных

Question

При попытке выбрать данные из большого набора данных (d1) на основе небольшого набора данных (d2) возникает ошибка. Ниже мой сценарий и проблема.

**d1 <- read.table("MSv25.txt",header=T) 

d2 <- read.table("Flairall.Gene.txt",header=T) 

d2$low <- d2$start-10000 ; d2$high <- d2$end+10000 

d1$matched <- apply(d1,1,function(p) which(p['POS'] >=d2[,'low'] & p['POS'] <= d2[,'high'] & p['CHR']==d2[,"chromosome"])) 

d3 <- cbind(d1[which(d1$matched >0),], d2[unlist(d1$matched[which(d1$matched>0)]),]) 

write.table(d3,file="Flairall.GOBSgene.txt",quote=FALSE,sep="\t",row.names=FALSE,col.names=TRUE)** 

Д1, как это:

SNP CHR POS A1 A2 OR P 
rs1000007 2 237416793 C T 0.9785 0.4868 
rs1000003 3 99825597 G A 0.9091 0.009774 
rs1000002 3 185118462 C T 1.0111 0.6765 
rs10000012 4 1347325 G C 1.0045 0.9087 
rs10000042 4 5288053 T C 1.0622 0.3921 
rs10000062 4 5305645 G C 1.0116 0.779 
rs10000132 4 7450570 T C 0.9734 0.4748 
rs10000081 4 16957461 C T 1.0288 0.3585 
rs10000100 4 19119591 A G 1.0839 0.1417 
rs10000010 4 21227772 C T 0.971 0.2693 
rs10000092 4 21504615 C T 1.0342 0.27 
rs10000068 4 36600682 T C 1.055 0.103 
rs10000013 4 36901464 C A 1.0198 0.5379 
rs10000037 4 38600725 A G 1.0249 0.4217 
rs10000017 4 84997149 T C 0.9576 0.1912 
rs10000109 4 91586292 A T 0.9956 0.9349 
rs10000023 4 95952929 T G 0.9998 0.9951 
rs10000030 4 103593179 A G 1.0839 0.04208 
rs10000111 4 107137517 A G 1.0812 0.3128 
rs10000124 4 109325900 A C 1.0642 0.1906 
rs10000064 4 128029071 C T 1.0388 0.1578 
rs10000029 4 138905074 C T 0.7832 0.14 
rs10000036 4 139438712 T C 0.9848 0.5683 
rs10000033 4 139819348 C T 0.9918 0.7542 
rs10000121 4 157793485 A G 1.0008 0.9769 
rs10000041 4 165841405 G T 1.0042 0.9146 
rs10000082 4 167529642 T C 0.9733 0.6612 

Д2, как это:

gene start end chromosome 
WFDC9 237416000 237418000 2 
SRGAP3 19119590 21504615 4 

В общем, я хочу, чтобы выбрать SNPs внутри генов, расширяя окно 10kb в начальной и конечной позиции.

Вот мои результаты сценария:

SNP CHR  POS A1 A2  OR  P matched gene  start  end chromosome  low  high 
    1 rs1000007 2 237416793 C T 0.9785 0.4868  1 WFDC9 237416000 237418000   2 237406000 237428000 

Который не является правильным, поскольку есть один ген отсутствует. Правильный один должен быть:

gene start end chromosome SNP CHR POS A1 A2 OR P 
WFDC9 237416000 237418000 2 rs1000007 2 237416793 C T 0.9785 0.4868 
SRGAP3 19119590 21504615 4 rs10000100 4 19119591 A G 1.0839 0.1417 
SRGAP3 19119590 21504615 4 rs10000010 4 21227772 C T 0.971 0.2693 
SRGAP3 19119590 21504615 4 rs10000092 4 21504615 C T 1.0342 0.27 

Может кто-нибудь помочь мне указать на то, что это неправильно .... Большое спасибо ...

Answer 1

Ваш код, кажется, полностью назад к тому, что вы пытаетесь достичь:

«Для каждого гена (в d2), какие SNP (из d1) находятся в пределах 10kb этого гена?»

Прежде всего, ваш код для d1$matched назад. Все ваши p и d2 s должны быть в обратном порядке (в настоящее время это не имеет смысла?), Давая вам список SNP, которые находятся в цис с каждым геном (+/- 10kb).

я бы подойти к нему так, как я уже Выражаясь вашим вопрос:

cisWindow <- 10000 # size of your +/- window, in this case 10kb. 
d3 <- data.frame() 
# For each gene, locate the cis-SNPs 
for (i in 1:nrow(d2)) { 
    # Broken down into steps for readability. 
    inCis <- d1[which(d1[,"CHR"] == d2[i, "chromosome"]),] 
    inCis <- inCis[which(inCis[,"POS"] >= (d2[i, "start"] - cisWindow)),] 
    inCis <- inCis[which(inCis[,"POS"] <= (d2[i, "end"] + cisWindow)),] 
    # Now we have the cis-SNPs, so lets build the data.frame for this gene, 
    # and grow our data.frame d3: 
    if (nrow(inCis) > 0) { 
    d3 <- rbind(d3, cbind(d2[i,], inCis)) 
    } 
} 

Я пытался найти решение, которое не затрагивало растет d3 в цикле, а потому, что вы прикрепление каждой строки d2 до 0 или более строк от d1 Мне не удалось найти решение, которое не является ужасно неэффективным.