### Create: Jianming Zeng
### Date: 2019-04-02 21:59:01
### Email: jmzeng1314@163.com
rm(list=ls())
options(stringsAsFactors = F)
# 注意，并不是說使用 RTCGA.miRNASeq包的數(shù)據(jù)是最佳選擇，只是因為這個演示起來最方便。
# 因為GDC官網(wǎng)下載數(shù)據(jù)具有一定門檻，也不是每個人都必須學(xué)會的。
getwd()
Rdata_dir='../Rdata/'
Figure_dir='../figures/'

# 如果開啟下面代碼，就會從RTCGA.miRNASeq包里面提取miRNA表達(dá)矩陣和對應(yīng)的樣本臨床信息。 
if(F){
  library(RTCGA.miRNASeq)  #加載包 TCGA的miRNASeq數(shù)據(jù)就在該包里      安裝見http://m.itdecent.cn/p/a3c64ee1f63b
  s=rownames(KIRC.miRNASeq)[seq(1,nrow(KIRC.miRNASeq),by=3)] #觀察KIRC.miRNASeq行名發(fā)現(xiàn)以3為規(guī)律，seq提取行名
  expr <- expressionsTCGA(KIRC.miRNASeq) #獲取表達(dá)矩陣
  dim(expr)
  expr[1:40,1:4]
  expr=as.data.frame(expr[seq(1,nrow(expr),by=3),3:ncol(expr)]) #根據(jù)之前提取的行名s過濾表達(dá)矩陣
  mi=colnames(expr)
  expr=apply(expr,1,as.numeric)  #將行名設(shè)為數(shù)值型
  colnames(expr)=s #將表達(dá)矩陣列名轉(zhuǎn)換為s
  rownames(expr)=mi #將表達(dá)矩陣列名轉(zhuǎn)換為s
  expr[1:4,1:4]
  expr=na.omit(expr) 
  dim(expr)  #此時有1046行，593列，說明沒有缺失值
  expr=expr[apply(expr, 1,function(x){sum(x>1)>10}),] 
  dim(expr)    # 552 593  這樣就獲得了593個樣本對應(yīng)的552個miRNA信息

  library(RTCGA.clinical) 
  meta <- KIRC.clinical
  tmp=as.data.frame(colnames(meta))
  meta[(grepl('patient.bcr_patient_barcode',colnames(meta)))]
  meta[(grepl('patient.days_to_last_followup',colnames(meta)))]
  meta[(grepl('patient.days_to_death',colnames(meta)))]
  meta[(grepl('patient.vital_status',colnames(meta)))]
  meta=as.data.frame(meta[c('patient.bcr_patient_barcode','patient.vital_status',
                            'patient.days_to_death','patient.days_to_last_followup',
                            'patient.race',
                            'patient.age_at_initial_pathologic_diagnosis',
                            'patient.gender' ,
                           'patient.stage_event.pathologic_stage')])
# 提取patient.bcr_patient_barcode，patient.vital_status等樣本信息
  #meta[(grepl('patient.stage_event.pathologic_stage',colnames(meta)))]
  ## 每次運(yùn)行代碼，就會重新生成文件。
  save(expr,meta,
       file = file.path(Rdata_dir,'TCGA-KIRC-miRNA-example.Rdata')
         )
}

## 我們已經(jīng)運(yùn)行了上面被關(guān)閉的代碼，而且保存了miRNA表達(dá)矩陣和對應(yīng)的樣本臨床信息
# 現(xiàn)在直接加載即可。
load( file = 
        file.path(Rdata_dir,'TCGA-KIRC-miRNA-example.Rdata')
)
dim(expr)
dim(meta)
# 可以看到是 537個病人，但是有593個樣本，每個樣本有 552個miRNA信息。
# 當(dāng)然，這個數(shù)據(jù)集可以下載原始測序數(shù)據(jù)進(jìn)行重新比對，可以拿到更多的miRNA信息

表達(dá)矩陣

臨床信息

參考來源：生信技能樹

友情鏈接：

課程分享
生信技能樹全球公益巡講
（https://mp.weixin.qq.com/s/E9ykuIbc-2Ja9HOY0bn_6g）
B站公益74小時生信工程師教學(xué)視頻合輯
（https://mp.weixin.qq.com/s/IyFK7l_WBAiUgqQi8O7Hxw）
招學(xué)徒：
（https://mp.weixin.qq.com/s/KgbilzXnFjbKKunuw7NVfw）

歡迎關(guān)注公眾號：青島生信菜鳥團(tuán)