加入我們用TCGAbiolinks下載之后我們可以看一下另外一種操作

image.png

data_download_from_gdc
打開后發(fā)現(xiàn)是1208個(gè)文件夾晾捏，為什么是1208個(gè)蒿涎，這跟我當(dāng)時(shí)學(xué)的數(shù)據(jù)有關(guān)系，可以不用管讓 每個(gè)文件夾里面還有個(gè)壓縮文件惦辛，我們現(xiàn)在的任務(wù)就是劳秋，如何把每個(gè)文件夾里面的壓縮文件放置在統(tǒng)一的文件夾中 假設(shè)我們的文件夾是這個(gè),
首先設(shè)置工作目錄到這一文件夾下：
00_data_read_in_one_file
并且他就是創(chuàng)建在 data_download_from_gdc中，所以現(xiàn)在文件夾的總體數(shù)目是應(yīng)該是1208加1就是1209個(gè) 我們在R語言中輸入

length(list.files())

這個(gè)命令中l(wèi)ist.files()的意思就是列舉當(dāng)前工作目錄中文件裙品，length的意思就是有有多少個(gè)俗批，一看，果真是1209個(gè) 其實(shí)

length(dir())

這個(gè)命令也能做這個(gè)事情市怎。 為什么我的文件夾要起一個(gè)這么詭異的名字岁忘？ 以00開頭？区匠？

這是因?yàn)槲倚枰屗幱谖募械谝粋€(gè)位置干像，我們來驗(yàn)證一下：

dir()[1]

發(fā)現(xiàn)第一個(gè)文件夾就是他，這樣我就每次循環(huán)訪問后面的1208個(gè)文件夾驰弄，每次都把看到的東西復(fù)制過來就可以了

你說我為什么這么事無巨細(xì)麻汰，因?yàn)槲铱梢院芎啙崳部梢院軉缕莞荩Q于我是否想要?jiǎng)e人聽的懂五鲫。

在R語言里面是可以直接創(chuàng)建文件夾的，鼠標(biāo)右擊也可以

dir.create("00_data_read_in_one_file") #創(chuàng)建新的文件夾,確保文件夾排在第一位

遍歷和復(fù)制,為什么從2開始岔擂，因?yàn)榈谝粋€(gè)文件夾你已經(jīng)知道了

for (dirname in dir()[2:length(dir())]){  
  file <- list.files(dirname,pattern = "*.counts")  #找到對應(yīng)文件夾中的內(nèi)容位喂，pattern可以是正則表達(dá)式
  file.copy(paste0(dirname,"/",file),"00_data_read_in_one_file")  #復(fù)制內(nèi)容到新的文件夾
}

運(yùn)行完了之后，我們可以打開文件夾看一下乱灵，確實(shí)在里面塑崖，這時(shí)候可以全部選擇，

將文件解壓到data_unzip文件夾,解壓數(shù)據(jù)1.42GB痛倚，

我們發(fā)現(xiàn)即使這個(gè)樣子规婆，文件的名稱也是怪怪的

那我們就來轉(zhuǎn)換，轉(zhuǎn)換的信息藏在metadata文件中蝉稳，這個(gè)要去看開始的那個(gè)帖子下載

切換一下工作目錄

setwd("~/skill_practice/BRCA_999")

注意一開始我就建立了 BRCA_999這個(gè)文件夾抒蚜， data_download_from_gdc是他的子文件夾， 00_data_read_in_one_file又是 data_download_from_gdc的子文件夾 metadata是json格式的

讀入json格式的文件,他是一個(gè)1208行耘戚，15列的數(shù)據(jù)框

metadata <- jsonlite::fromJSON("metadata.cart.2017-11-15T09_56_59.722935.json")

我們轉(zhuǎn)換的信息就是兩列filename和associatedentities削锰，我們把它選出來

require(dplyr)
metadata_id <- metadata %>% 
  dplyr::select(c(file_name,associated_entities)) 

我需要的是 file_name和樣本名稱，樣本名稱藏在了 associated_entities 列表中里面包括了 entity_id, case_id, entity_submitter_id, entity_type這四個(gè)項(xiàng)目毕莱，查看第一個(gè)了解一下

