今日のDBCLS

統計解析ソフト「R」での立廻り 統合TV(togotv)｜生命科学系DB・ツール使い倒し系チャンネルを完成させました.

### mCGrate関数
mCGrate <- function(dat, begin=0, end=48129895, window.size=1000000) {
  dat.CG <- dat[dat$class == "CG", ] # classがCGの行のみを選択
  separators <- seq(begin, end, window.size) # 区切りを作る separators : [begin, begin + window.size, begin + 2*window.size, ...]
  mapply(function(sep) { # sepにseparatorsの要素を順次あてはめて、下の4行を実行
    window <- dat.CG[  sep < dat.CG$position
                     & dat.CG$position < sep + window.size, ] # ウィンドウへ切り出し 
    rate <- sum(window$mc) / sum(window$h) # メチル化率の割り算
    if (is.nan(rate)) NA else rate # 分母が0だとNaNになるが、これをNAに変換
  }, separators)
}

### ファイルの読み込み・メチル化率の計算
ADS <- mCGrate(read.csv("ADS", sep = "\t"))
ADS_adipose <- mCGrate(read.csv("ADS-adipose", sep = "\t"))
ADS_iPSC <- mCGrate(read.csv("ADS-iPSC", sep = "\t"))
FF_iPSC_19.11 <- mCGrate(read.csv("FF-iPSC-19.11", sep = "\t"))
FF_iPSC_19.11_BMP4 <- mCGrate(read.csv("FF-iPSC-19.11-BMP4", sep = "\t"))
FF_iPSC_19.7 <- mCGrate(read.csv("FF-iPSC-19.7", sep = "\t"))
FF_iPSC_6.9 <- mCGrate(read.csv("FF-iPSC-6.9", sep = "\t"))
FF <- mCGrate(read.csv("FF", sep = "\t"))
H1 <- mCGrate(read.csv("H1", sep = "\t"))
H1_BMP4 <- mCGrate(read.csv("H1-BMP4", sep = "\t"))
H9 <- mCGrate(read.csv("H9", sep = "\t"))
IMR90 <- mCGrate(read.csv("IMR90", sep = "\t"))
IMR90_iPSC <- mCGrate(read.csv("IMR90-iPSC", sep = "\t"))

### それぞれのメチル化率の列ベクトルを結合して行列に
dat <- cbind(ADS, ADS_adipose, ADS_iPSC, FF_iPSC_19.11, FF_iPSC_19.11_BMP4, FF_iPSC_19.7, FF_iPSC_6.9, FF, H1, H1_BMP4, H9, IMR90, IMR90_iPSC)

### 1 - 相関係数を距離として階層的クラスタリング・樹状図の描画
plot(hclust(as.dist(1 - cor(dat, use="complete.obs"))))