一文解决大量基因的生存分析并作图
前言
生信文章中，比较常见的套路是批量寻找与生存（一般是总体生存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')本文参与 腾讯云自媒体分享计划，分享自作者个人站点/博客。
原始发表：2019.05.31 ，如有侵权请联系 cloudcommunity@tencent.com 删除
评论
登录后参与评论
推荐阅读
目录