本文首發(fā)于公眾號(hào)：醫(yī)學(xué)和生信筆記
醫(yī)學(xué)和生信筆記，專注R語(yǔ)言在臨床醫(yī)學(xué)中的使用，R語(yǔ)言數(shù)據(jù)分析和可視化。主要分享R語(yǔ)言做醫(yī)學(xué)統(tǒng)計(jì)學(xué)、meta分析、網(wǎng)絡(luò)藥理學(xué)、臨床預(yù)測(cè)模型、機(jī)器學(xué)習(xí)、生物信息學(xué)等。
現(xiàn)在使用TCGAbiolinks下載轉(zhuǎn)錄組數(shù)據(jù)后，直接是一個(gè)SummarizedExperiment對(duì)象，這個(gè)對(duì)象非常重要且好用。因?yàn)槔锩?直接包含了表達(dá)矩陣、樣本信息、基因信息，可以非常方便的通過內(nèi)置函數(shù)直接提取想要的數(shù)據(jù)，再也不用手扒了！!
這個(gè)對(duì)象的結(jié)構(gòu)是這樣的：
上次我們下載了常見的組學(xué)數(shù)據(jù)，今天學(xué)習(xí)下怎么提取數(shù)據(jù)，就以TCGA-READ的轉(zhuǎn)錄組數(shù)據(jù)為例。
分別提取mRNA和lncRNA的表達(dá)矩陣，還要添加gene symbol的那種！
加載數(shù)據(jù)和R包
加載之前下載好的數(shù)據(jù)。
rm(list = ls())
library(SummarizedExperiment)
load("TCGA-mRNA/TCGA-READ_mRNA.Rdata")
se <- data
這個(gè)se就是你的對(duì)象，含有coldata, rowdata， meta-data，以及最重要的assay，共有6個(gè)assay
探索SummarizedExperiment對(duì)象
se
## class: RangedSummarizedExperiment
## dim: 60660 177
## metadata(1): data_release
## assays(6): unstranded stranded_first ... fpkm_unstrand fpkm_uq_unstrand
## rownames(60660): ENSG00000000003.15 ENSG00000000005.6 ...
## ENSG00000288674.1 ENSG00000288675.1
## rowData names(10): source type ... hgnc_id havana_gene
## colnames(177): TCGA-AG-3580-01A-01R-0821-07
## TCGA-AF-2692-11A-01R-A32Z-07 ... TCGA-AG-3894-01A-01R-1119-07
## TCGA-AG-3574-01A-01R-0821-07
## colData names(107): barcode patient ... paper_vascular_invasion_present
## paper_vital_status
看看這個(gè)對(duì)象，它告訴你：
- 類型是：RangedSummarizedExperiment
- 維度：60660行，177列
- 6個(gè)assay以及它們的名字
- 表達(dá)矩陣的行名
- 行信息（也就是基因信息），比如gene id,gene name,gene type
- 表達(dá)矩陣的列名（也就是樣本名）
- 列信息，也就是樣本信息，比如生存時(shí)間、生存狀態(tài)這些
太齊全了有沒有?。?br> 每個(gè)assay你可以理解為一個(gè)表達(dá)矩陣，我們需要的counts矩陣、TPM矩陣、FPKM矩陣就是其中一個(gè)~
# 查看每個(gè)assay的名字
names(assays(se))
## [1] "unstranded" "stranded_first" "stranded_second" "tpm_unstrand"
## [5] "fpkm_unstrand" "fpkm_uq_unstrand"
每個(gè)基因?qū)儆趍RNA還是lncRNA存儲(chǔ)在rowData中，這個(gè)rowData你可以理解為一個(gè)包含基因各種信息的數(shù)據(jù)框。
其中gene_type是基因類型，幫助我們區(qū)分到底是lncRNA還是mRNA，當(dāng)然還包括很多其他類型。
# 提取rowData
rowdata <- rowData(se)
# 看看rowData包括哪些內(nèi)容，可以看到里面有我們需要的gene_name和gene_type
names(rowdata)
## [1] "source" "type" "score" "phase" "gene_id"
## [6] "gene_type" "gene_name" "level" "hgnc_id" "havana_gene"
# gene_type是基因類型，看看有哪些
table(rowdata$gene_type)
##
## IG_C_gene IG_C_pseudogene
## 14 9
## IG_D_gene IG_J_gene
## 37 18
## IG_J_pseudogene IG_pseudogene
## 3 1
## IG_V_gene IG_V_pseudogene
## 145 187
## lncRNA miRNA
## 16901 1881
## misc_RNA Mt_rRNA
## 2212 2
## Mt_tRNA polymorphic_pseudogene
## 22 48
## processed_pseudogene protein_coding
## 10167 19962
## pseudogene ribozyme
## 18 8
## rRNA rRNA_pseudogene
## 47 497
## scaRNA scRNA
## 49 1
## snoRNA snRNA
## 943 1901
## sRNA TEC
## 5 1057
## TR_C_gene TR_D_gene
## 6 4
## TR_J_gene TR_J_pseudogene
## 79 4
## TR_V_gene TR_V_pseudogene
## 106 33
## transcribed_processed_pseudogene transcribed_unitary_pseudogene
## 500 138
## transcribed_unprocessed_pseudogene translated_processed_pseudogene
## 939 2
## translated_unprocessed_pseudogene unitary_pseudogene
## 1 98
## unprocessed_pseudogene vault_RNA
## 2614 1
gene_name就是gene_symbol，我們的id轉(zhuǎn)換就用這一列信息。
# gene_name就是gene_symbol，就是我們需要的
head(rowdata$gene_name)
## [1] "TSPAN6" "TNMD" "DPM1" "SCYL3" "C1orf112" "FGR"
length(rowdata$gene_name)
## [1] 60660
還有一個(gè)重要的知識(shí)點(diǎn)是：SummarizedExperiment對(duì)象可以取子集，就像對(duì)數(shù)據(jù)框取子集那樣，選擇符合條件的行和列，并且子集也是SummarizedExperiment對(duì)象！
rowdata <- rowData(se)
# 分別提取mRNA和lncRNA的SummarizedExperiment對(duì)象
# 根據(jù)gene_type取子集，太簡(jiǎn)單了！
se_mrna <- se[rowdata$gene_type == "protein_coding",]
se_lnc <- se[rowdata$gene_type == "lncRNA"]
se_mrna #還是一個(gè)SummarizedExperiment對(duì)象，神奇！
## class: RangedSummarizedExperiment
## dim: 19962 177
## metadata(1): data_release
## assays(6): unstranded stranded_first ... fpkm_unstrand fpkm_uq_unstrand
## rownames(19962): ENSG00000000003.15 ENSG00000000005.6 ...
## ENSG00000288674.1 ENSG00000288675.1
## rowData names(10): source type ... hgnc_id havana_gene
## colnames(177): TCGA-AG-3580-01A-01R-0821-07
## TCGA-AF-2692-11A-01R-A32Z-07 ... TCGA-AG-3894-01A-01R-1119-07
## TCGA-AG-3574-01A-01R-0821-07
## colData names(107): barcode patient ... paper_vascular_invasion_present
## paper_vital_status
se_lnc
## class: RangedSummarizedExperiment
## dim: 16901 177
## metadata(1): data_release
## assays(6): unstranded stranded_first ... fpkm_unstrand fpkm_uq_unstrand
## rownames(16901): ENSG00000082929.8 ENSG00000083622.8 ...
## ENSG00000288667.1 ENSG00000288670.1
## rowData names(10): source type ... hgnc_id havana_gene
## colnames(177): TCGA-AG-3580-01A-01R-0821-07
## TCGA-AF-2692-11A-01R-A32Z-07 ... TCGA-AG-3894-01A-01R-1119-07
## TCGA-AG-3574-01A-01R-0821-07
## colData names(107): barcode patient ... paper_vascular_invasion_present
## paper_vital_status
提取表達(dá)矩陣
有了這些東西，就可以提取表達(dá)矩陣了，直接使用assay()搞定！
# mRNA的counts矩陣
expr_counts_mrna <- assay(se_mrna,"unstranded")
# mRNA的tpm矩陣
expr_tpm_mrna <- assay(se_mrna,"tpm_unstrand")
# mRNA的fpkm矩陣
expr_fpkm_mrna <- assay(se_mrna,"fpkm_unstrand")
# lncRNA的counts矩陣
expr_counts_lnc <- assay(se_lnc,"unstranded")
# lncRNA的tpm矩陣
expr_tpm_lnc <- assay(se_lnc,"tpm_unstrand")
# lncRNA的fpkm矩陣
expr_fpkm_lnc <- assay(se_lnc,"fpkm_unstrand")
簡(jiǎn)單！方便！快捷！
隨便展示下：
expr_counts_mrna[1:10,1:2]
## TCGA-AG-3580-01A-01R-0821-07 TCGA-AF-2692-11A-01R-A32Z-07
## ENSG00000000003.15 3199 5839
## ENSG00000000005.6 6 91
## ENSG00000000419.13 828 1867
## ENSG00000000457.14 386 639
## ENSG00000000460.17 228 289
## ENSG00000000938.13 130 452
## ENSG00000000971.16 277 5170
## ENSG00000001036.14 1648 2946
## ENSG00000001084.13 823 2414
## ENSG00000001167.14 619 1487
添加gene_symbol
添加gene_symbol也就非常簡(jiǎn)單了，只要提取gene_name這一列，然后和原來(lái)的表達(dá)矩陣合并即可！
# 先提取gene_name
symbol_mrna <- rowData(se_mrna)$gene_name
head(symbol_mrna)
## [1] "TSPAN6" "TNMD" "DPM1" "SCYL3" "C1orf112" "FGR"
symbol_lnc <- rowData(se_lnc)$gene_name
head(symbol_lnc)
## [1] "LINC01587" "AC000061.1" "AC016026.1" "IGF2-AS" "RRN3P2"
## [6] "AC087235.1"
和你喜歡的表達(dá)矩陣合并就行了：
expr_counts_mrna_symbol <- cbind(data.frame(symbol_mrna),
as.data.frame(expr_counts_mrna))
非常順利~
但是呢，此時(shí)gene_symbol是有重復(fù)的，看上圖中就有2個(gè)CD99，需要去重！
去重復(fù)也很簡(jiǎn)單，這里我們保留最大的那個(gè)。
suppressPackageStartupMessages(library(tidyverse))
expr_read <- expr_counts_mrna_symbol %>%
as_tibble() %>% # tibble不支持row name，我竟然才發(fā)現(xiàn)！
mutate(meanrow = rowMeans(.[,-1]), .before=2) %>%
arrange(desc(meanrow)) %>%
distinct(symbol_mrna,.keep_all=T) %>%
select(-meanrow) %>%
column_to_rownames(var = "symbol_mrna") %>%
as.data.frame()
不過還是要注意，gene_symbol是有重復(fù)的，需要去重復(fù)哦~
結(jié)果就變成大家最熟悉的表達(dá)矩陣了：
這樣一個(gè)表達(dá)矩陣就搞定了！
本文首發(fā)于公眾號(hào)：醫(yī)學(xué)和生信筆記
醫(yī)學(xué)和生信筆記，專注R語(yǔ)言在臨床醫(yī)學(xué)中的使用，R語(yǔ)言數(shù)據(jù)分析和可視化。主要分享R語(yǔ)言做醫(yī)學(xué)統(tǒng)計(jì)學(xué)、meta分析、網(wǎng)絡(luò)藥理學(xué)、臨床預(yù)測(cè)模型、機(jī)器學(xué)習(xí)、生物信息學(xué)等。
