藥物基因組學_個體化實驗分析_實驗報告

一、實驗目的

1、（學會）收集特定癌種特定藥物的包含藥物響應與否相關的基因表達譜資料（TCGA/GEO）

2、（掌握）根據響應資訊將樣本劃分為耐藥或者敏感兩類

3、（掌握）構建分類器

4、（學會）利用構建的分類器實現個體化耐藥預測

二、實驗流程

（1）蒐集資料【資料使用的是平時研究的癌型】

從GEO網站上搜尋 Colorecal Cancer FOLFOX【結直腸癌藥物治療聯合方案】

在series_matrix檔案中查詢是否有藥物相應的資料，如下圖。

根據自己感興趣的方向所選擇的資料集GSE19860。同時還下載了GSE3964（5-FU），GSE72970（FOLFOX），但是由於在做T檢驗篩選差異基因的時候，後面兩套資料根據藥物是否敏感資訊進行分類得到的t檢驗結果並不顯著，因此，最後選擇了GSE19860進行以下的個體化耐藥分析實驗。

（2）根據響應資訊將樣本劃分為耐藥或者敏感兩類

在這套資料中，將FL_Responder定義為敏感標籤為0，FL_Non_Responder定義為耐藥標籤為1。

（3）構建分類器【以下使用的分類器是weka軟體中所構建的】

①篩選差異基因

程式碼：

setwd("E:\\...\\experiment")

rm(list=ls())

\######篩選差異基因

exp<-read.table("GSE19860_series_matrix.txt",header=T,sep="\t",stringsAsFactors = F)

\####補充缺失值

library(DMwR)

exp_new<-knnImputation(exp[,-1], k = 15, scale = T, meth = "weighAvg", distData = NULL)

\#exp<-exp[-which(is.na(exp[,1])),]  ####觀察到有缺失值，因此採用k均值補缺失值

\#length(which(is.na(exp[,1])))

labell<-read.csv("label.txt",sep="\t",header=F,stringsAsFactors = F)

geneid<-exp[,1] #####響應response敏感為0，no response不響應耐藥為1

exp_normal=exp_new[,which(labell[1,]==0)]

exp_cancer=exp_new[,which(labell[1,]==1)]

label<-c(rep(0,dim(exp_normal)[2]),rep(1,dim(exp_cancer)[2]))

\##rep()表示重複，重複1這個標籤exp_normal的列數（[2]表示列，[1]表示行）

exprs<-cbind(exp_normal,exp_cancer)

exp_case=exprs[,which(label==1)];

exp_control=exprs[,which(label==0)];

t_result=matrix(,length(geneid),4)

for(i in 1:nrow(exp)){

 T_test=t.test(exp_case[i,],exp_control[i,],alternative="two.sided",paired=FALSE)

 t_result[i,1:3]=c(geneid[i],T_test$statistic,T_test$p.value);

} 

t_result[,4]=p.adjust(t_result[,3],method="BH");

t_final<-t_result[t_result[,3]<0.05,1:4]  ###觀察到fdr<0.05的情況下，無結果，因此採用p<0.05進行篩選差異

②構建放入weka軟體的訓練集和驗證集檔案

#####構建放入weka的檔案

\####對於篩選出的差異的探針，將這些探針對應的在樣本中的表達譜提取出來

loc<-match(t_final[,1],exp[,1])

new<-exp[loc,]

new<-t(new) ####將資料設定為行為樣本，列為基因

colnames(new)<-new[1,]

new<-new[-1,]

label_new<-as.data.frame(t(labell))

all_differ<-cbind(new,label_new)

\####訓練集和測試集資料

traindata<-all_differ[c(1:20),]

testdata<-all_differ[c(21:40),]

\###選取20個樣本作為訓練集，20個樣本作為測試集。

write.csv(traindata,"traindata.csv")

write.csv(testdata,"testdata.csv")

#####接下來是在excel表格中進行的操作

訓練集和測試集的資料最後一列都是label的標籤，其中測試機的label要都換成英文的?,如下圖。

####放入weka將.csv的檔案儲存為.arff檔案

步驟：Explorer -> open file [開啟所要的訓練集和測試集檔案] ->save[儲存為.arrf檔案格式]

####對.arrf檔案在notepad+進行操作

將label後面的numeric都換成英文的{0,1}，這一步很重要，不然後續會報錯。

③接下來就可以放入weka進行分析

以下分析採取20個樣本作為訓練集，20個樣本作為測試集，效果是與其他不同樣本比例劃分所得到的結果相比較後是較好的了，因此採用1:1構建訓練集和測試集。

注：因為採取訓練集個數大於測試集個數，會導致分類器的分類結果過度擬合，而導致分類準確率下降，並且準確率都在50%。

\1. 隨機森林

在cmd（命令提示符）輸入如下程式碼：

java -Xmx1024m -classpath .;D:\Weka-3-6\weka.jar weka.classifiers.trees.RandomForest -t C:\Users\11833\Desktop\class_data\traindata.csv.arff -d colonRandomForest.model > colonRandomForest.out

以下操作也如上圖進行輸入，建議將weka裝在D盤，路徑都為英文，防止出錯。

\2. svm支援向量機

在cmd（命令提示符）輸入如下程式碼：

svm模型的構建

java -Xmx1024m -classpath .;D:\Weka-3-6\weka.jar;D:\Weka-3-6\libsvm\java\libsvm.jar weka.classifiers.functions.LibSVM  -t C:\Users\11833\Desktop\class_data\traindata.csv.arff  -d colonLibSVM.model > colonLibSVM.out

（4）利用構建的分類器實現個體化耐藥預測

根據前面的操作，我們會得到相應的分類器模型，接下來就是對測試結果進行驗證。

\1. 隨機森林

程式碼：

java -Xmx1024m -classpath .;D:\Weka-3-6\weka.jar weka.classifiers.trees.RandomForest -l colonRandomForest.model -T C:\Users\11833\Desktop\class_data\testdata.csv.arff -p 0 > PredictionResultRandomForest.txt

%%%%這個語句是用構建的模型來預測驗證集-l是載入儲存的模型，-T是驗證集 -p 0這個引數是給出每個樣本分類的標籤，和分為該類的可能性

%%%%驗證集的arff檔案和訓練集不一樣，他的類標籤相當於是未知的，在類標籤那列就用？表示，但是屬性那裡要注意，你這裡要和訓練集一樣

結果：

該資料是選擇了訓練集30個樣本，測試集10個樣本，根據比較，該分類器準確率只有50%。

\2. svm支援向量機

程式碼：

svm驗證集的驗證結果

java -Xmx1024m -classpath .;D:\Weka-3-6\weka.jar;D:\Weka-3-6\libsvm\java\libsvm.jar weka.classifiers.functions.LibSVM -l colonLibSVM.model -T C:\Users\11833\Desktop\class_data\testdata.csv.arff -p 0 > PredictionResultSVM.txt

其中分類正確的共有13個，分類錯誤的有7個，該分類器的正確率為65%。

三、討論

下載的三套資料都沒有得出較為差異的結果，可能是因為資料的分類未能很好將差異的基因篩選出來，且我同時注意到大概只剩下1000多個探針沒有對應多個基因ID，而其他探針大多都比對到了多個探針，因此仍然對資料結果有所保留，也可能因此導致了分類器不夠robust，所以分類效果不夠準確。

歡迎關注我的微信公眾號呀~