下載流程

數(shù)據(jù)較多時，用gdc-client，下載后配置環(huán)境變量

數(shù)據(jù)整理

在R中設(shè)置工作目錄
setwd("C:\\Users\\Administrator.DESKTOP-4UQ3Q0K\\Desktop\\TCGA")
獲取所有文件的地址
filepath<- dir(path = ".\\gdc_download_433",full.names = TRUE)
用這個循環(huán)將gdc_download_433中的所有文件復(fù)制到SampleFiles中

for(wd in filepath){
  files<-dir(path=wd,pattern = "gz$")
  fromfilepath<- paste(wd,"\\",files,sep="")
  tofilepath<-paste(".\\SampleFiles\\",files,sep="")
  file.copy(fromfilepath,tofilepath)
}

更改工作目錄到樣本文件夾
setwd("C:\\Users\\Administrator.DESKTOP-4UQ3Q0K\\Desktop\\TCGA\\SampleFiles")
查看滿足要求的文件

countsFiles<- dir(path = ".\\",pattern="gz$")
length(countsFiles)
# [1] 433

全部解壓并刪除原有文件

library(R.utils)
sapply(countsFiles,gunzip)

處理json文件，加兩個.代表工作路徑的上一級目錄

library(rjson)
metadata_json_File<-fromJSON(file="..\\metadata.cart.2022-01-19.json")  

將文件名稱和TCGA的Barcode碼對應(yīng)

json_file_Info<- data.frame(filesName=c(),TCGA_Barcode=c())
for(i in 1:length(metadata_json_File)){
  TCGA_Barcode<-metadata_json_File[[i]][["associated_entities"]][[1]][["entity_submitter_id"]]
  file_name<-metadata_json_File[[i]][["file_name"]]
  json_file_Info<-rbind(json_file_Info,data.frame(filesName=file_name,TCGA_Barcode=TCGA_Barcode))
}

將文件名作為行名
rownames(json_file_Info)<-json_file_Info[,1]
寫入數(shù)據(jù)
write.csv(json_file_Info,file="..\\json_file_Info.csv")

數(shù)據(jù)融合

將數(shù)據(jù)讀取融合成數(shù)據(jù)框。行名是Ensembl ID,列名是樣品名稱（TCGA的barcode）

刪掉第一列（只有barcode一列了）
filesName_To_TCGA_BarcodeFile<- json_file_Info[-1]
讀入樣本文件夾中所有文件
countsFileNames<- dir(pattern = "counts$")
定義一個儲存數(shù)據(jù)的數(shù)據(jù)框
allSampleRawCounts<-data.frame()
把讀入的數(shù)據(jù)合并為數(shù)據(jù)框

for(txtFile in countsFileNames){
      Samplecounts <- read.table(txtFile,header = FALSE)
      rownames(Samplecounts) <- Samplecounts[,1]
      Samplecounts <- Samplecounts[-1]
      colnames(Samplecounts) <- filesName_To_TCGA_BarcodeFile[paste(txtFile,".gz",sep = ""),]
      if (dim(allSampleRawCounts)[1]==0){
          allSampleRawCounts <- Samplecounts
        }
      else
        {allSampleRawCounts<- cbind(allSampleRawCounts,Samplecounts)}
    }

寫入本地
write.csv(allSampleRawCounts,file = "..\\allSampleRawCounts.csv")
生成的allSampleRawCounts里面，ENS……點(diǎn)后面的數(shù)字是版本號，可以去掉

ensembl_id <- substr(row.names(allSampleRawCounts),1,15)
rownames(allSampleRawCounts) <- ensembl_id
write.csv(allSampleRawCounts,file = "..\\RawCounts.csv")

RawCounts.csv和allSampleRawCounts.csv不同在于去掉了版本號

ID轉(zhuǎn)換

