前言

生信文章中，比较常见的套路是批量寻找与生存（一般是总体生存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')