一键解锁！用R脚本轻松统计SUPPA2识别的可变剪接（AS）事件结果

关注我们并在后台回复 “进群”，即可加入农心生信工作室学习交流群，群内不定时分享源代码及示例文件，并在线交流答疑。我们在农心生信，等你！

由于微信改版，乱序推送，很多朋友反映收不到公众号推文。快跟着图片步骤，将农心生信公众号设为星标，不错过每一条精彩内容！！

本公众号已开通留言功能，欢迎大家直接在文末留言交流！

写在前面

可变剪接（Alternative splicing，AS）是真核生物中普遍存在的转录后调控过程，使得同一基因能够产生不同的信使RNA（mRNA）异构体。AS可以增强生物基因的功能多样性，并且参与调节基因表达。SUPPA2（https://github.com/comprna/SUPPA）是一个可以根据转录本结构注释文件（GTF）识别AS事件的软件，基于python编写，它的安装和使用都非常方便，这里不再赘述，大家可以详细阅读SUPPA的GitHub。今天我们要讨论的问题是，如何在利用SUPPA完成AS的鉴定后，对鉴定结果进行统计，包括发生了多少AS事件，不同类别的AS事件分别有多少，多少基因受AS事件影响，产生了多少不同的转录本。下面，我们将用一个R脚本来完成AS事件统计。

代码详解

首先，我们需要用awk来合并SUPPA产生的所有结果：

awk 'FNR==1 && NR!=1 { while (/^<header>/) getline; } 1 {print}' *.ioe > all.events.ioe

接下来，我们基于R包tidyverse来完成AS事件的统计与可视化，首先我们加载R包，读取all.events.ioe文件，并对数据框进行处理：

# 加载tidyverse库，这是一个包含了多个数据处理和可视化工具的R包集合  
library(tidyverse)

# 从指定路径读取一个以制表符分隔的文件，并设置header为True表示文件的第一行是列名  
df <- read.delim("all.events.ioe", header =T)

# 使用管道操作符对df进行处理  
# 1. 使用separate函数将event_id列拆分为gID和event_message两列，分隔符是";"  
# 2. 使用mutate函数和正则表达式从event_message中提取事件类别到新的列event_class  
# 3. 使用select函数选择gene_id, event_class, alternative_transcripts三列  
df <- df %>% separate(col = event_id,into=c("gID","event_message"), sep =";")%>%
mutate(event_class = str_extract(event_message,"^[^:]+"))%>%
select(gene_id, event_class, alternative_transcripts)

统计可变剪接影响了多少基因：

# 打印出发生可变剪接的基因数量  
print(paste("Number of AS genes is", length(unique(df$gene_id)))) 

统计由于可变剪接产生了多少转录异构体：

# 将df中的alternative_transcripts列转换为数据框，并重命名列为asline  
# 然后使用separate_rows函数将asline列中的每一行按","分隔成多行  
df2 <-as.data.frame(df$alternative_transcripts)%>%
`colnames<-`("asline")%>%
  separate_rows(asline, sep =",")

# 打印出经过可变剪接的转录本数量  
print(paste("Number of AS transcripts is",length(unique(df2$asline)))) 

统计发生两种及以上AS事件的基因数量：

# 对df进行处理，按gene_id分组，然后统计每个基因的不同事件类别数量  
# 过滤出那些具有至少两种不同事件类别的基因  
multi_group_genes <- df %>%
  group_by(gene_id)%>%
  summarise(group_count = n_distinct(event_class))%>%
  filter(group_count >=2)

# 打印出具有多种组合可变剪接事件的基因数量  
print(paste("genes had isoforms resulting from multiple combinatory AS events: ", nrow(multi_group_genes)))  

统计不同的可变剪接事件类别的频数和百分比：

# 统计df中event_class的频数，并转换为数据框  
# 然后添加Freq的整数版本和百分比版本  
classdf <- as.data.frame(table(df$event_class))  
classdf <- classdf %>% mutate(Freq = as.integer(Freq)) %>% mutate(percent = Freq / sum(classdf$Freq))  

# 打印出事件类别的频数和百分比  
print(classdf)

可视化不同的可变剪接事件类别的频数和百分比：

# 使用ggplot2绘制事件类别的柱状图  
ggplot(classdf, aes(Var1,Freq/1000, fill =Var1))+
  geom_bar(stat ='identity', width =0.8)+
  scale_fill_manual(values =c("#88d2d2","#bae6fa","#4c98ce","#c4c4d3","#ffef96","#eb9d9e","#fbcac5"))+
  ylab(expression("No. of AS events (10"^"4)"))+
  xlab("")+
  labs(fill="")+
  theme_classic()

完整代码

# 加载tidyverse库，这是一个包含了多个数据处理和可视化工具的R包集合  
library(tidyverse)

# 从指定路径读取一个以制表符分隔的文件，并设置header为True表示文件的第一行是列名  
df <- read.delim("all.events.ioe", header =T)

# 使用管道操作符对df进行处理  
# 1. 使用separate函数将event_id列拆分为gID和event_message两列，分隔符是";"  
# 2. 使用mutate函数和正则表达式从event_message中提取事件类别到新的列event_class  
# 3. 使用select函数选择gene_id, event_class, alternative_transcripts三列  
df <- df %>% separate(col = event_id,into=c("gID","event_message"), sep =";")%>%
  mutate(event_class = str_extract(event_message,"^[^:]+"))%>%
select(gene_id, event_class, alternative_transcripts)

# 打印出发生可变剪接的基因数量  
print(paste("Number of AS genes is",length(unique(df$gene_id))))

# 将df中的alternative_transcripts列转换为数据框，并重命名列为asline  
# 然后使用separate_rows函数将asline列中的每一行按","分隔成多行  
df2 <-as.data.frame(df$alternative_transcripts)%>%
`colnames<-`("asline")%>%
  separate_rows(asline, sep =",")

# 打印出经过可变剪接的转录本数量  
print(paste("Number of AS transcripts is",length(unique(df2$asline))))

# 对df进行处理，按gene_id分组，然后统计每个基因的不同事件类别数量  
# 过滤出那些具有至少两种不同事件类别的基因  
multi_group_genes <- df %>%
  group_by(gene_id)%>%
  summarise(group_count = n_distinct(event_class))%>%
  filter(group_count >=2)

# 打印出具有多种组合可变剪接事件的基因数量  
print(paste("genes had isoforms resulting from multiple combinatory AS events: ", nrow(multi_group_genes)))

# 统计df中event_class的频数，并转换为数据框  
# 然后添加Freq的整数版本和百分比版本  
classdf <-as.data.frame(table(df$event_class))
classdf <- classdf %>% mutate(Freq=as.integer(Freq))%>% mutate(percent =Freq/sum(classdf$Freq))

# 打印出事件类别的频数和百分比  
print(classdf)

# 使用ggplot2绘制事件类别的柱状图  
ggplot(classdf, aes(Var1,Freq/1000, fill =Var1))+
  geom_bar(stat ='identity', width =0.8)+
  scale_fill_manual(values =c("#88d2d2","#bae6fa","#4c98ce","#c4c4d3","#ffef96","#eb9d9e","#fbcac5"))+
  ylab(expression("No. of AS events (10"^"4)"))+
  xlab("")+
  labs(fill="")+
  theme_classic()

写在后面

同门说，他总觉得这个年纪，抬头看见的不应该是水泥天花板，而应该是满天的星空。或许科研欺骗了他，生活欺骗了他，但他还是想，为自己，努把力，别自己欺骗自己。大概，是能有那么一天吧，他会砸开头顶的天花板，不论外面是否有漫天的星空。