metadata$associatedentities[1] [[1]]
 entityid caseid 1 52033f64-1e6f-4657-a4fb-7cfeffc61951 39de7761-e762-4811-b95c-8216b79ae06b entitysubmitterid entitytype 1 TCGA-AN-A0XW-01A-11R-A109-07 aliquot

現(xiàn)在的想法是我把filename和 associated_entities中的 entity_submitter_id提取出來器贩，做成一個(gè)數(shù)據(jù)框颅夺，然后我批量對應(yīng)轉(zhuǎn)換

naid_df <- data.frame()
for (i in 1:1208){
  naid_df[i,1] <- substr(metadata_id$file_name[i],1,nchar(metadata_id$file_name[i])-3)
  naid_df[i,2] <- metadata_id$associated_entities[i][[1]]$entity_submitter_id
}

現(xiàn)在把1208個(gè)小文件讀入一個(gè)矩陣文件，并且給每一個(gè)文件加上 filename和 entity_submitter_id 論壇有道題目就是處理的這個(gè)問題

生信編程直播第四題：多個(gè)同樣的行列式文件合并起來

http://www.biotrainee.com:8080/thread-603-1-1.html

我自己也給出了R語言的解法

在R語言中將多個(gè)同樣的行列式文件合并起來

http://guoshipeng.com/2017/11/10/14-R-for-binding-colums/

但是當(dāng)時(shí)不知道蛹稍，TCGA的單個(gè)文件是沒有列名的吧黄，導(dǎo)致無法合并，所以本次要復(fù)雜一點(diǎn)

讀入所有解壓的文件 1208個(gè)

nameList <- list.files("data_unzip/")
location <- which(naid_df==nameList[1],arr.ind = TRUE) ##which函數(shù)有一個(gè)已知value返回坐標(biāo)的功能
TCGA_id <- as.character(naid_df[location[1],2]) ##通過坐標(biāo)唆姐，獲取TCGA_id
expr_df<- read.table(paste0("data_unzip/",nameList[1]),stringsAsFactors = F, header = F) #讀入第一個(gè)文件拗慨，保存為data.frame
names(expr_df) <- c("gene_id",TCGA_id) #給剛才數(shù)據(jù)庫命名

這邊開始批量作業(yè)

for (i in 2:length(nameList)){
  location <- which(naid_df==nameList[i],arr.ind = TRUE)
  TCGA_id <- as.character(naid_df[location[1],2])
  dfnew <- read.table(paste0("data_unzip/",nameList[i]),stringsAsFactors = F,header = F)
  names(dfnew) <- c("gene_id",TCGA_id)
  expr_df <- inner_join(expr_df,dfnew,by="gene_id")
}

晚上走的時(shí)候沒運(yùn)行完，早上來的時(shí)候已經(jīng)完畢奉芦，限速環(huán)節(jié)應(yīng)該是read.table赵抢，早上再來嘗試運(yùn)行一次總是說內(nèi)存不夠 我嘗試了一下fread來解決這個(gè)問題：

require(data.table)
nameList <- list.files("data_unzip/")
location <- which(naid_df==nameList[1],arr.ind = TRUE)
TCGA_id <- as.character(naid_df[location[1],2])
expr_df<- fread(paste0("data_unzip/",nameList[1]))
names(expr_df) <- c("gene_id",TCGA_id)
for (i in 2:length(nameList)){
  location <- which(naid_df==nameList[i],arr.ind = TRUE)
  TCGA_id <- as.character(naid_df[location[1],2])
  dfnew <- fread(paste0("data_unzip/",nameList[i]))
  names(dfnew) <- c("gene_id",TCGA_id)
  expr_df <- inner_join(expr_df,dfnew,by="gene_id")
}

結(jié)果大概2分鐘搞定，速度喜人ＩΑ７橙础！先巴！ 總共60488行,查看最后幾行發(fā)現(xiàn)有5行不是我們要的

tail(expr_df$gene_id,10)

去掉最后五行

expr_df <- expr_df[1:(length(expr_df$gene_id)-5),]

