零基础入门转录组分析-数据处理（TCGA数据库）教程配套的代码，原始数据以及最终处理好的数据

rm(list = ls()) # 删除工作空间中所有的对象 setwd('/data/nas1/xuguangxin/test/教程/') # 设置工作路径 if(!dir.exists('./Lesson_1')){ dir.create('./Lesson_1') } # 判断该工作路径下是否存在名为00_rawdata的文件夹，如果不存在则创建，如果存在则pass setwd('./Lesson_1/') # 设置路径到刚才新建的00_rawdata下 library(tidyverse) library(dplyr) library(rtracklayer) ### 读入下载的数据 tcga_count <- read_tsv(file = './TCGA-LUAD.htseq_counts.tsv.gz') # count数据 tcga_fpkm <- read_tsv(file = "./TCGA-LUAD.htseq_fpkm.tsv.gz") # fpkm数据 tcga_survival <- read_tsv(file = "./TCGA-LUAD.survival.tsv") # 患者生存状态 tcga_phenotype <- read_tsv(file = "./TCGA-LUAD.GDC_phenotype.tsv.gz") # 患者临床信息 ### count数据处理 tcga_count <- as.data.frame(tcga_count) # 将导入的文件转成数据框格式 tcga_count <- tcga_count[!duplicated(tcga_count$Ensembl_ID), ] tcga_count <- column_to_rownames(tcga_count, var = 'Ensembl_ID') # 将文件第一列转为行名 tcga_count <- 2^tcga_count-1 # xena下载的数据经过了log2+1转化，需要将其还原 gene_annotation <- import('./gencode.v22.annotation.gtf.gz') gene_annotation <- as.data.frame(gene_annotation)#将文件转换为数据框格式 gene_annotation <- gene_annotation [gene_annotation$type == 'gene', ] # 筛选为gene的类型 gene_annotation <- dplyr::select(gene_annotation, gene_id, gene_name, seqnames, start, end, width, strand, gene_type) gene_symbol <- gene_annotation[,(1:2)] # 筛选gene_annotation文件中的第一列和第二列 colnames(gene_symbol) <- c("ID", "symbol") # 将gene_symbol文件的列名改成ID和symbol data_tcga <- tcga_count data_tcga$ID <- rownames(data_tcga) # 给data_tcga添加一列，列名为ID data_tcga$ID <- as.character(data_tcga$ID) # 将data_tcga文件中ID列从数值类型数据转化为字符串类型数据 gene_symbol$ID <- as.character(gene_symbol$ID) # 将gene_symbol文件中ID列从数值类型数据转化为字符串类型数据 data_tcga <- inner_join(data_tcga, gene_symbol, by = "ID") # 将data_tcga文件和gene_symbol文件根据ID列进行合并（这样就能获得ID对应的基因名，但是都排在最后） data_tcga <- dplyr::select(data_tcga, -ID) # 删除ID列 data_tcga <- dplyr::select(data_tcga, symbol, everything()) # 重新排列,将最后一列的基因名放到第一列（如果不加everything只会选择symbol一列） data_tcga <- mutate(data_tcga, rowMean=rowMeans(data_tcga[grep('TCGA',names(data_tcga))])) # 添加一列为表达量的平均值 data_tcga <- arrange(data_tcga, desc(rowMean)) # 把表达量的平均值从大到小排序 data_tcga <- distinct(data_tcga, symbol, .keep_all = T) # distinct函数被用于去重，.keep_all参数表示去重后返回数据框的所有列向量 data_tcga <- dplyr::select(data_tcga, -rowMean) # 去除rowMean这一列 data_tcga <- tibble::column_to_rownames(data_tcga, var = "symbol") ## 把第一列变成行名并删除 dim(data_tcga) ### 筛选癌症和对照样本。01-09为肿瘤，10-19为正常对照 mete <- data.frame(colnames(data_tcga)) for (i in 1:length(mete[, 1])) { num <- as.numeric(as.character(substring(mete[i, 1], 14, 15))) # 提取每行第一列中第14，15字符 if(num %in% seq(1, 9)){mete[i, 2] = "T"} # 判断：如果提取的数字在0-9之内就在每行第二列加上T表示肿瘤样本 if(num %in% seq(10, 29)){mete[i, 2] = "N"} # 判断：如果提取的数字在10-29之内就在每行第二列加上N表示正常对照 } colnames(mete) <- c("id", "group") table(mete$group) mete$group <- as.factor(mete$group) # 将mete中group列转为因子模式 mete <- subset(mete, mete$group == "T") # subset函数，从数据框中选择出group组为T的行 exp_tumor <- data_tcga[, which(colnames(data_tcga) %in% mete$id)] # 从大样本中选出组为T的样本 exp_tumor <- as.data.frame(exp_tumor) exp_tumor <- exp_tumor[, colnames(exp_tumor) %in% tcga_survival$sample] # 进一步筛选具有生存信息的样本 exp_control <- data_tcga[, which(!colnames(data_tcga) %in% mete$id)] # 反向从大样本中选出组为N的样本 exp_control <- as.data.frame(exp_control) exp_control <- exp_control[, colnames(exp_control) %in% tcga_survival$sample] # 进一步筛选具有生存信息的样本 data_final <- cbind(exp_control, exp_tumor) # 将肿瘤样本文件和正常样本文件合并 data_final <- na.omit(data_final) # 去除含有NA的行 table(tcga_phenotype$sample_type.samples) # 查看样本类型，一共四种：FFPE Scrolls，Primary Tumor，Recurrent Tumor，Solid Tissue Normal，这里面我们选择Primary Tumor——原发肿瘤和Solid Tissue Normal——正常组织 table(tcga_phenotype$primary_site) # 查看癌症发生的部位 clinical_data <- subset(tcga_phenotype, sample_type.samples == 'Primary Tumor' | sample_type.samples == 'Solid Tissue Normal') clinical_data <- clinical_data[clinical_data$submitter_id.samples %in% colnames(data_final), ] clinical_data <- clinical_data[order(clinical_data$sample_type.samples, decreasing = T), ] Group <- data.frame(sample = clinical_data$submitter_id.samples, group = clinical_data$sample_type.samples) Group$group <- ifelse(Group$group == 'Solid Tissue Normal', 'control', 'disease') data_final <- subset(data_final, select = Group$sample) # 确保前面整理的count数据与整理的分组信息的样本是一致的 table(Group$group) ##control:59 disease:511 tcga_survival <- tcga_survival[, -3] tcga_survival <- tcga_survival[tcga_survival$sample %in% Group$sample, ] write.csv(tcga_survival, './data_survival.csv') write.csv(clinical_data, './data_clinical.csv') write.csv(data_final, file = './data_count.csv') write.csv(Group, file = './data_group.csv') ### fpkm数据整理 tcga_fpkm <- as.data.frame(tcga_fpkm) tcga_fpkm <- tcga_fpkm[!duplicated(tcga_fpkm$Ensembl_ID), ] tcga_fpkm <- column_to_rownames(tcga_fpkm, var = 'Ensembl_ID') tcga_fpkm <- 2^tcga_fpkm-1 dat_fpkm <- tcga_fpkm dat_fpkm$ID <- rownames(dat_fpkm) dat_fpkm$ID <- as.character(dat_fpkm$ID) gene_symbol$ID <- as.character(gene_symbol$ID) dat_fpkm<-dat_fpkm %>% inner_join(gene_symbol,by='ID')%>% select(-ID)%>% ## 去除多余信息 select(symbol,everything())%>% ## 重新排列 mutate(rowMean=rowMeans(.[grep('TCGA',names(.))]))%>% ## 求出平均数 arrange(desc(rowMean))%>% ## 把表达量的平均值从大到小排序 distinct(symbol,.keep_all = T)%>% ## symbol留下第一个 select(-rowMean)%>% ## 反向选择去除rowMean这一列 tibble::column_to_rownames(colnames(.)[1]) ## 把第一列变成行名并删除 dim(dat_fpkm) dat_fpkm <- dat_fpkm[rownames(data_final), ] dat_fpkm <- dat_fpkm[, colnames(data_final)] dat_fpkm <- na.omit(dat_fpkm) dat_fpkm <- log2(dat_fpkm + 1) write.csv(dat_fpkm, file = './data_fpkm.csv')