附近含有 m6A 修饰的 Stop Codon 序列提取

背景

“

m6A 修饰存在于 RNA 上，在 RNA 的出生，转运和成熟、降解、翻译等过程发挥着重要的作用。

”

m6A 修饰分布于 RNA 上的 5'UTR、CDS、3'UTR、Exon 和 Intron 里，主要富集在 stop codon 附近，那么 stop codon 附近的 m6A 修饰具体发挥着怎样的功能和分子机制研究的不是很清楚，在不同物种间的、不同组织间和不同生理状态下 m6A 修饰的水平也会随之动态变化。

AI 画个示意图如下：

m6A 位点的预测和实验研究验证的信息都有很多在线数据库可以查询，对于数据库里数据的分析可以是多样化和个性化的，讲述一个合理的生物学过程，结合实验验证，是文章的主题框架和思路。

人生和生命科学一样，在不断折腾和求真知的路途中找到真理和生命的本质！

准备

整合实验数据的 m6A 修饰位点注释信息的数据库有好几个，如 RMBase、RNAmod、REPIC 和 m6A-atals 等，可以查询单个基因的 m6A 位点信息或者直接下载全部的信息。此外还可查询下载其它 RNA 修饰的信息。

1、数据下载
直接去数据库下载，我下载了 m6A-atals 数据库的 Human 的 m6A 的 Basic_Site_Information 信息，下载方法是进入网址后点击 Download 进入下载界面，下载相应文件即可。

• http://180.208.58.66/m6A-Atlas/index.html

或者 linux 里下载

$ wget http://180.208.58.66/m6A-Atlas/Download/m6A_Human_Basic_Site_Information.txt

下载好后的文件打开长这样：
包括染色体、起始位置、终止位置、分布区域、数据来源、实验类型等等

$  less -S m6A_Human_Basic_Site_Information.txt
        chr     start   end     pos     strand  name    gene_id dist    type    seq                     data    exper   cell
human_m6A_1     chr1    564474  564474  1       +       LOC101928626    ENSG00000225972 ncRNA_exonic    unprocessed_pseudogene  CAAGGTCAAAGGGAGTCCGAACTAGTCTCAGGCTTCAACAT     No      25491922        GSE54921        PA-m6A-seq      HeLa    Control
human_m6A_2     chr1    564545  564545  1       +       LOC101928626    ENSG00000225972 ncRNA_exonic    unprocessed_pseudogene  CTTCATAGCCGAATACACAAACATTATTATAATAAACACCC     No      30867593        GSE121942       miCLIP  HepG2   Control|25491922
human_m6A_3     chr1    564560  564560  1       +       LOC101928626    ENSG00000225972 ncRNA_exonic    unprocessed_pseudogene  CACAAACATTATTATAATAAACACCCTCACCACTACAATCT     No      30867593        GSE121942       miCLIP  HepG2   SETD2   knockdown|25491922

然后在 excel 里打开对数据进行筛选，假如筛选细胞系 cell 为 CD8T、type 为 protein coding、dist 为 exonic，intronic 和 3UTR，然后提取 chr、start、end、strand、name、dist 列：

$ less m6a.txt
chr     start   end     strand  name    dist
chr1    949587  949587  +       ISG15   exonic
chr1    1168512 1168512 +       B3GALT6 exonic
chr1    1169034 1169034 +       B3GALT6 UTR3

2、提取 stop codon 信息

要提取 stop codon 信息首先得知道数据库用的注释文件或基因组版本是哪有个，不然序列位置和注释信息是对不上的。去 m6A-atals 找了一下也没看到有写，那就先下一个 hg38 的 gtf 文件看看，结果发现序列位置差的太远了，肯定不是，于是又下载了 hg19 的才对的上，毕竟早期的文献数据分析都是用的较早的基因组和注释文件，hg19 基因组是 2009 年出来的。

下载 UCSC hg19 注释文件：

$ wget https://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/genes/hg38.ncbiRefSeq.gtf.gz

