一文解决大量基因的生存分析并作图

前言

生信文章中,比较常见的套路是批量寻找与生存(一般是总体生存OS)相关的标志物。

这两篇纯生信文章都是对单个基因或者所有单个marker做生存分析,目的是找到其中能够影响患者生存的marker或者基因(包括miRNA,lncRNA,mRNA等等)。这也是目前非常常见的筛选基因或者marker的方法。

下载TCGA数据(表达量)

# ===============================================================

#load package

# ===============================================================

library(TCGAbiolinks)
library(SummarizedExperiment)
library(dplyr)
library(DT)
library(tibble)
library(tidyr)
rm(list=ls())


# ===============================================================

#通过TCGAbiolinks下载表达量数据

# ===============================================================

setwd('D:\\train\\lasso')


exp <- GDCquery(project = "TCGA-COAD", 
                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))



# ===============================================================

#注释基因并计算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")

得到的count_matrix矩阵即为我们所需要的表达量矩阵。其中每一列为一个样本,每一行为一个基因

此外,我们对下载表达矩阵(FPKM格式)进行转化成TPM格式的表达量数。这一步的理由是因为目前认为TPM是最能够反映基因真实表达量的指标(相对于raw count、FPKM、RPKM等等)

TPM解释

Transcripts Per Kilobase of exonmodel per Million mapped reads (每千个碱基的转录每百万映射读取的Transcripts),优化的RPKM计算方法,可以用于同一物种不同组织的比较。TPM 的计算公式:

TPMi={( Ni/Li )*1000000 } / sum( Ni/Li+……..+ Nm/Lm ) Ni:mapping到基因i上的read数; Li:基因i的外显子长度的总和。 在一个样本中一个基因的TPM:先对每个基因的read数用基因的长度进行校正,之后再用校正后的这个基因read数(Ni/Li)与校正后的这个样本的所有read数(sum(Ni/Li+……..+ Nm/Lm))求商。由此可知,TPM概括了基因的长度、表达量和基因数目。TPM可以用于同一物种不同组织间的比较,因为sum值总是唯一的。

提取mRNA

#提取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")

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

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

mRNA_exprSet2 <- mRNA_exprSet1[,which(colnames(mRNA_exprSet1) %in% metadata$id)]





#提取非重复样本
metadata <- data.frame(colnames(mRNA_exprSet2))

for (i in 1:length(metadata[,1])) {
  chr <-  as.character(substring(metadata[i,1],16,16))
  if ( chr == 'A' ) {metadata[i,2] <- "T"}
  if ( chr!= 'A' ) {metadata[i,2] <- "un"}
}

names(metadata) <- c("id","group")

metadata$group <- as.factor(metadata$group)

metadata <- subset(metadata,metadata$group == "T")

mRNA_exprSet3 <- mRNA_exprSet2[,which(colnames(mRNA_exprSet2) %in% metadata$id)]

setwd('D:\\train\\lasso')

save(mRNA_exprSet3,file='mRNA_exprSet3.Rdata')

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

发表于

我来说两句

0 条评论
登录 后参与评论

扫码关注云+社区

领取腾讯云代金券