trace cell lineage - JSD

幫師姐算一個用poly-G trace cell lineage的數據（實際上就是在細胞裡面畫phylogeny tree），接觸到了JSD，Jensen–Shannon divergence。覺得這個統計量很有趣，拿出來分享一下。

簡單來說，他們實驗室用了一個辦法，給每個cell加上「追蹤器」，這個追蹤器稱為poly-G,在細胞複製的時候，GG...G的長度會隨機發生變化。比如題頭裡面的背景圖，展示了從一個zygote開始到不斷的二分產生的子代細胞，可以通過GG...G pattern的變化，來對細胞進行聚類。在腫瘤中的應用就是，可以知道不同tumor之間的「親緣」關係，比如轉移tumo是從哪一個原發tumor轉移出來的。

數據大概長這個樣子：

有6個sample和4個marker（mxx1, p-52, tol-9 和 hatyo）。

每個marker可能佔據多行，比如mxx1有4種allele type，那麼就佔據了4行。

每行中，後面6列，對應了6個sample中，這種allele的數目。

現在的問題是，從這個數據出發，最後構建出6個sample的phylogeny tree.

有很多種方法，可以將這面這個表格

轉化成一個distance matrix，然後再根據distance matrix來用hclust進行聚類或者是nj鄰接法進行畫樹。

根據PNAS（2014）和Science（2017）這兩篇paper，簡述一下他們的處理方法：

因為6個sample中有一個是normal tissue，其餘5個是mutate or tumor tissue（其實貌似這裡的normal tissue是「pseudo」的，是把其餘所有tumor tissue數據combine一起，「製造」出來的一個normal tissue）

1，對每一個marker，根據allele freq的分布，計算這個marker，每個mutate sample和normal tissue之間的相似度。

這樣就轉化成為這樣一個表格(假設sample1是normal or reference tissue）：

2，對於上面的表格(記為jsd.txt)構建距離矩陣，有很多種方法，比如R中的dist(),cov(),cor()等，

m<-as.matrix(cov(out,method="pearson"))

m <- as.matrix(dist(t(out),method = "euclidean"));

3，對得到的distance matrix進行聚類or畫樹

hclust(),或者ape包中的nj()

可以看到，2、3步都是常規分析，step1是分析的重頭戲。

在PNAS和Science的兩篇paper中，計算mutate sample和normal tissue之間的相似度時，採用了JSD 統計量。

JSD簡單來講，就是衡量兩個distribution之間相似度的統計量，因為用到了「熵Entropy」這個概念，顯得很fancy。

Stackoverflow上面有一個很好的計算JSD的R code: Jensen Shannon divergence in R

比如兩個分布p 和 q，現在想計算他倆之間的相似度，

# p & q are distributions so their elements should sum up to 1p <- c(0.00029421, 0.42837957, 0.1371827, 0.00029419, 0.00029419, 0.40526004, 0.02741252, 0.00029422, 0.00029417, 0.00029418)q <- c(0.00476199, 0.004762, 0.004762, 0.00476202, 0.95714168, 0.00476213, 0.00476212, 0.00476202, 0.00476202, 0.00476202)

先計算p和q單獨的熵，然後依據權重加和（如果兩個分布等權重，那就是 0.5*H(p) + 0.5*H(q), H(.)計算一個分布的熵的函數。

H <- function(v) { v <- v[v > 0] return(sum(-v * log(v)))}m <- 0.5 * (p + q)JS <- 0.5 * (sum(p * log(p / m)) + sum(q * log(q / m)))> JS[1] 0.6457538

然後按權重合併兩個p和q得到一個新的分布m：

m <- 0.5 * (p + q)

再對這個m計算熵值，JSD就是這兩個熵值之間的差。

JSD <- function(w, m) { return(H(m %*% w) - apply(m, 2, H) %*% w)}> JSD(w = c(1/3, 1/3, 1/3), m = cbind(p, q, m)) [,1][1,] 0.4305025

JSD作為一個量取兩個分布相似度的統計量，有蠻多不錯的性質，詳細的就不展開啦，這裡只是結合應用案例，簡單介紹一下這個統計量。

具體的分析流程可以看下面這兩篇文章，會有更多的技術細節。

Ref：

1， Naxerova, Kamila, et al. "Hypermutable DNA chronicles the evolution of human colon cancer."Proceedings of the National Academy of Sciences111.18 (2014): E1889-E1898.

2，Naxerova, Kamila, et al. "Origins of lymphatic and distant metastases in human colorectal cancer."Science357.6346 (2017): 55-60.