保存數(shù)據(jù)其爵，大概是3個(gè)G左右

save(expr_df,file = "expr_df.Rda")

下面開始id轉(zhuǎn)換，信息在GTF文件中表達(dá)矩陣?yán)锩娴膅ene_id有小數(shù)點(diǎn)伸蚯， 而GTF文件中沒有,調(diào)整一下摩渺，先以“.”分列，在去掉小數(shù)點(diǎn)后的列

require(dplyr)
require(tidyr)
expr_df_nopoint <- expr_df %>% 
  tidyr::separate(gene_id,into = c("gene_id","drop"),sep="\\.") %>% 
  dplyr::select(-drop)
save(expr_df_nopoint,file = "expr_df_nopoint.Rda")
load(file = "expr_df_nopoint.Rda")

下載GTF文件來注釋 ftp://ftp.ensembl.org/pub/release-90/gtf/homo_sapiens

安裝包：

source("https://bioconductor.org/biocLite.R")
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
biocLite("rtracklayer")
biocLite("SummarizedExperiment")

讀入GTF數(shù)據(jù)

gtf1 <- rtracklayer::import('Homo_sapiens.GRCh38.90.chr.gtf')
gtf_df <- as.data.frame(gtf1)

保存數(shù)據(jù)

save(gtf_df,file = "gtf_df.Rda")

讀入27個(gè)變量剂邮，2612129個(gè)觀測摇幻，測試一下顯示的不錯(cuò)

test <- gtf_df[1:5,]
View(test)

進(jìn)展的很快，我們現(xiàn)在可以提取mRNA的表達(dá)矩陣?yán)玻?以gtf文件中的 gene_biotype為標(biāo)準(zhǔn)挥萌，里面寫 protein_coding的就是編碼基因

首先要把這些基因提取出來绰姻，然后與表達(dá)譜融合，我在這個(gè)例子還提取了 gene_name, gene_id瑞眼，所以最后的時(shí)候，我把三種表達(dá)方式合在了一起 這樣棵逊，以后我無論用什么方式都可以選出我要的基因了

require(dplyr)
require(tidyr)
mRNA_exprSet <- gtf_df %>% 
  dplyr::filter(type=="gene",gene_biotype=="protein_coding") %>% #篩選gene,和編碼指標(biāo)
  dplyr::select(c(gene_name,gene_id,gene_biotype)) %>% 
  dplyr::inner_join(expr_df_nopoint,by ="gene_id") %>% 
  tidyr::unite(gene_id,gene_name,gene_id,gene_biotype,sep = " | ")

得到19688行伤疙，跟我們的認(rèn)知很吻合

保存數(shù)據(jù)

save(mRNA_exprSet,file = "mRNA_exprSet.Rda")

下面進(jìn)行標(biāo)準(zhǔn)化

suppressPackageStartupMessages(library(DESeq2))
suppressPackageStartupMessages(library(edgeR))

本次使用的方法是DESeq2中vst法 至于有哪幾種方法，還有該選什么方法辆影，明天會(huì)有一個(gè)帖子介紹

制作metadata徒像，區(qū)別腫瘤和正常組，這時(shí)候需要對TCGA的id有個(gè)了解蛙讥，從左往右數(shù)數(shù)锯蛀，第14,15上的數(shù)字很重要 其中01-09是tumor，也就是癌癥樣本次慢；其中10-29是normal旁涤，也就是癌旁

metadata <- data.frame(names(mRNA_exprSet)[-1])
for (i in 1:length(metadata[,1])) {
  num <- as.numeric(substring(metadata[i,1],14,15))
  if (num %in% seq(1,9)) {metadata[i,2] <- "T"}
  if (num %in% seq(10,29)) {metadata[i,2] <- "N"}
}

加個(gè)名稱

names(metadata) <- c("TCGA_id","sample")

將sample轉(zhuǎn)化成factor

metadata$sample <- as.factor(metadata$sample)

我們可總結(jié)一下,有112例normal翔曲，1096例腫瘤組織

metadata %>% dplyr::group_by(sample) %>% summarise(n())

