作業(yè)要求:
根據(jù)給定的兩個(gè)數(shù)據(jù)文件genefile.txt(基因表達(dá)譜)和geneid(基因ID號(hào))颅痊,利用熵值評(píng)價(jià)各個(gè)基因的表達(dá)穩(wěn)定性吕粹,并找出前10個(gè)表達(dá)最不穩(wěn)定的基因(可能差異表達(dá))描睦。
基本上就是對(duì)熵的理解疙描,和R中entropy包的使用
↓該加載加載管行,該讀數(shù)據(jù)讀數(shù)據(jù)
library(entropy)
gene_exp <- read.table("genefile.txt",sep=" ",header=T)
gene_exp <- as.matrix(gene_exp)
gene_id <- read.table("geneid.txt",header=T)
數(shù)據(jù)觀察
表達(dá)數(shù)據(jù)是54675*22格带,id是54675與表達(dá)數(shù)據(jù)行數(shù)匹配
再看看課件里有啥別的信息
表達(dá)數(shù)據(jù)是11列加11列
大體思路
1. 表達(dá)穩(wěn)定性(熵)如何得到
我是一個(gè)喜歡手寫的原始人
2. 將各自的熵和gene_id對(duì)上號(hào)
略
3. 提取前10不穩(wěn)定
按熵降序排序
代碼和結(jié)果
1. 表達(dá)穩(wěn)定性(熵)獲取
giveMeYourEntropy <- function(x){
n <- length(x)
dis1 <- discretize(x[1:(n/2)],numBins = 3,r=range(x[1:(n/2)]))
dis2 <- discretize(x[(n/2):n],numBins = 3,r=range(x[(n/2):n]))
freqs2d <- rbind(dis1,dis2)
en <- mi.plugin(freqs2d,unit="log2")
return (en)
}
雖然我是已知11+11啦囊蓝,但是寫成長(zhǎng)度一劈二去算兩個(gè)離散化后頻率其實(shí)是想體現(xiàn)我尊重代碼饿悬!在函數(shù)里寫一些已知的數(shù)字不會(huì)很沒有美感嘛哈哈哈。
2. 將各自的熵和gene_id對(duì)上號(hào)
gene_processed <- gene_id
gene_processed$entropy <- apply(gene_exp,1,giveMeYourEntropy)
眾所周知我是一條懶狗
沒必要把熵的信息再加到表達(dá)那個(gè)大大大矩陣屁股上聚霜。既然id和表達(dá)數(shù)據(jù)的行一一對(duì)應(yīng)狡恬,我為啥不直接把熵加載一列的id后面呢珠叔?
3. 提取前10不穩(wěn)定
sorted <- gene_processed[order(gene_processed$entropy,decreasing=T),]
head(sorted,10)
前10不穩(wěn)定