抓出那些重复的基因
我们的生信入门班和数据挖掘线上直播课程已经有了三年多的历史，培养了一波又一波优秀的生信人才。课堂上设置的练习题代表着以目前学到的知识所能完成实战，学员们对待练习题的认真程度非常值得学习。虽然有基础的人来看，可能会觉得解法冗余，但人人都是这样学过来的，只不过学完以后忘记了怎么学会的，思考和知识内化的过程值的记录~
下面来看优秀学员棋魂的分享：
R语言练习题
(生信技能树优秀学员棋魂)
- 数据挖掘（GEO,TCGA,单细胞）2022年6月场，快速了解一些生物信息学应用图表
- 生信入门课-2022年6月场，你的生物信息学第一课
题目：g是一个长度为100的字符型向量，元素为基因名，有重复值。有多少个基因在向量g中出现了不止一次？要求：用函数计算出重复基因的具体个数
先看一下g的内容
>load("g.Rdata")
>g
[1] "WNT9B" "CRAMP1L" "MYL12B" "PRSS8" "GFM2"
[6] "CRAMP1L" "TUBA4A" "SLCO1C1" "NYNRIN" "COMMD1"
[11] "COMMD1" "CCT4" "AC017081.1" "RAB7A" "CCT4"
[16] "ZDHHC16" "CASKIN2" "MYL12B" "GGT7" "SNRPE"
[21] "RGPD3" "ZNF586" "COMMD1" "GGT7" "URB1"
[26] "RAB7A" "MPP2" "AFG3L2" "URB1" "AC104581.1"
[31] "IL19" "MPP2" "SYT6" "ATP2A2" "IL19"
[36] "SNRPE" "ARHGAP1" "PRSS8" "PNMT" "ZNF461"
[41] "OR2D3" "CECR5" "SPDL1" "CLEC17A" "ZNF461"
[46] "ATG10" "ATG10" "ATG10" "ZDHHC16" "SLC25A25"
[51] "TCP10" "KRTAP4-3" "SLC30A9" "SLCO1C1" "UBAC1"
[56] "GGT7" "CASKIN2" "GSTP1" "PRY" "UBAC1"
[61] "MPP2" "NYNRIN" "INTS12" "MYL12B" "MPP2"
[66] "KCND1" "RGPD3" "RGPD3" "SLC30A9" "C10orf128"
[71] "HBD" "SLC30A9" "MYL12B" "GGT7" "HEPH"
[76] "TUBA4A" "RP5-1021I20.4" "KLHDC8A" "HBD" "HBP1"
[81] "CCT4" "MARC2" "ZNF586" "LCP1" "CECR5"
[86] "OR2D3" "CRAMP1L" "LIPE" "INTS12" "LIPE"
[91] "NETO2" "CANX" "SPDL1" "ATP6V1B2" "SLCO1C1"
[96] "MARC2" "GGT7" "LCP1" "CECR5" "HOOK2"
解法1：不成熟的解法
因为有重复值，所以先用table()查看有多少重复元素，即出现次数大于1的基因。
###下面为table(g)的返回结果，可以看出每一个基因都给出了出现次数
> table(g)
g
AC017081.1 AC104581.1 AFG3L2 ARHGAP1 ATG10 ATP2A2
1 1 1 1 3 1
ATP6V1B2 C10orf128 CANX CASKIN2 CCT4 CECR5
1 1 1 2 3 3
CLEC17A COMMD1 CRAMP1L GFM2 GGT7 GSTP1
1 3 3 1 5 1
HBD HBP1 HEPH HOOK2 IL19 INTS12
2 1 1 1 2 2
KCND1 KLHDC8A KRTAP4-3 LCP1 LIPE MARC2
1 1 1 2 2 2
MPP2 MYL12B NETO2 NYNRIN OR2D3 PNMT
4 4 1 2 2 1
PRSS8 PRY RAB7A RGPD3 RP5-1021I20.4 SLC25A25
2 1 2 3 1 1
SLC30A9 SLCO1C1 SNRPE SPDL1 SYT6 TCP10
3 3 2 2 1 1
TUBA4A UBAC1 URB1 WNT9B ZDHHC16 ZNF461
2 2 2 1 2 2
ZNF586
2
###查看运行的结果，发现并没有筛选出来，反而将g_3中的重复的基因赋值为TRUE，只出现一次的基因赋值为FALSE
g_3 <- table(g)>1 #将重复出现的元素挑选出来赋值给一个新的向量g_3> g_3
g
AC017081.1 AC104581.1 AFG3L2 ARHGAP1 ATG10 ATP2A2
FALSE FALSE FALSE FALSE TRUE FALSE
ATP6V1B2 C10orf128 CANX CASKIN2 CCT4 CECR5
FALSE FALSE FALSE TRUE TRUE TRUE
CLEC17A COMMD1 CRAMP1L GFM2 GGT7 GSTP1
FALSE TRUE TRUE FALSE TRUE FALSE
HBD HBP1 HEPH HOOK2 IL19 INTS12
TRUE FALSE FALSE FALSE TRUE TRUE
KCND1 KLHDC8A KRTAP4-3 LCP1 LIPE MARC2
FALSE FALSE FALSE TRUE TRUE TRUE
MPP2 MYL12B NETO2 NYNRIN OR2D3 PNMT
TRUE TRUE FALSE TRUE TRUE FALSE
PRSS8 PRY RAB7A RGPD3 RP5-1021I20.4 SLC25A25
TRUE FALSE TRUE TRUE FALSE FALSE
SLC30A9 SLCO1C1 SNRPE SPDL1 SYT6 TCP10
TRUE TRUE TRUE TRUE FALSE FALSE
TUBA4A UBAC1 URB1 WNT9B ZDHHC16 ZNF461
TRUE TRUE TRUE FALSE TRUE TRUE
ZNF586
TRUE
于是通过取子集的方式从g_3中选出值为TRUE的元素，即为重复出现的基因
g_3 <- g_3[g_3==TRUE] # 其实不写==TRUE也一样，就是代码会变得不好理解
> g_3
g
ATG10 CASKIN2 CCT4 CECR5 COMMD1 CRAMP1L GGT7 HBD IL19 INTS12 LCP1
TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
LIPE MARC2 MPP2 MYL12B NYNRIN OR2D3 PRSS8 RAB7A RGPD3 SLC30A9 SLCO1C1
TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
SNRPE SPDL1 TUBA4A UBAC1 URB1 ZDHHC16 ZNF461 ZNF586
TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
统计重复基因的个数，为30个
length(g_3)
> length(g_3)
[1] 30
解法2：解决解法1的冗余
解法1没有直接通过table(g)去取子集，反而中间转化了一步逻辑值，不仅多此一举，而且最后的结果里也不能查看每个基因出现了几次。于是去掉中间的转化，直接取子集。
g_3 <- table(g)[table(g)>1]
> g_3
g
ATG10 CASKIN2 CCT4 CECR5 COMMD1 CRAMP1L GGT7 HBD IL19 INTS12 LCP1
3 2 3 3 3 3 5 2 2 2 2
LIPE MARC2 MPP2 MYL12B NYNRIN OR2D3 PRSS8 RAB7A RGPD3 SLC30A9 SLCO1C1
2 2 4 4 2 2 2 2 3 3 3
SNRPE SPDL1 TUBA4A UBAC1 URB1 ZDHHC16 ZNF461 ZNF586
2 2 2 2 2 2 2 2
> length(g_3)
[1] 30
解法3：继续解决解法1和解法2的问题
虽然优化了的解法2能够挑选出来重复基因，并可以统计每个这样的基因出现的次数，但是g_3（无论解法1还是解法2里）都不是一个简单的向量。仔细观察每个元素都跟着一个逻辑值（解法1中）或者数值（解法2）中。那么有没有可能通过指定重复基因所在的index（即位置）而把它们挑选出来呢？
首先通过as.numeric()尝试将table(g)的出现次数转化为一个单纯的数值型向量。成功！这样数值大于1的元素的index就可以用来进行后面的提取。
as.numeric(table(g))
> as.numeric(table(g))
[1] 1 1 1 1 3 1 1 1 1 2 3 3 1 3 3 1 5 1 2 1 1 1 2 2 1 1 1 2 2 2 4 4 1 2 2 1 2 1 2 3 1 1 3 3
[45] 2 2 1 1 2 2 2 1 2 2 2
接着再次观察table(g)的结果，发现不仅给出了每个基因出现的次数，而且将这些基因按照名称的英文字母顺序进行了排序。这就提示能不能用sort(unqiue(g))再构造一个排好序的单纯的字符型向量呢。尝试后发现成功！sort(unituqe(g))返回的基因排列顺序与table(g)的完全一致。注意：直接使用unqiue(g)返回的元素是按照在g中的顺序排列的。
> sort(unique(g))
[1] "AC017081.1" "AC104581.1" "AFG3L2" "ARHGAP1" "ATG10"
[6] "ATP2A2" "ATP6V1B2" "C10orf128" "CANX" "CASKIN2"
[11] "CCT4" "CECR5" "CLEC17A" "COMMD1" "CRAMP1L"
[16] "GFM2" "GGT7" "GSTP1" "HBD" "HBP1"
[21] "HEPH" "HOOK2" "IL19" "INTS12" "KCND1"
[26] "KLHDC8A" "KRTAP4-3" "LCP1" "LIPE" "MARC2"
[31] "MPP2" "MYL12B" "NETO2" "NYNRIN" "OR2D3"
[36] "PNMT" "PRSS8" "PRY" "RAB7A" "RGPD3"
[41] "RP5-1021I20.4" "SLC25A25" "SLC30A9" "SLCO1C1" "SNRPE"
[46] "SPDL1" "SYT6" "TCP10" "TUBA4A" "UBAC1"
[51] "URB1" "WNT9B" "ZDHHC16" "ZNF461" "ZNF586"
#对比table(g)的结果，发现排列顺序完全吻合。
> table(g) #table()不仅能统计出现次数，还能排序！这是不是一个意外收获？哈哈！
g
AC017081.1 AC104581.1 AFG3L2 ARHGAP1 ATG10 ATP2A2
1 1 1 1 3 1
ATP6V1B2 C10orf128 CANX CASKIN2 CCT4 CECR5
1 1 1 2 3 3
CLEC17A COMMD1 CRAMP1L GFM2 GGT7 GSTP1
1 3 3 1 5 1
HBD HBP1 HEPH HOOK2 IL19 INTS12
2 1 1 1 2 2
KCND1 KLHDC8A KRTAP4-3 LCP1 LIPE MARC2
1 1 1 2 2 2
MPP2 MYL12B NETO2 NYNRIN OR2D3 PNMT
4 4 1 2 2 1
PRSS8 PRY RAB7A RGPD3 RP5-1021I20.4 SLC25A25
2 1 2 3 1 1
SLC30A9 SLCO1C1 SNRPE SPDL1 SYT6 TCP10
3 3 2 2 1 1
TUBA4A UBAC1 URB1 WNT9B ZDHHC16 ZNF461
2 2 2 1 2 2
ZNF586
2
通过前面两步，两个单纯的向量的已经准备好了，下面就可以通过向量按照位置取子集的方式挑出重复出现的基因了。任务完成！
> sort(unique(g))[as.numeric(table(g))>1]
[1] "ATG10" "CASKIN2" "CCT4" "CECR5" "COMMD1" "CRAMP1L" "GGT7" "HBD"
[9] "IL19" "INTS12" "LCP1" "LIPE" "MARC2" "MPP2" "MYL12B" "NYNRIN"
[17] "OR2D3" "PRSS8" "RAB7A" "RGPD3" "SLC30A9" "SLCO1C1" "SNRPE" "SPDL1"
[25] "TUBA4A" "UBAC1" "URB1" "ZDHHC16" "ZNF461" "ZNF586"
> length(sort(unique(g))[as.numeric(table(g))>1])
[1] 30
对比这三种解法，可以看出如果单纯追求满足题目的要求，那么解法2无疑是最好的。但是，如果考虑到后续处理数据的需求，那么我会使用解法3。因为我以为使用as.character(table(g))可以忽略基因出现的次数，而生成一个单纯的基因名的字符型向量时，发现返回的是下面的结果。
> as.character(table(g))
[1] "1" "1" "1" "1" "3" "1" "1" "1" "1" "2" "3" "3" "1" "3" "3" "1" "5" "1" "2" "1" "1" "1"
[23] "2" "2" "1" "1" "1" "2" "2" "2" "4" "4" "1" "2" "2" "1" "2" "1" "2" "3" "1" "1" "3" "3"
[45] "2" "2" "1" "1" "2" "2" "2" "1" "2" "2" "2"
这就说明table(g)仅仅能满足“看一看”的需求，我们不能用它进行后面的数据处理。因此在延续性上导致解法2很弱。
解法4：通过数据框的方式解决
还有一个方法是使用sort(unique(g)和as.character(table(g))生成一个两列的数据框，然后对数据框进行操作。
df <- data.frame(
gene=sort(unique(g)),
times=as.numeric(table(g))
)
#生成的数据框如下，第一列是基因名，第二列是出现次数
> df
gene times
1 AC017081.1 1
2 AC104581.1 1
3 AFG3L2 1
4 ARHGAP1 1
5 ATG10 3
6 ATP2A2 1
7 ATP6V1B2 1
8 C10orf128 1
9 CANX 1
10 CASKIN2 2
11 CCT4 3
12 CECR5 3
13 CLEC17A 1
14 COMMD1 3
15 CRAMP1L 3
16 GFM2 1
17 GGT7 5
18 GSTP1 1
19 HBD 2
20 HBP1 1
21 HEPH 1
22 HOOK2 1
23 IL19 2
24 INTS12 2
25 KCND1 1
26 KLHDC8A 1
27 KRTAP4-3 1
28 LCP1 2
29 LIPE 2
30 MARC2 2
31 MPP2 4
32 MYL12B 4
33 NETO2 1
34 NYNRIN 2
35 OR2D3 2
36 PNMT 1
37 PRSS8 2
38 PRY 1
39 RAB7A 2
40 RGPD3 3
41 RP5-1021I20.4 1
42 SLC25A25 1
43 SLC30A9 3
44 SLCO1C1 3
45 SNRPE 2
46 SPDL1 2
47 SYT6 1
48 TCP10 1
49 TUBA4A 2
50 UBAC1 2
51 URB1 2
52 WNT9B 1
53 ZDHHC16 2
54 ZNF461 2
55 ZNF586 2
另外一个生成数据框的操作方法是直接把table(g）转化为数据框，table函数自带统计和排序功能，可以节省时间。转化后，下面的操作完全相同，但是因为列名不一样，所以不要直接套用。
df_1 <- data.frame(table(g))
#查看一下转化而来的数据框
> head(df_1)
g Freq
1 AC017081.1 1
2 AC104581.1 1
3 AFG3L2 1
4 ARHGAP1 1
5 ATG10 3
6 ATP2A2 1
#可以把列名改一下，显得更正式些
colnames(df_2)[1] <- "Gene"
> colnames(df_2)
[1] "Gene" "Freq"
生成数据框之后，只需要把times这一列（或者Freq列）里大于1的行筛选出来，就可以满足要求了。
df_1 <- df1[df1$times>1,]#生成一个中间变量df_1是为了便于观察和对比
#查看df_1，发现满足要求的基因全部挑出来了
> df_1
gene times
5 ATG10 3
10 CASKIN2 2
11 CCT4 3
12 CECR5 3
14 COMMD1 3
15 CRAMP1L 3
17 GGT7 5
19 HBD 2
23 IL19 2
24 INTS12 2
28 LCP1 2
29 LIPE 2
30 MARC2 2
31 MPP2 4
32 MYL12B 4
34 NYNRIN 2
35 OR2D3 2
37 PRSS8 2
39 RAB7A 2
40 RGPD3 3
43 SLC30A9 3
44 SLCO1C1 3
45 SNRPE 2
46 SPDL1 2
49 TUBA4A 2
50 UBAC1 2
51 URB1 2
53 ZDHHC16 2
54 ZNF461 2
55 ZNF586 2
到上一步其实已经满足要求了，不过最后还可以做下面两步来个心满意足的收尾。
#先检查一下是不是挑出了30个重复出现的基因
> dim(df_1)
[1] 30 2
#再对数据框按照基因出现次数从少到多排序，看着更舒服些
> df_1[order(df_1$times),]
gene times
10 CASKIN2 2
19 HBD 2
23 IL19 2
24 INTS12 2
28 LCP1 2
29 LIPE 2
30 MARC2 2
34 NYNRIN 2
35 OR2D3 2
37 PRSS8 2
39 RAB7A 2
45 SNRPE 2
46 SPDL1 2
49 TUBA4A 2
50 UBAC1 2
51 URB1 2
53 ZDHHC16 2
54 ZNF461 2
55 ZNF586 2
5 ATG10 3
11 CCT4 3
12 CECR5 3
14 COMMD1 3
15 CRAMP1L 3
40 RGPD3 3
43 SLC30A9 3
44 SLCO1C1 3
31 MPP2 4
32 MYL12B 4
17 GGT7 5
本文参与 腾讯云自媒体同步曝光计划，分享自微信公众号。
原始发表：2022-05-21，如有侵权请联系 cloudcommunity@tencent.com 删除
评论
登录后参与评论
推荐阅读
目录
相关产品与服务
云直播
云直播（Cloud Streaming Services，CSS）为您提供极速、稳定、专业的云端直播处理服务，根据业务的不同直播场景需求，云直播提供了标准直播、快直播、云导播台三种服务，分别针对大规模实时观看、超低延时直播、便捷云端导播的场景，配合腾讯云视立方·直播 SDK，为您提供一站式的音视频直播解决方案。
腾讯云开发者