保存數(shù)據(jù)，我想問一下劈愚，為什么你老是要保存數(shù)據(jù)瞳遍，因?yàn)槲已劾锍：瑴I花

save(metadata,file = "metadata.Rda")

制作countData

mycounts <- mRNA_exprSet
keepGene=rowSums(edgeR::cpm(mycounts[-1])>0) >=2
table(keepGene);dim(mycounts)

這是數(shù)據(jù)維度19668 和1209

dim(mycounts[keepGene,])

剩余19546 1209

mycounts <-mycounts[keepGene,]
構(gòu)建 DESeq對象

dds <-DataSetFromMatrix(countData=mycounts, 
                              colData=metadata, 
                              design=~sample,
                              tidy=TRUE)

標(biāo)準(zhǔn)化vst法，5分鐘

vsd <- vst(dds, blind = FALSE)
mRNA_exprSet_vst <- as.data.frame(assay(vsd))
save(mRNA_exprSet_vst,file = "mRNA_exprSet_vst.Rda")

load(file = "dds.Rda")
下面就是完整的過程

  dds <- DESeq(dds)
  save(dds,file = "dds_DEseq.Rda")
  res <- results(dds, tidy=TRUE) #獲取結(jié)果
  res <- as_tibble(res)
  require(dplyr)
  res <- res %>%
    separate(row,c("symbol","ensemble","genetype"),sep = " \\| ") %>%
    dplyr::select(- c(ensemble,genetype)) %>%
    arrange(desc(abs(log2FoldChange))) %>% #排序菌羽。為了去重
    distinct(symbol,.keep_all = TRUE) %>%
    arrange(desc(log2FoldChange))#再次按照log2FoldChange從大到小排序
  save(res,file = "result.Rda")

既然做了差異分析掠械，那對下游基因注釋就是順手 的事情啦，我們沒有選擇注祖，就用Y叔的包猾蒂，clusterprofiler 在這之前我們需要獲得差異基因列表 這一步，就是基本數(shù)據(jù)框的基本操作是晨，按照pvalue和fc來排序肚菠，選擇自己一定數(shù)量的基因，數(shù)量你來定署鸡，最終得到基因列表gene(我沒有演示案糙，需要自己做)

加載包，加載之前確定自己裝了

library(DOSE)
library(GO.db)
library(org.Hs.eg.db)
library(topGO)
library(GSEABase)
library(clusterProfiler)

基因名稱轉(zhuǎn)換靴庆，返回的是數(shù)據(jù)框

gene = bitr(gene, fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db")
head(gene)

GO分析 細(xì)胞組分

ego_CC <- enrichGO(gene = gene$ENTREZID,
                  OrgDb= org.Hs.eg.db,
                  ont = "CC",
                  pAdjustMethod = "BH",
                  minGSSize = 1,
                  pvalueCutoff = 0.01,
                  qvalueCutoff = 0.01,
                  readable = TRUE)
生物過程

ego_BP <- enrichGO(gene = gene$ENTREZID,
                  OrgDb= org.Hs.eg.db,
                  ont = "BP",
                  pAdjustMethod = "BH",
                  minGSSize = 1,
                  pvalueCutoff = 0.01,
                  qvalueCutoff = 0.01,
                  readable = TRUE)
分子功能：

ego_MF <- enrichGO(gene = gene$ENTREZID,
                  OrgDb= org.Hs.eg.db,
                  ont = "MF",
                  pAdjustMethod = "BH",
                  minGSSize = 1,
                  pvalueCutoff = 0.01,
                  qvalueCutoff = 0.01,
                  readable = TRUE)
作圖

barplot(ego_CC, showCategory=20,title="EnrichmentGO_CC")#條狀圖时捌，按p從小到大排的
dotplot(ego_BP,title="EnrichmentGO_BP_dot")#點(diǎn)圖，按富集的數(shù)從大到小的
KEGG分析

kk <- enrichKEGG(gene = gene$ENTREZID,
                organism ="human",
                pvalueCutoff = 0.01,
                qvalueCutoff = 0.01,
                minGSSize = 1,
                #readable = TRUE ,
                use_internal_data =FALSE)
作圖

barplot(kk)
dotplot(kk)