因?yàn)橄馝NSG00000000003這樣的ensembl_id不知道是什么基因，故進(jìn)行轉(zhuǎn)換
用clusterprofiler匹配的太少了，故用gtf文件進(jìn)行轉(zhuǎn)換
基因組注釋文件（gtf）的下載方法：
網(wǎng)址：https://www.gencodegenes.org/human/release_39lift37.html
gencode下載的文件解壓出錯，用http://genome.ucsc.edu/cgi-bin/hgTables的話之后的Ensembl_ID_TO_Genename是空的，第三種方法可以成功:
網(wǎng)站：ftp://ftp.ensembl.org/pub/current_gtf/homo_sapiens/
下載Homo_sapiens.GRCh38.105.gtf.gz，解壓為Homo_sapiens.GRCh38.105.gtf
要把gtf文件放在TCAG這個文件夾

添加一列Ensemble_ID到RawCounts數(shù)據(jù)框中

RawCounts<-allSampleRawCounts
Ensembl_ID<- data.frame(Ensembl_ID=row.names(RawCounts))
rownames(Ensembl_ID)<-Ensembl_ID[,1]
RawCounts<-cbind(Ensembl_ID,RawCounts)

一個函數(shù)，通過gtf文件獲取Ensemble_ID與基因名稱的對應(yīng)關(guān)系

get_map=function(input){
  if(is.character(input)){
    if(!file.exists(input)) stop("Bad input file.")
    message('Treat input as file')
    input=data.table::fread(input,header=FALSE)
  }else{
    data.table::setDT(input)
  }
  input=input[input[[3]]=='gene',]
  pattern_id=".*gene_id \"([^;]+)\";.*"
  pattern_name=".*gene_name \"([^;]+)\";.*"
  gene_id=sub(pattern_id,"\\1",input[[9]])
  gene_name=sub(pattern_name,"\\1",input[[9]])
  Ensembl_ID_TO_Genename<-data.frame(gene_id=gene_id,
                                     gene_name=gene_name,
                                     stringsAsFactors = FALSE)
  return(Ensembl_ID_TO_Genename)
}

如果pattern_id=".*gene_id之后沒有空格的話，得到的Ensembl_ID_TO_Genename就不會是干凈的id和name的注釋信息，找出這個空格的問題用了兩天

Ensembl_ID_TO_Genename<-get_map("..\\gencode.v39.annotation.gtf.gz")

gtf_Ensembl_ID<-substr(Ensembl_ID_TO_Genename[,1],1,15)
Ensembl_ID_TO_Genename<-data.frame(gtf_Ensembl_ID,Ensembl_ID_TO_Genename[,2])
colnames(Ensembl_ID_TO_Genename)<-c("Ensembl_ID","gene_id")

write.csv(Ensembl_ID_TO_Genename,file = "..\\Ensembl_ID_TO_Genename.csv")

#融合數(shù)據(jù)
mergeRawCounts<-merge(Ensembl_ID_TO_Genename,RawCounts,by="Ensembl_ID")

#按照gene_id列進(jìn)行排序,發(fā)現(xiàn)Ensembl_ID不同但gene ID有重復(fù)
mergeRawCounts <- mergeRawCounts[order(mergeRawCounts[,"gene_id"]),]
#根據(jù)gene_id列建立索引
index<-duplicated(mergeRawCounts$gene_id)
#我們想要的那一行為FALSE，所以要取反。
mergeRawCounts <-mergeRawCounts[!index,]
#利用基因名稱作為行名
rownames (mergeRawCounts)<-mergeRawCounts[,"gene_id"]
# 刪除前2列
BLCA_Counts_expMatrix<-mergeRawCounts[,-c(1:2)]
#保存文件
write.csv(BLCA_Counts_expMatrix,file ="..\\BLCA_Counts_expMatrix.csv")

差異分析

counts<- read.csv("..\\BLCA_Counts_expMatrix.csv",header = T,row.names = 1)


BiocManager::install("TCGAbiolinks")
library(TCGAbiolinks)
# 請求數(shù)據(jù)，查詢 Illumina HiSeq 平臺的數(shù)據(jù)
query <- GDCquery(project ="TCGA-BLCA",
                  data.category = "Transcriptome Profiling",
                  data.type ="Gene Expression Quantification" ,
                  workflow.type="HTSeq - Counts")

