1. 整理表达矩阵
下载的文件是按样本存放的，每个tsv文件中都记录着一个样本的基因表达量，需要将所有tsv文件合并，得到所有样本的基因表达量的表格。
转录组数据合并
rm(list = ls())
library(stringr)
library(jsonlite)
library(progress)
options(stringsAsFactors = F)
cancer_type <- "TCGA-COAD" # 设置肿瘤类型（需修改）
setwd("D:/BioInformationAnalyze/TCGAdata/CRC") # 设置工作目录（需修改）
## 获取每个样品id和相应tsv文件名的对应关系：
json <- jsonlite::fromJSON("metadata.cart.2022-12-24.json") # 读取Metadata文件中每个样品的信息（需修改matadata文件的名称）
View(json)
sample_id <- sapply(json$associated_entities, function(x){x[,1]}) # 提取所有样品id，如：TCGA-W5-AA2X-01A-11R-A41I-07
file_sample <- data.frame(sample_id, file_name = json$file_name) # 获取样品id和相应tsv文件名的对应关系
## 合并所有样品id的基因表达量：
count_file <- list.files("expdata", pattern = "*.tsv$", recursive = TRUE) # 获取expdata文件夹下的所有tsv表达文件的“路径+文件名”
count_file_name <- strsplit(count_file, split="/") %>%
sapply(function(x){x[2]}) # 在count_file中分割出文件名
## 新建空白矩阵，行数与tsv文件中基因名行数形同（需查看并修改）
matrix = data.frame(matrix(nrow = 60660, ncol = 0))
## 取每个文件中的一列数据（count或tpm），然后合并
pb <- progress_bar$new(format = "(:spin) [:bar] :percent 已用时::elapsedfull 预计剩余时间::eta",
total = length(count_file), clear = FALSE, width = 80) # 设置进度条
for (i in 1:length(count_file)){
data<- paste0("expdata/", count_file[i]) %>%
read.delim(fill = TRUE, header = FALSE, row.names = 1)
data <- data[-c(1:6),]
data <- data[3] # [3]是count，[6]是tpm，[7]是fpkm。一般用tpm，但差异分析R包需要用count。生存分析用tpm。保存3和6两套数据。
colnames(data) <- file_sample$sample_id[which(file_sample$file_name == count_file_name[i])]
matrix <- cbind(matrix, data)
pb$tick() # 显示进度条
Sys.sleep(1 / 100) # 用于计算时间
} # 得到了表达矩阵 matrix ：行名为样品id，列名为Ensembl ID，值为基因的表达量
表达矩阵 matrix 长这样：
转录组数据拆分，mRNA，lncRNA
表达矩阵是很多种基因类型的混合，我们需要将分析常用的mRNA（protein_coding）、lncRNA分别挑出来。
虽然表达矩阵中也有miRNA，但其gene_name不是常规命名。若要分析mirna，需在TCGA-GDC网站另外下载，其数量与表达矩阵中是一致的。
data <- paste0("expdata/", count_file[1]) %>%
read.delim(fill = TRUE, header = FALSE, row.names = 1) %>%
as.matrix() # 随便读取一个tsv文件
gene_name <- data[-c(1:6),1] # 挑出gene_name
gene_type <- data[-c(1:6),2] # 挑出gene_type
table(gene_type)
matrix <- cbind(gene_name, gene_type, matrix) # 将 gene_name 和 gene_type 添加到表达矩阵 mateix 的最前面
mRNA_matrix <- matrix[matrix$gene_type == "protein_coding",] ; dim(mRNA_matrix) # 拆分得到mRNA的表达矩阵 mRNA_matrix
lnc_matrix <- matrix[matrix$gene_type == "lncRNA",] ; dim(lnc_matrix) # 拆分得到lncRNA的表达矩阵 lncRNA_matrixm
mRNA_matrix 和 lncRNA_matrix 长这样：
设置Gene Symbol为表达矩阵的列名（上一步得到的列名是Ensembl ID）
## 按gene_name列去除重复的基因，保留每个基因最大表达量结果 matrix0 <- aggregate( . ~ gene_name, data = lnc_matrix, max) # 需要lnc就写lnc，需要mRNA就写mRNA。这里以lnc为例。这一步运行时间较长。 # 函数解析：aggregate函数中，“. ~ gene_name”表示以gene_name为分类，分别计算其他列的最大值。 rownames(matrix0) <- matrix0[,1] # 将“gene_name”列设为列名 matrix0 <- matrix0[, -c(1,2)] matrix0 <- as.data.frame(lapply(matrix0, as.numeric), row.names = rownames(matrix0), check.names = F)
matrix0 长这样：
过滤掉在很多样本种表达量都为0基因
表达矩阵整理完成，需要过滤掉在很多样本中表达量都为0的基因。过滤标准不唯一。
判断“matrix0”中的每个基因，是否在至少44个样本中>1。保留符合条件的基因，剔除不符合条件的基因。样本数“>44”根据总样本量修改，可以是数值，比如大于正常样本数“>=40”，也可以是百分比，如大于总样本数的75%“>=0.75*(ncol(matrix0))”
exp <- matrix0[apply(matrix0, 1, function(x)sum(x >= 1) >= 44), ] # counts用这个。x = 1表示对 matrix0 的每一行都执行 function 函数 # exp = matrix0[apply(matrix0, 1, function(x)sum(x > 0.1) >= 44), ] # tpm用这个。
过滤后的表达矩阵 exp 长这样：
定义表达矩阵分组信息
group_list <- ifelse(as.numeric(str_sub(colnames(matrix0), 14, 15)) < 10, "tumor", "normal") %>%
factor(levels = c("normal","tumor")) # 将正常样品和肿瘤样品分开，并设置因子，便于后续分析
table(group_list)
保存转录组数据的整理结果
save(exp, group_list, cancer_type, file = paste0(cancer_type, "_01.1.lncRNA_Expdata_count.Rdata")) # save(exp, group_list, cancer_type, file = paste0(cancer_type, "_01.2.lncRNA_Expdata_tpm.Rdata")) # 换成 tpm 再运行、保存一次
2. 整理临床信息
rm(list = ls())
library(stringr)
library(jsonlite)
library(XML)
options(stringsAsFactors = F)
cancer_type <- "TCGA-COAD" # 设置肿瘤类型（需修改）
setwd("D:/BioInformationAnalyze/TCGAdata/CRC") # 设置工作目录（需修改）
先对单个文件进行检查（备用）
确认文件无误时可不检查
result <- xmlParse( "./clinical/00b9eb65-04c6-4d62-9cf6-6b77d75ab79b/nationwidechildrens.org_clinical.TCGA-DM-A285.xml") # 读取单个xml文件，需修改文件目录和文件名 rootnode <- xmlRoot(result) # 将根节点从xml文件中移除 rootsize <- xmlSize(rootnode) # 查找根中的节点数 print(rootnode[1]) # 查看第一个节点 print(rootnode[2]) # 查看第二个节点，可以确认第二个节点为临床信息的保存位置 xmldataframe <- xmlToDataFrame(rootnode[2]) # 将xml中第二个节点转换为数据框 head(t(xmlToDataFrame(rootnode[2]))) # 转置并查看该数据框
读取所有临床数据文件
## 读取所有临床数据文件路径
xmls <- dir("clinical/", pattern = "*.xml$", recursive = T) # 读取目录为“clinical/”，读取模式为以“*.xml”结尾的文件，recursive表示列表是否应重新进入目录
## 构建读取函数
td = function(x){
result <- xmlParse(file.path("clinical/", x))
rootnode <- xmlRoot(result)
xmldataframe <- xmlToDataFrame(rootnode[2])
return(t(xmldataframe))
}
## 读取并整合所有临床信息
cl = lapply(xmls, td) # 对于xmls中的每一个文件，执行td函数，输出每一个病人的临床信息，并整合到列表中
cl_df <- t(do.call(cbind,cl)) # 合并列表中的每一个元素为数据框，并将数据框转置（转置的输出结果为矩阵）
cl_df[1:3,1:3]
clinical = data.frame(cl_df)# 将转置的结果重新转变为数据框
clinical[1:4,1:4]
保存临床数据的整理结果
save(clinical, cancer_type, file = paste0(cancer_type,"_02.lncRNA_Clinical.Rdata"))
临床数据整理补充方法：
用perl整理数据（补充方法二，推荐） ----
## 将”getClinical.pl“文件复制到clinical文件夹下
## 在clinical文件夹下运行terminal，运行命令“perl getClinical.pl”。
## 即可在clinical文件夹下得到“clinical.xls”文件
## 其中包含：Id futime fustat Age Gender Grade Stage T M N
## 以下为getClinical.pl的perl脚本，复制到txt文件中重命名为.pl即可
use strict;
#use warnings;
use XML::Simple;
opendir(RD, ".") or die $!;
my @dirs=readdir(RD);
closedir(RD);
open(WF,">clinical.xls") or die $!;
print WF "Id\tfutime\tfustat\tAge\tGender\tGrade\tStage\tT\tM\tN\n";
foreach my $dir(@dirs){
#print $dir . "\n";
next if($dir eq '.');
next if($dir eq '..');
#print $dir . "\n";
if(-d $dir){
opendir(RD,"$dir") or die $!;
while(my $xmlfile=readdir(RD)){
if($xmlfile=~/\.xml$/){
#print "$dir\\$xmlfile\n";
my $userxs = XML::Simple->new(KeyAttr => "name");
my $userxml="";
if(-f "$dir/$xmlfile"){
$userxml = $userxs->XMLin("$dir/$xmlfile");
}else{
$userxml = $userxs->XMLin("$dir\$xmlfile");
}
# print output
#open(WF,">dumper.txt") or die $!;
#print WF Dumper($userxml);
#close(WF);
my $disease_code=$userxml->{'admin:admin'}{'admin:disease_code'}{'content'}; #get disease code
my $disease_code_lc=lc($disease_code);
my $patient_key=$disease_code_lc . ':patient'; #ucec:patient
my $follow_key=$disease_code_lc . ':follow_ups';
my $patient_barcode=$userxml->{$patient_key}{'shared:bcr_patient_barcode'}{'content'}; #TCGA-AX-A1CJ
my $gender=$userxml->{$patient_key}{'shared:gender'}{'content'}; #male/female
my $age=$userxml->{$patient_key}{'clin_shared:age_at_initial_pathologic_diagnosis'}{'content'};
my $race=$userxml->{$patient_key}{'clin_shared:race_list'}{'clin_shared:race'}{'content'}; #white/black
my $grade=$userxml->{$patient_key}{'shared:neoplasm_histologic_grade'}{'content'}; #G1/G2/G3
my $clinical_stage=$userxml->{$patient_key}{'shared_stage:stage_event'}{'shared_stage:clinical_stage'}{'content'}; #stage I
my $clinical_T=$userxml->{$patient_key}{'shared_stage:stage_event'}{'shared_stage:tnm_categories'}{'shared_stage:clinical_categories'}{'shared_stage:clinical_T'}{'content'};
my $clinical_M=$userxml->{$patient_key}{'shared_stage:stage_event'}{'shared_stage:tnm_categories'}{'shared_stage:clinical_categories'}{'shared_stage:clinical_M'}{'content'};
my $clinical_N=$userxml->{$patient_key}{'shared_stage:stage_event'}{'shared_stage:tnm_categories'}{'shared_stage:clinical_categories'}{'shared_stage:clinical_N'}{'content'};
my $pathologic_stage=$userxml->{$patient_key}{'shared_stage:stage_event'}{'shared_stage:pathologic_stage'}{'content'}; #stage I
my $pathologic_T=$userxml->{$patient_key}{'shared_stage:stage_event'}{'shared_stage:tnm_categories'}{'shared_stage:pathologic_categories'}{'shared_stage:pathologic_T'}{'content'};
my $pathologic_M=$userxml->{$patient_key}{'shared_stage:stage_event'}{'shared_stage:tnm_categories'}{'shared_stage:pathologic_categories'}{'shared_stage:pathologic_M'}{'content'};
my $pathologic_N=$userxml->{$patient_key}{'shared_stage:stage_event'}{'shared_stage:tnm_categories'}{'shared_stage:pathologic_categories'}{'shared_stage:pathologic_N'}{'content'};
$gender=(defined $gender)?$gender:"unknow";
$age=(defined $age)?$age:"unknow";
$race=(defined $race)?$race:"unknow";
$grade=(defined $grade)?$grade:"unknow";
$clinical_stage=(defined $clinical_stage)?$clinical_stage:"unknow";
$clinical_T=(defined $clinical_T)?$clinical_T:"unknow";
$clinical_M=(defined $clinical_M)?$clinical_M:"unknow";
$clinical_N=(defined $clinical_N)?$clinical_N:"unknow";
$pathologic_stage=(defined $pathologic_stage)?$pathologic_stage:"unknow";
$pathologic_T=(defined $pathologic_T)?$pathologic_T:"unknow";
$pathologic_M=(defined $pathologic_M)?$pathologic_M:"unknow";
$pathologic_N=(defined $pathologic_N)?$pathologic_N:"unknow";
my $survivalTime="";
my $vital_status=$userxml->{$patient_key}{'clin_shared:vital_status'}{'content'};
my $followup=$userxml->{$patient_key}{'clin_shared:days_to_last_followup'}{'content'};
my $death=$userxml->{$patient_key}{'clin_shared:days_to_death'}{'content'};
if($vital_status eq 'Alive'){
$survivalTime="$followup\t0";
}
else{
$survivalTime="$death\t1";
}
for my $i(keys %{$userxml->{$patient_key}{$follow_key}}){
eval{
$followup=$userxml->{$patient_key}{$follow_key}{$i}{'clin_shared:days_to_last_followup'}{'content'};
$vital_status=$userxml->{$patient_key}{$follow_key}{$i}{'clin_shared:vital_status'}{'content'};
$death=$userxml->{$patient_key}{$follow_key}{$i}{'clin_shared:days_to_death'}{'content'};
};
if($@){
for my $j(0..5){ #假设最多有6次随访
my $followup_for=$userxml->{$patient_key}{$follow_key}{$i}[$j]{'clin_shared:days_to_last_followup'}{'content'};
my $vital_status_for=$userxml->{$patient_key}{$follow_key}{$i}[$j]{'clin_shared:vital_status'}{'content'};
my $death_for=$userxml->{$patient_key}{$follow_key}{$i}[$j]{'clin_shared:days_to_death'}{'content'};
if( ($followup_for =~ /\d+/) || ($death_for =~ /\d+/) ){
$followup=$followup_for;
$vital_status=$vital_status_for;
$death=$death_for;
my @survivalArr=split(/\t/,$survivalTime);
if($vital_status eq 'Alive'){
if($followup>$survivalArr[0]){
$survivalTime="$followup\t0";
}
}
else{
if($death>$survivalArr[0]){
$survivalTime="$death\t1";
}
}
}
}
}
my @survivalArr=split(/\t/,$survivalTime);
if($vital_status eq 'Alive'){
if($followup>$survivalArr[0]){
$survivalTime="$followup\t0";
}
}
else{
if($death>$survivalArr[0]){
$survivalTime="$death\t1";
}
}
}
print WF "$patient_barcode\t$survivalTime\t$age\t$gender\t$grade\t$pathologic_stage\t$pathologic_T\t$pathologic_M\t$pathologic_N\n";
}
}
close(RD);
}
}
close(WF);
