2016-07-21 3 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の染色体がすべて含まれていますが、ここでは簡単にするために1つの染色体に焦点を当てることができます。 最終データには、「スポット」とその近傍の突然変異の数、理想的にはデータフレームまたはマトリックスが含まれている必要があります。

アドバイスやアドバイスを事前にいただきありがとうございます。


は、ここで最初の試みだが、私はかなりSHURE午前、それを行う方法よりエレガントな方法があります。

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

有用なセット操作を提供するゲノム範囲を見てください: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のベッドの操作に慣れていないので、私はグランジや他のRのバイオインフォマティクスのライブラリに変換しようとすることができ、ここでbedtoolsと人との答えを提案するつもりです。

本質的に、2つのベッドファイルがあります.1つはあなたのスポットに、もう1つは突然変異に伴います(後者ではそれぞれ1bpの座標を仮定しています)。この場合、closestBedを使用して、各突然変異の最も近いスポットおよび距離をbpで取得し、スポットから10KBのものをフィルタリングします。コラム9($9)は、最も近い地点からbpの距離になります

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

:UNIX環境でのコードは次のようになります。より具体的な方法に応じて、マニュアルページのhttp://bedtools.readthedocs.io/en/latest/content/tools/closest.htmlを確認することができます。私はRにbedtoolsのようなパッケージが少なくとも1つあると確信しています。機能が似ている場合は、まったく同じソリューションを適用できます。

希望に役立ちます!

関連する問題