2016-07-21 1 views
1

특정 지점 또는 오히려 범위에서 발생하는 게놈의 돌연변이를 계산해야합니다. 돌연변이는 게놈 위치 (염색체 및 염기쌍, 예를 들어 Chr1, 10658324)를 갖는다. 범위 또는 지점은 각각 게놈의 주어진 위치의 상향 및 하향 (+ -) 10000 염기 쌍으로 정의됩니다. 돌연변이의 위치와 "스팟"의 위치는 모두 데이터 프레임에 저장됩니다.데이터 프레임에 주어진 게놈 영역 주변의 출현 횟수

예 :

그래서 요구하고있는 질문은
set.seed(1) 

Chr <- 1 
Pos <- as.integer(runif(5000 , 0, 1e8)) 
mutations <- data.frame(Pos, Chr) 

Chr <- 1 
Pos <- as.integer(runif(50 , 0, 1e8)) 
spots <- data.frame(Pos, Chr) 

: "장소"에 주어진 위치 주변 -10k + 존재하는 염기쌍이 얼마나 많은 돌연변이. (예 : 지점이 100k 인 경우 범위는 90k-110k입니다.) 실제 데이터에는 24 개의 염색체가 모두 포함되지만 당연히 하나의 염색체에만 집중할 수 있습니다. 최종 데이터에는 "spot"과 그 근처에있는 돌연변이 수, 이상적으로는 데이터 프레임 또는 매트릭스가 있어야합니다.

미리 제안 해 주셔서 감사합니다.


여기에는 첫 번째 시도가 있지만, 나는 그것을하는 더 우아한 방법이 있습니다.

w <- 10000 #setting range to 10k basepairs 
loop <- spots$Pos #creating vector of positions to loop through 
out <- data.frame(0,0) 
colnames(out) <- c("Pos", "Count") 

for (l in loop) { 
    temp <- nrow(filter(mutations, Pos>=l-w, Pos<=l+w)) 
    temp2 <- cbind(l,temp) 
    colnames(temp2) <- c("Pos", "Count") 
    out <- rbind(out, temp2) 
} 
out <- out[-1,] 
+0

는 R 커뮤니티의 도움을 얻으려면이 매우 특정이, 다음은 – Learner

+0

왜 연속 분포에서 의사 난수를 사용하여 이산 (정수) 분포에서 발생하고있는 것을 시뮬레이션하고 있습니까? "올바른"대답을 줄 수있는 예를 게시해야합니다. –

+1

유용한 설정 작업을 제공하는 Genomic Ranges를 살펴보십시오. https://bioconductor.org/packages/release/bioc/html/GenomicRanges.html – Drey

답변

3

, 다음 집계 data.table foverlaps를 사용 : 당신은 입력 및 예상 출력 예를 들어, 다음 사람을 제공하는 더 나은

library(data.table) 
#set the flank 
myFlank <- 100000 

#convert to ranges with flank 
spotsRange <- data.table(
    chr = spots$Chr, 
    start = spots$Pos - myFlank, 
    end = spots$Pos + myFlank, 
    posSpot = spots$Pos, 
    key = c("chr", "start", "end")) 

#convert to ranges start end same as pos 
mutationsRange <- data.table(
    chr = mutations$Chr, 
    start = mutations$Pos, 
    end = mutations$Pos, 
    key = c("chr", "start", "end")) 

#merge by overlap 
res <- foverlaps(mutationsRange, spotsRange, nomatch = 0) 

#count mutations 
resCnt <- data.frame(table(res$posSpot)) 
colnames(resCnt) <- c("Pos", "MutationCount") 
merge(spots, resCnt, by = "Pos") 
#   Pos Chr MutationCount 
# 1 3439618 1   10 
# 2 3549952 1   15 
# 3 4375314 1   11 
# 4 7337370 1   13 
# ... 
2

나는 R 침대 조작에 익숙하지 않은, 그래서 나는 GRANGES 또는 다른 R 생물 정보학 라이브러리로 변환 할 수 있습니다 여기 bedtools 누군가와 답변을 제안하겠습니다.

본질적으로, 당신은 두개의 침대 파일을 가지고 있습니다. 하나는 당신의 스팟들과 다른 하나는 돌연변이들입니다 (나는 후자에서 각각에 대해 1bp 좌표를 가정합니다). 이 경우 closestBed을 사용하여 가장 가까운 지점을 얻고 각 돌연변이의 거리를 bp로 지정한 다음 지점에서 10KB 인 것을 필터링합니다. UNIX 환경의 코드는 다음과 같이 보일 것입니다 : 열 아홉 ($9가) 가장 가까운 지점에서 BP의 거리가 될 것입니다

# Assuming 4-column file structure (chr start end name) 
closestBed -d -a mutations.bed -b spots.bed | awk '$9 <= 10000 {print}' 

. 좀 더 구체적으로하고 싶다면 매뉴얼 페이지 http://bedtools.readthedocs.io/en/latest/content/tools/closest.html을 확인하십시오. 나는 적어도 R에 bedtools-like 패키지가 하나 있다고 확신합니다. 기능이 비슷하다면 똑같은 솔루션을 적용 할 수 있습니다.

희망 하시겠습니까?