我们需要提取 stop codon 的位置信息：
chr、stop_codon、start、end、strand、gene_name 这么几列

$ zless -S hg19.ncbiRefSeq.gtf.gz | grep -w 'stop_codon' |sed 's/"/\t/g' |sed 's/;/\t/g' |awk '{print $1"\t"$3"\t"$4"\t"$5"\t"$7"\t"$18}' | wc -l
67419

有 67419 个 stop codon，一个基因有多个转录本，如果这多个转录本的 stop codon 位置不一样的话那么一个基因就会有多个不一样的 stop codon，但是这多个转录本的 stop codon 位置也可能是一样，所以这 67419 个终止密码子是有冗余或者重复的，我们需要去一下重复。

查看重复的 stop codon：

$ zless -S hg19.ncbiRefSeq.gtf.gz | grep -w 'stop_codon' |sed 's/"/\t/g' |sed 's/;/\t/g' |awk '{print $1"\t"$3"\t"$4"\t"$5"\t"$7"\t"$18}' |less
chr22   stop_codon      38689293        38689295        -       CSNK1E
chr22   stop_codon      38689293        38689295        -       CSNK1E
chr22   stop_codon      38689293        38689295        -       TPTEP2-CSNK1E
chr22   stop_codon      38617476        38617478        -       TMEM184B
chr22   stop_codon      38617476        38617478        -       TMEM184B
chr22   stop_codon      38617476        38617478        -       TMEM184B
chr22   stop_codon      38610883        38610885        +       MAFF
chr22   stop_codon      38610883        38610885        +       MAFF
chr22   stop_codon      38610883        38610885        +       MAFF
chr22   stop_codon      38610883        38610885        +       MAFF

去除重复：

$ zless -S hg19.ncbiRefSeq.gtf.gz | grep -w 'stop_codon' |sed 's/"/\t/g' |sed 's/;/\t/g' |awk '{print $1"\t"$3"\t"$4"\t"$5"\t"$7"\t"$18}' |sort -k3,3n |uniq |wc -l
27734
# 输出到文件保存结果
$ zless -S hg19.ncbiRefSeq.gtf.gz | grep -w 'stop_codon' |sed 's/"/\t/g' |sed 's/;/\t/g' |awk '{print $1"\t"$3"\t"$4"\t"$5"\t"$7"\t"$18}' |sort -k3,3n |uniq > stopc.txt

最后我们把 stopc.txt 第二列放到最后一列，给它加上和 m6a.txt 一样的列名就可以了。

$ less stopc.txt
chr     start   end     strand  name    dist
chr1    70006   70008   +       OR4F5   stop_codon
chr1    368595  368597  +       OR4F29  stop_codon
chr1    879531  879533  +       SAMD11  stop_codon

算法思路

计算每个 m6A 修饰位点距离 stop codon 的距离，满足一定距离（100bp）的保存对应的 stop codon 位置信息。

分析：
1、m6A 修饰位点可位于 stop codon 上游也能是下游
2、一个基因存在多个 stop codon 位点 （一样位置的stop codon已经去除）
3、一个基因存在多个 m6A 修饰位点

YTHDF3 有 8 个转录本：

要是一个个去 excel 里找到每个基因对应的 m6A 位点和 stop codon 位点再互相减，挑选出在 100bp 以内的 stop codon 位点信息。如果有多个 stop codon 得更加多次的相减，几万个基因得减到猴年马月！此时对于这样一次次条件判断的重复过程我们可以使用编程的循环语句就能轻松实现。

解决思路
1、循环每个基因，取出相应的 m6A 和 stop codon 信息
2、判断该基因的 stop codon 的数量
3、如果 stop codon ==1，则与每个 m6A 位置相减，取出绝对值小于 100bp 的 stop codon
3、如果 stop codon >1，循环该基因每个 stop codon，与每个 m6A 位置相减，取出绝对值小于 100bp 的 stop codon

和师妹共同讨论了文件和解决思路，一起完成了一个循环。以下是代码部分：

# 加载R包
library(tidyverse)

