原因
在平常科研工作中，经常有师兄师姐师弟师妹问我：我现在有一个单基因，我该怎么开展生信研究？出现这个问题的原因是：（1）目前生信研究火热也逐渐受到认可（2）许多医学生在开展实验研究的同时，如果结合生信，则自己的结论和工作量更加吸引到编辑和手审稿人（3）现有的geo、TCGA或者其他免费公开数据库确实是很多研究者的第一选择。
思路
（1）下载整理临床数据、TCGA表达量
（2）单基因的生存分析与临床参数相关分析
（3）单基因的差异分析或者相关分析
（4）单基因的下游通路分析包括GO、KEGG或者通过GSEA
第一节（TCGA生存数据下载）
本节主要下载透明细胞癌KIRC的生存数据
- 加载R包
library(TCGAbiolinks)
library(SummarizedExperiment)
library(dplyr)
library(DT)
rm(list=ls())
setwd('D:\\train\\single_gene')
- 下载生存数据
## ----results='hide', echo=TRUE, message=FALSE, warning=FALSE-------------
clinical <- GDCquery_clinic(project = "TCGA-KIRC", type = "clinical")
## ----echo=TRUE, message=FALSE, warning=FALSE-----------------------------
datatable(clinical, filter = 'top',
options = list(scrollX = TRUE, keys = TRUE, pageLength = 5),
rownames = FALSE)
我们可以看到从上到下共计有537个样本，而且该临床数据有37列。当然我们这里主要关注生存相关的信息比如生存时间和生存状态。
- 整理TCGA肾透明细胞癌的生存时间和生存状态
rm(list=ls())
query <- GDCquery(project = "TCGA-KIRC",
data.category = "Clinical",
file.type = "xml")
GDCdownload(query)
clinical <- GDCprepare_clinic(query, clinical.info = "patient")
clinical_trait <- clinical %>%
dplyr::select(bcr_patient_barcode,gender,vital_status,
days_to_death,days_to_last_followup,race_list,
person_neoplasm_cancer_status,age_at_initial_pathologic_diagnosis,
laterality,neoplasm_histologic_grade,stage_event_pathologic_stage,
stage_event_tnm_categories ) %>%
distinct( bcr_patient_barcode, .keep_all = TRUE)
#整理死亡患者的临床信息
dead_patient <- clinical_trait %>%
dplyr::filter(vital_status == 'Dead') %>%
dplyr::select(-days_to_last_followup) %>%
reshape::rename(c(bcr_patient_barcode = 'Barcode',
gender = 'Gender',
vital_status = 'OS',
days_to_death='OS.Time',
race_list = 'Race',
person_neoplasm_cancer_status='cancer_status',
age_at_initial_pathologic_diagnosis = 'Age',
neoplasm_histologic_grade = 'Grade',
stage_event_pathologic_stage = 'Stage',
stage_event_tnm_categories = 'TNM' )) %>%
mutate(OS=ifelse(OS=='Dead',1,0))%>%
mutate(OS.Time=OS.Time/365)
#整理生存患者的临床信息
alive_patient <- clinical_trait %>%
dplyr::filter(vital_status == 'Alive') %>%
dplyr::select(-days_to_death) %>%
reshape::rename(c(bcr_patient_barcode = 'Barcode',
gender = 'Gender',
vital_status = 'OS',
days_to_last_followup='OS.Time',
race_list = 'Race',
person_neoplasm_cancer_status='cancer_status',
age_at_initial_pathologic_diagnosis = 'Age',
neoplasm_histologic_grade = 'Grade',
stage_event_pathologic_stage = 'Stage',
stage_event_tnm_categories = 'TNM' )) %>%
mutate(OS=ifelse(OS=='Dead',1,0))%>%
mutate(OS.Time=OS.Time/365)
#合并两类患者，得到肾透明细胞癌的临床信息
survival_data <- rbind(dead_patient,alive_patient)
write.csv(survival_data , file = 'KIRC_survival.csv')
最终得到的生存信息，其中包含样本ID，性别，OS（生存状态）、OS（生存时间）、种族、年龄、位置、分级、分期、TNM等信息
image.png
第二节 TCGA表达量下载
- 我这里以肾透明细胞KIRC为例，下载其表达量数据。
exp <- GDCquery(project = "TCGA-KIRC",
data.category = "Transcriptome Profiling",
data.type = "Gene Expression Quantification",
workflow.type = "HTSeq - FPKM")
GDCdownload(exp)
expdat <- GDCprepare(query =exp)
count_matrix= as.data.frame(assay(expdat))
得到的count_matrix矩阵即为我们所需要的KIRC表达量矩阵。
其中每一列为一个样本，每一行为一个基因，我们看到共计56602个基因（包括mRNA和lncRNA等），611个样本（包括肿瘤和癌旁样本）
第三节：TCGA数据库的TPM计算
setwd("D:\\Originaldata\\GRCH\\Homo_sapiens.GRCh38.90")
load("gtf_df.Rda")
test <- gtf_df[1:5,]
View(test)
# =======================================================
setwd('D:\\train\\single_gene')
fpkmToTpm <- function(fpkm)
{
exp(log(fpkm) - log(sum(fpkm)) + log(1e6))
}
expr <- as.data.frame (apply(count_matrix , 2, fpkmToTpm))
expr <- expr %>% rownames_to_column("gene_id")
第四节：表达矩阵中提取mRNA表达矩阵
# =======================================================
mRNA_exprSet <- gtf_df %>%
dplyr::filter(type=="gene",gene_biotype=="protein_coding") %>%
dplyr::select(c(gene_name,gene_id,gene_biotype)) %>%
dplyr::inner_join(expr,by ="gene_id") %>%
tidyr::unite(gene_id,gene_name,gene_id,gene_biotype,sep = " | ")
save(mRNA_exprSet,file = "mRNA_exprSet.Rda")
mRNA_exprSet <- mRNA_exprSet %>%
tidyr::separate(gene_id, c("gene_name","gene_id","gene_biotype"),
sep = " \\| ")
mRNA_exprSet <- mRNA_exprSet[,-(2:3)]
index <- duplicated(mRNA_exprSet$gene_name)
mRNA_exprSet <- mRNA_exprSet[!index,]
row.names(mRNA_exprSet) <- mRNA_exprSet$gene_name
mRNA_exprSet$gene_name <- NULL
第五节：删除癌旁样本和二次测序的样本
#=======================================================
metadata <- data.frame(colnames(mRNA_exprSet))
for (i in 1:length(metadata[,1])) {
num <- as.numeric(as.character(substring(metadata[i,1],14,15)))
if (num == 1 ) {metadata[i,2] <- "T"}
if (num != 1) {metadata[i,2] <- "N"}
}
names(metadata) <- c("id","group")
metadata$group <- as.factor(metadata$group)
metadata <- subset(metadata,metadata$group == "T")
metadata
mRNA_exprSet1 <- mRNA_exprSet[,which(colnames(mRNA_exprSet) %in% metadata$id)]
metadata <- data.frame(colnames(mRNA_exprSet1))
for (i in 1:length(metadata[,1])) {
chr <- as.character(substring(metadata[i,1],22,25))
if ( chr == 'A277' ) {metadata[i,2] <- "Rep"}
if ( chr!= 'A277' ) {metadata[i,2] <- "T"}
}
names(metadata) <- c("id","group")
metadata$group <- as.factor(metadata$group)
metadata <- subset(metadata,metadata$group == "T")
metadata
第六节：保存mRNA表达矩阵
mRNA_exprSet2 <- mRNA_exprSet1[,which(colnames(mRNA_exprSet1) %in% metadata$id)]
write.csv(mRNA_exprSet2,file="KIRC_mRNA_exprSet.csv")
网友评论