高表達的PVT1(lncRNA)能夠獨立且有效地預測葡萄膜黑色素瘤生存情況

生物_醫(yī)藥_科研 2019-04-16

展開全文

你現(xiàn)在看到的，仍然是學徒大作，上個星期我?guī)淼耐诰蜻^程展現(xiàn)非常的受歡迎，所以再接再厲，希望大家喜歡

5分的生信數(shù)據(jù)挖掘文章全部圖表復現(xiàn)（純R代碼）

文章鏈接

https:///10.1371/journal.pone.0189675

背景介紹

葡萄膜黑色素瘤是成人最常見的原發(fā)性眼內惡性腫瘤，在國外其發(fā)病率占眼內腫瘤的首位，在國內則僅次于視網(wǎng)膜母細胞瘤，居眼內腫瘤的第二位。主要經血循環(huán)轉移，全身轉移主要見于肝臟，肝轉移所致死亡率高達50%。

文章結論

高表達的lncRNA-PVT1能夠獨立且有效地預測葡萄膜黑色素瘤（uveal melanoma, UVM）患者的差預后。

由于最近接觸TCGA數(shù)據(jù)比較多，每次都要去XENA官網(wǎng)上下載很麻煩，就解析了一下官網(wǎng)。

保存一下Xena的網(wǎng)頁你會發(fā)現(xiàn)，源代碼全是js，找不要想要的數(shù)據(jù)，搜索發(fā)現(xiàn)了“Selenium”這個強大的包，在正式處理這次數(shù)據(jù)之前，先介紹一下提取的結果。

首先安裝以下軟件：