研究的是BLCA，研究其他方面就替換成其他腫瘤名稱，使用getGDCprojects()$project_id得到各個癌種的項(xiàng)目id，總共有52個ID值。

使用TCGAbiolinks:::getProjectSummary(project)查看project中有哪些數(shù)據(jù)類型，如查詢"TCGA-BLCA"，有8種數(shù)據(jù)類型，case_count為病人數(shù)，file_count為對應(yīng)的文件數(shù)。要下載表達(dá)譜，可以設(shè)置data.category="Transcriptome Profiling"

# 從query中獲取結(jié)果表，它可以選擇帶有cols參數(shù)的列，并使用rows參數(shù)返回若干行。
#433個barcode
samplesDown <-getResults(query,cols=c("cases"))
#篩選
#414個腫瘤樣本的barcode
dataSmTP <-TCGAquery_SampleTypes(barcode = samplesDown,typesample ="TP")
#19個正常組織的barcode
dataSmNT<-TCGAquery_SampleTypes(barcode=samplesDown,typesample ="NT")

#重新排序樣本順序，正常組織樣本在前，腫瘤樣本在后，即前19列為正常樣本
Counts<-data.frame(c(BLCA_Counts_expMatrix[,dataSmNT],BLCA_Counts_expMatrix[,dataSmTP]))
rownames(Counts)<-row.names(BLCA_Counts_expMatrix)
colnames(Counts) <-c(dataSmNT,dataSmTP)

# 在edger中，1代表contro1樣本，2代表case樣本。 
#原始數(shù)據(jù)中有433個樣本，對照有19個和實(shí)驗(yàn)組各414個。所以我們可以創(chuàng)建一個分組向量。
###################方法一 dger 
#包的安裝
BiocManager::install("DESeq2")
BiocManager::install("edgeR")
library(DESeq2)
library(edgeR)
library("edgeR")
#創(chuàng)建分組
group <- c(rep(1,19),rep(2,414))

#創(chuàng)建DGEList類型變量
# 這一步相當(dāng)于創(chuàng)建一列表。注意group中的順序和counts中行名要對應(yīng)，
#也就是對照組和實(shí)驗(yàn)組要指定正確。這里前19列為contro1，后414列為case。
y<-DGEList(counts=Counts,group=group)

#數(shù)據(jù)過濾（去除表達(dá)量低的數(shù)據(jù)）
keep <- filterByExpr(y)
y<-y[keep,,keep.lib.sizes=FALSE]
# 計算標(biāo)準(zhǔn)化因子
y<-calcNormFactors(y)
#計算離散度
y<-estimateDisp(y)
#顯著性檢驗(yàn)
et <-exactTest(y)
# 獲取排名靠前的基因，這里設(shè)置n=100000是為了輸出所有基因
et <- topTags(et,n=100000)
#轉(zhuǎn)換為數(shù)據(jù)框類型
et <- as.data.frame(et)
# 將行名粘貼為數(shù)據(jù)框的第一列 
et<- cbind(rownames(et),et)
# 指定列名
colnames(et)<-c("gene_id","log2FoldChange","log2CPM","PValue","FDR")
# 保存數(shù)據(jù)到本地
write.table(et,"..\\all_BLCA_DEG.xls",sep ="\t",col.names=TRUE, 
            row.names = FALSE, quote =FALSE, na="")
#差異基因篩選
etSig <-et[which(et$PValue<0.05&abs(et$log2FoldChange)>1),]
# 加入一列，up_down 體現(xiàn)上下調(diào)信息
etSig[which(etSig$log2Foldchange>0),"up_down"]<-"Up"
etSig[which(etSig$log2Foldchange<0),"up_down"]<-"Down" 
# 保存文件 
write.table(etSig,"..\\BLCA_DEG.x1s",sep="t",col.names =TRUE,row.names =FALSE,quote =FALSE,na="")

主要是在基因注釋那里花了很長時間，也可以使用DESeq2包來分析差異基因。詳情TCGA數(shù)據(jù)挖掘之基因表達(dá)差異分析_嗶哩嗶哩_bilibili