TCGA的ID转换可以一步到位了

TCGA的ID转换可以一步到位了 by 生信星球

今天是生信星球陪你的第674天

大神一句话，菜鸟跑半年。我不是大神，但我可以缩短你走弯路的半年~

就像歌儿唱的那样，如果你不知道该往哪儿走，就留在这学点生信好不好~

这里有豆豆和花花的学习历程，从新手到进阶，生信路上有你有我!

0.背景知识

TCGA或TCGA+GTEx的表达矩阵，行名都是ensamble id，因为TCGA数据的参考基因组版本是genecode V22，xena重新分析的TCGA+GTEx数据参考基因组版本则是genecode V23。

代码复制太多次了，于是我写了一个函数，将ensamble id表达矩阵直接转换为symbol。

仍然是tinyarray包，今天说的函数是新写的，到Github下载最新版本的包吧：

https://github.com/xjsun1221/tinyarray

1.输入数据

从xena网站下载的胰腺癌TCGA+GTEx数据，是log过的。

load("Pancreas_expr_ph.Rdata")
Pancreas_expr = 2^Pancreas_expr -1
Pancreas_expr[1:4,1:4]

##                    GTEX-S33H-1226-SM-4AD69
## ENSG00000242268.2                 2.000078
## ENSG00000259041.1                 0.000000
## ENSG00000270112.3                 0.000000
## ENSG00000167578.16              307.002873
##                    GTEX-VJYA-0826-SM-4KL1M
## ENSG00000242268.2                    1.000
## ENSG00000259041.1                    0.000
## ENSG00000270112.3                    0.000
## ENSG00000167578.16                 328.812
##                    TCGA-FB-A545-01
## ENSG00000242268.2           1.0000
## ENSG00000259041.1           0.0000
## ENSG00000270112.3           0.0000
## ENSG00000167578.16        624.6014
##                    GTEX-ZF3C-2026-SM-4WWB5
## ENSG00000242268.2                 0.000000
## ENSG00000259041.1                 0.000000
## ENSG00000270112.3                 3.999903
## ENSG00000167578.16              603.416820

根据数据的列名生成tumor-normal分组信息

library(stringr)
k1 = str_starts(colnames(Pancreas_expr),"TCGA")
k2 = as.numeric(str_sub(colnames(Pancreas_expr),14,15))<10

## Warning: NAs introduced by
## coercion

table(k1&k2)

## 
## FALSE  TRUE 
##   171   179

group_list = ifelse(k1&k2,"tumor","normal")
group_list = factor(group_list,levels = c("normal","tumor"))

2.ID转换

就一个函数，trans_exp搞定。

library(tinyarray)
exp2 = trans_exp(Pancreas_expr)

## 19797 of genes successfully mapping to mRNA,14852 of genes successfully mapping to lncRNA

exp2[1:4,1:4]

##       GTEX-S33H-1226-SM-4AD69
## RAB4B              307.002873
## TIGAR               36.999275
## RNF44             1161.004548
## DNAH3                8.999805
##       GTEX-VJYA-0826-SM-4KL1M
## RAB4B              328.811981
## TIGAR               47.001248
## RNF44             1132.996129
## DNAH3                2.000078
##       TCGA-FB-A545-01
## RAB4B       624.60143
## TIGAR       217.00237
## RNF44      1692.97868
## DNAH3        32.00013
##       GTEX-ZF3C-2026-SM-4WWB5
## RAB4B              603.416820
## TIGAR               51.999249
## RNF44             2222.011968
## DNAH3                2.000078

还有两个可选的参数,mrna_only 和 lncrna_only, mrna_only 设为T则是只挑选mrna，以此类推。

3.指定基因的表达量t检验

这个需求是因为我在重复一篇文章时，看到作者将表达量上调和预后不良联系到一起，给基因取交集。也是做了很多次，所以写了个函数，给定一个基因列表，直接基于T检验返回表达量增加/减少/变化的基因。
相比于使用R包进行差异分析，这个只是最简单的统计学方法，只是定性到底是表达量增加还是减少。

set.seed(12345)
hubgenes = sample(rownames(exp2),20)
#这里是随便挑了20个基因，实际使用应该换成感兴趣的基因
t_choose(hubgenes,exp2,group_list)

##  [1] "SETDB2"       
##  [2] "ZFYVE19"      
##  [3] "C16orf87"     
##  [4] "PIH1D2"       
##  [5] "IQCF2"        
##  [6] "CLIP2"        
##  [7] "ENPP3"        
##  [8] "GSG2"         
##  [9] "VSX1"         
## [10] "HRG"          
## [11] "MGC16275"     
## [12] "RP11.407N17.4"
## [13] "RP11.159D23.2"
## [14] "AP000704.5"   
## [15] "AC107218.3"

t_choose(hubgenes,exp2,group_list,up_only = T)

##  [1] "ENPP3"        
##  [2] "AP000997.1"   
##  [3] "CLIP2"        
##  [4] "SETDB2"       
##  [5] "SYP-AS1"      
##  [6] "CTD-2532K18.1"
##  [7] "RP11-407N17.4"
##  [8] "AP000704.5"   
##  [9] "ZFYVE19"      
## [10] "TET2-AS1"     
## [11] "MGC16275"     
## [12] "AC107218.3"   
## [13] "VSX1"

t_choose(hubgenes,exp2,group_list,down_only = T)

## [1] "RP11-159D23.2"
## [2] "PIH1D2"

插个小广告！
生信零基础入门学习小组
全国巡讲全球听（生信线上直播课）
数据挖掘线上班第4期（两天变三周，实力加量，长期开班）
一起来学单细胞吗？
答疑公告：生信星球答疑公告-2020年全年有效

ixxmu / mp_duty