1 “JDK”[ https://download.oracle.com/otn-pub/java/jdk/12+33/312335d836a34c7c8bba9d963e26dc23/jdk-12_windows-x64_bin.exe ]

2 “selenium” [https://selenium-release.storage.googleapis.com/3.9/selenium-server-standalone-3.9.1.jar ]

3 瀏覽器驅動任意選擇一個(我的瀏覽器都是最新版，所以驅動我也選擇了最新的)：

Chrome: http://npm./mirrors/chromedriver/2.9/

Firefox: https://github.com/mozilla/geckodriver/releases/download/v0.24.0/geckodriver-v0.24.0-win64.zip

注：下載得到的都是exe，需要移動到你的瀏覽器安裝目錄下，并且將瀏覽器安裝目錄添加到你的環(huán)境變量。

4 “RSelenium”是一個R包，直接install.packages就可以安裝。

然后就可以在R里處理頁面了

library(RCurl)
library(XML)
library(RSelenium)
## 不想專門打開一個cmd，就在R里直接運行了，這些是我安裝軟件的路徑，改成你自己的就可以啦
system(paste(''C:\Program Files\Java\jdk-10.0.2\bin\java.exe'',
             '-Dwebdriver.chrome.driver='C:\Program Files (x86)\Google\Chrome\Application\chromedriver.exe'',
             '-jar 'C:\Program Files\Java\selenium-server-standalone-3.9.1.jar''),
       wait = FALSE)

## 最開始沒有找到 “RSelenium”包，所以是用 XML解析下面的網(wǎng)址，F(xiàn)12可以看到元素，我就

## 直接復制出了源代碼，在文件開頭加上“<?xml version='1.0'?>”## https:///datapages/?hub=https://tcga.:44

if (!file.exists( './data/xena_url.Rdata' )) {
  webpage <- readLines('./Xena.htm')
  pagetree <- htmlTreeParse(webpage)
  # parse the tree by tables
  a <- getNodeSet(pagetree$children$html, '/html/body/li/a')
  href <- sapply(a, function(el) xmlGetAttr(el, 'href'))
  href <- paste('https:///datapages/', href, sep = '')
  abbr <- getNodeSet(pagetree$children$html, '/html/body/li/a/abbr')
  abbr <- as.data.frame(unlist(abbr))
  abbr <- abbr[seq(3, nrow(abbr), by = 3), ]
  xena_url <- data.frame(url = href, tumor = abbr)
  save( xena_url, file = './data/xena_url.Rdata' )
}else{
  load('./data/xena_url.Rdata')
}

## 得到的對應關系如下：

如果上面這一段教程你們沒有看懂，就放棄吧，是我刻意在秀技巧，總之，前面我們講解過了

太多的TCGA數(shù)據(jù)庫下載，你肯定學會了！

5分的生信數(shù)據(jù)挖掘文章全部圖表復現(xiàn)（純R代碼）

## 下面就要開始下載我們這次的UVM數(shù)據(jù)了

tumor_type <- 'UVM'
url <- xena_url[grep(tumor_type, xena_url$tumor), 1]

remDr <- remoteDriver(browserName = 'firefox',
                      remoteServerAddr = 'localhost',
                      port = 4444L)
remDr$open()
remDr$navigate(url)

webElems <- remDr$findElements(using = 'css selector', 'a')
sub_type <- unlist(lapply(webElems, function(x) {x$getElementText()}))
sub_type
targetElem <- webElems[[grep('Phenotypes', sub_type)]]
targetElem$clickElement()
Sys.sleep(5)
webElems <- remDr$findElements(using = 'css selector', 'a')
sub_url <- unlist(lapply(webElems, function(x) {x$getElementText()}))
phenoData_url <- sub_url[grep('http', sub_url)[1]]
phenoData_url
remDr$goBack()

webElems <- remDr$findElements(using = 'css selector', 'a')
sub_type <- unlist(lapply(webElems, function(x) {x$getElementText()}))
sub_type
targetElem <- webElems[[grep('IlluminaHiSeq', sub_type)[2]]]
targetElem$clickElement()
Sys.sleep(5)
webElems <- remDr$findElements(using = 'css selector', 'a')
sub_url <- unlist(lapply(webElems, function(x) {x$getElementText()}))
raw_data_url <- sub_url[grep('http', sub_url)[1]]
raw_data_url

## 最終得到了兩個URL，臨床數(shù)據(jù)的下載地址phenoData_url和表達矩陣的下載地址raw_data_url

表達矩陣和臨床信息

接下來就是一般的流程了

tumor_type <- 'UVM'
Rdata_file <- paste('./data/TCGA.', tumor_type, '.Rdata', sep = '')
if (!file.exists( Rdata_file )) {
  url <- raw_data_url
  url_word <- strsplit(url, split = '/')
  gz_file <- url_word[[1]][length(url_word[[1]])]
  destfile <- paste('./raw_data/TCGA.', tumor_type, '.', gz_file, sep = '')
  download.file(url, destfile)
  library(R.utils)
  gunzip(destfile, remove = F)
  library(data.table)
  raw_data <- fread( destfile, sep = '	', header = T)
  raw_data <- as.data.frame( raw_data )
  raw_data[1:5, 1:6] 
  rownames( raw_data ) <- raw_data[, 1]
  raw_data <- raw_data[, -1]
  raw_data[1:5, 1:6]
  raw_data <- 2^raw_data - 1
  raw_data <- ceiling( raw_data )
  raw_data[1:5, 1:6]
  pick_row <- apply( raw_data, 1, function(x){
    sum(x == 0) < 10
  })
  raw_data <- raw_data[pick_row, ]
  dim(raw_data  )
  save( raw_data, file = Rdata_file )
}else{
  load( Rdata_file )
}


Rdata_file <- paste('./data/TCGA.', tumor_type, '.phenoData.Rdata', sep = '')
if (!file.exists( Rdata_file )) {
  url <- phenoData_url
  url_word <- strsplit(url, split = '/')
  gz_file <- url_word[[1]][length(url_word[[1]])]
  destfile <- paste('./raw_data/TCGA.', gz_file, sep = '')
  download.file(url, destfile)
  phenoData <- read.table( destfile,
                           header = T,
                           sep = '	',
                           quote = '' )
  name <- phenoData[ , 1]
  phenoData <- phenoData[ , -1]
  rownames( phenoData ) <- name
  phenoData[1:5, 1:3]
  save( phenoData, file = Rdata_file )
}else{
  load( Rdata_file )
}

## Category : PVT1 gene expression

picked_gene_data <- as.data.frame(t(raw_data[rownames(raw_data) == 'PVT1',]))
picked_gene_data[,2] <- rownames(picked_gene_data)
picked_gene_data <- picked_gene_data[order(picked_gene_data[,1]),]
PVT1_group <- c(rep('Low', floor(length(rownames(picked_gene_data))/2)),
                 rep('High', ceiling(length(rownames(picked_gene_data))/2)))
names(PVT1_group) <- factor(picked_gene_data[,2])

## OS
library(survival)  
os.model <- phenoData[, c('OS.time', 'OS')]

os.model <- os.model[rownames(os.model) %in% names(PVT1_group),]
os.model <- os.model[names(PVT1_group),]
os.model$type <- PVT1_group

fit <- survfit(Surv(OS.time, OS) ~ type, data = os.model)

library(survminer)
ggsurvplot(fit,
           pval = TRUE, conf.int = TRUE,
           linetype = 'type',
           surv.median.line = 'hv',
           ggtheme = theme_bw(),
           palette = c('#E7B800', '#2E9FDF')
)