TCGA

TCGA数据挖掘

2022-01-20  本文已影响0人  千容安

下载流程

数据较多时，用gdc-client，下载后配置环境变量


将之前manifest下载的txt文件找到：“gdc_manifest_20220119_144759.txt”

在命令行输入gdc-client download -m C:/Users/Administrator.DESKTOP-4UQ3Q0K/Desktop/TCGA/gdc_manifest_20220119_144759.txt
gdc-client download -m +路径+文件名
下载完成后在运行目录（C:\Users\Administrator.DESKTOP-4UQ3Q0K）里找到对应文件
可以提前创建输出文件夹方便整理

数据整理

在R中设置工作目录
setwd("C:\\Users\\Administrator.DESKTOP-4UQ3Q0K\\Desktop\\TCGA")
获取所有文件的地址
filepath<- dir(path = ".\\gdc_download_433",full.names = TRUE)
用这个循环将gdc_download_433中的所有文件复制到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码对应

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]
写入数据
write.csv(json_file_Info,file="..\\json_file_Info.csv")

数据融合

将数据读取融合成数据框。行名是Ensembl ID,列名是样品名称（TCGA的barcode）

删掉第一列（只有barcode一列了）
filesName_To_TCGA_BarcodeFile<- json_file_Info[-1]
读入样本文件夹中所有文件
countsFileNames<- dir(pattern = "counts$")
定义一个储存数据的数据框
allSampleRawCounts<-data.frame()
把读入的数据合并为数据框

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……点后面的数字是版本号，可以去掉

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

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

ID转换

因为像ENSG00000000003这样的ensembl_id不知道是什么基因，故进行转换
用clusterprofiler匹配的太少了，故用gtf文件进行转换
基因组注释文件（gtf）的下载方法：
网址：https://www.gencodegenes.org/human/release_39lift37.html
gencode下载的文件解压出错，用http://genome.ucsc.edu/cgi-bin/hgTables的话之后的Ensembl_ID_TO_Genename是空的，第三种方法可以成功:
网站：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数据框中

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

一个函数，通过gtf文件获取Ensemble_ID与基因名称的对应关系

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")

#融合数据
mergeRawCounts<-merge(Ensembl_ID_TO_Genename,RawCounts,by="Ensembl_ID")

#按照gene_id列进行排序,发现Ensembl_ID不同但gene ID有重复
mergeRawCounts <- mergeRawCounts[order(mergeRawCounts[,"gene_id"]),]
#根据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)
# 请求数据，查询 Illumina HiSeq 平台的数据
query <- GDCquery(project ="TCGA-BLCA",
                  data.category = "Transcriptome Profiling",
                  data.type ="Gene Expression Quantification" ,
                  workflow.type="HTSeq - Counts")

研究的是BLCA，研究其他方面就替换成其他肿瘤名称，使用getGDCprojects()$project_id得到各个癌种的项目id，总共有52个ID值。

使用TCGAbiolinks:::getProjectSummary(project)查看project中有哪些数据类型，如查询"TCGA-BLCA"，有8种数据类型，case_count为病人数，file_count为对应的文件数。要下载表达谱，可以设置data.category="Transcriptome Profiling"

# 从query中获取结果表，它可以选择带有cols参数的列，并使用rows参数返回若干行。
#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样本。 
#原始数据中有433个样本，对照有19个和实验组各414个。所以我们可以创建一个分组向量。
###################方法一 dger 
#包的安装
BiocManager::install("DESeq2")
BiocManager::install("edgeR")
library(DESeq2)
library(edgeR)
library("edgeR")
#创建分组
group <- c(rep(1,19),rep(2,414))

#创建DGEList类型变量
# 这一步相当于创建一列表。注意group中的顺序和counts中行名要对应，
#也就是对照组和实验组要指定正确。这里前19列为contro1，后414列为case。
y<-DGEList(counts=Counts,group=group)

#数据过滤（去除表达量低的数据）
keep <- filterByExpr(y)
y<-y[keep,,keep.lib.sizes=FALSE]
# 计算标准化因子
y<-calcNormFactors(y)
#计算离散度
y<-estimateDisp(y)
#显著性检验
et <-exactTest(y)
# 获取排名靠前的基因，这里设置n=100000是为了输出所有基因
et <- topTags(et,n=100000)
#转换为数据框类型
et <- as.data.frame(et)
# 将行名粘贴为数据框的第一列 
et<- cbind(rownames(et),et)
# 指定列名
colnames(et)<-c("gene_id","log2FoldChange","log2CPM","PValue","FDR")
# 保存数据到本地
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 体现上下调信息
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数据挖掘之基因表达差异分析_哔哩哔哩_bilibili

上一篇下一篇

猜你喜欢

热点阅读