# 1、读取筛选过后的m6a文件和stop codon信息文件
stopcodon <- read.delim("C:\\Users\\admin\\Desktop\\stopc.txt")
m6A_293 <- read.delim("C:\\Users\\admin\\Desktop\\m6a.txt")

# 2、取共有基因交集
m6A_gene <- unique(m6A_293$name)
stop_gene <- unique(stopcodon$name)
oveloap <- intersect(m6A_gene,stop_gene)

# 3、去除没有m6a或者stop codon信息的基因
m6A_gene_clean <- m6A_293[which(m6A_293$name %in% oveloap ),]
stopcodon_clean <-stopcodon[which(stopcodon$name %in% oveloap),]

# 4、在stop codon clean信息文件添加一列，为终止密码子中间位置
stopcodon_clean$pos <- stopcodon_clean$start+1

# 5、循环筛选符合m6a位点在密码子上下游15nt 以内符合条件的stop codon信息

for(i in oveloap){
  temp_m6A <- m6A_gene_clean %>% filter(name == i)
  temp_stop <- stopcodon_clean %>% filter(name == i)
  if(length(temp_stop$pos) == 1){

    test_16 <-  abs(temp_m6A$start - temp_stop$pos) <= 100

    if("TRUE" %in% test_16){
      tar_stop1 <- temp_stop
      file1 <- paste(i,'s_tar_stop',sep="_")
      assign(file1,temp_stop)
    }else {
    }
  }else{
    for(m in temp_stop$pos){

      test_16m <-  abs(temp_m6A$start - m) <= 100
      # print(test_16m)
      # print(i)
      if("TRUE" %in% test_16m){
        tar_stop2 <- temp_stop %>% filter(pos == m)
        file2 <- paste(i,'m_tar_stop',sep="_")
        assign(file2,tar_stop2)
      }
      else{
      }
    }
  }
}

# 6、查看含有单个密码子和多个密码子的筛选结果数量
length(ls(pattern = '*s_tar_stop'))
[1] 725
length(ls(pattern = '*m_tar_stop'))
[1] 279

# 7、储存总的结果变量名
res_name <- ls(pattern ='*_[s,m]_tar_stop' )

# 8、结果数量
length(res_name)
[1] 1004

# 9、合并所有结果到一个数据表里
nm <- map(res_name,as.name)
final_res <- do.call('rbind',nm)
# 10、整理成bed格式的文件便于提取序列
final_res_bed <- final_res[,c(1,2,3,5,6,4,7)]

# 11、提取stop codon序列的文件
final_res_bed2fa <- final_res_bed
final_res_bed2fa$start <- final_res_bed2fa$start-1
head(final_res_bed)
    chr     start       end    name       dist strand       pos
1 chr10 101611386 101611388   ABCC2 stop_codon      + 101611387
2 chr21  43716500  43716502   ABCG1 stop_codon      +  43716501
3  chr3  52015032  52015034 ABHD14A stop_codon      +  52015033
4  chr6  24701841  24701843  ACOT13 stop_codon      +  24701842
5  chr4   2930248   2930250    ADD1 stop_codon      +   2930249
6  chr2 228419209 228419211   AGFG1 stop_codon      + 228419210

# 12、保存结果文件
write.table(final_res,file="C:\\Users\\admin\\Desktop\\final_res.csv",quote = F,sep=",",row.names = F)
write.table(final_res_bed2fa,file="C:\\Users\\admin\\Desktop\\final_res_bed2fa.csv",quote = F,sep=",",row.names = F)

这里产生的结果文件会以单个文件保存在环境变量中，环境变量会生成很多文件名。

代码优化：

想了一下其实不用分 stop codon 的数量，不管有几个，直接与每个 m6A 位置相减，取出绝对值小于 100bp 的 stop codon 就行了，其次把每次循环产生的结果保存追加到一个 list 中，这样环境变量就不会生成很多文件了。

以下是优化后的代码部分，循环前的代码是一样的，就不放了，这里放循环及后面代码：