下载流程
数据较多时，用gdc-client，下载后配置环境变量
将之前manifest下载的txt文件找到：“gdc_manifest_20220119_144759.txt”
在命令行输入
gdc-client download -m C:/Users/Administrator.DESKTOP-4UQ3Q0K/Desktop/TCGA/gdc_manifest_20220119_144759.txtgdc-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
网友评论