单基因生信分析流程(1)一文解决TCGA数据下载整理问题

原因

在平常科研工作中,经常有师兄师姐师弟师妹问我:我现在有一个单基因,我该怎么开展生信研究?出现这个问题的原因是:(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")

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

扫码关注云+社区

领取腾讯云代金券