• 主页
  • 课程

    关于课程

    • 课程归档
    • 成为一名讲师
    • 讲师信息
    同等学历教学

    同等学历教学

    免费
    阅读更多
  • 特色
    • 展示
    • 关于我们
    • 问答
  • 事件
  • 个性化
  • 博客
  • 联系
  • 站点资源
    有任何问题吗?
    (00) 123 456 789
    weinfoadmin@weinformatics.cn
    注册登录
    恒诺新知
    • 主页
    • 课程

      关于课程

      • 课程归档
      • 成为一名讲师
      • 讲师信息
      同等学历教学

      同等学历教学

      免费
      阅读更多
    • 特色
      • 展示
      • 关于我们
      • 问答
    • 事件
    • 个性化
    • 博客
    • 联系
    • 站点资源

      老俊俊的生信笔记

      • 首页
      • 博客
      • 老俊俊的生信笔记
      • 附近含有 m6A 修饰的 Stop Codon 序列提取

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

      • 发布者 weinfoadmin
      • 分类 老俊俊的生信笔记
      • 日期 2021年9月10日
      • 评论 0评论

      感谢老俊俊的大力支持。我们会每日跟新,欢迎您关注老俊俊的生信笔记。

      附近含有 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:UsersadminDesktopstopc.txt")
      m6A_293 <- read.delim("C:UsersadminDesktopm6a.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:UsersadminDesktopfinal_res.csv",quote = F,sep=",",row.names = F)
      write.table(final_res_bed2fa,file="C:UsersadminDesktopfinal_res_bed2fa.csv",quote = F,sep=",",row.names = F)

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

      代码优化:

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

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

      请关注“恒诺新知”微信公众号,感谢“R语言“,”数据那些事儿“,”老俊俊的生信笔记“,”冷🈚️思“,“珞珈R”,“生信星球”的支持!

      • 分享:
      作者头像
      weinfoadmin

      上一篇文章

      提取代表转录本之 gencode
      2021年9月10日

      下一篇文章

      Introduction of m6A
      2021年9月10日

      你可能也喜欢

      8-1651542331
      跟着Nature学绘图(2) 箱线图-累积分布曲线图
      2 5月, 2022
      9-1651542322
      Julia 笔记之字符串
      2 5月, 2022
      0-1651542343
      Julia 笔记之数学运算和初等函数
      1 5月, 2022

      留言 取消回复

      要发表评论,您必须先登录。

      搜索

      分类

      • R语言
      • TCGA数据挖掘
      • 单细胞RNA-seq测序
      • 在线会议直播预告与回放
      • 数据分析那些事儿分类
      • 未分类
      • 生信星球
      • 老俊俊的生信笔记

      投稿培训

      免费

      alphafold2培训

      免费

      群晖配置培训

      免费

      最新博文

      Nature | 单细胞技术揭示衰老细胞与肌肉再生
      301月2023
      lncRNA和miRNA生信分析系列讲座免费视频课和课件资源包,干货满满
      301月2023
      如何快速批量修改 Git 提交记录中的用户信息
      261月2023
      logo-eduma-the-best-lms-wordpress-theme

      (00) 123 456 789

      weinfoadmin@weinformatics.cn

      恒诺新知

      • 关于我们
      • 博客
      • 联系
      • 成为一名讲师

      链接

      • 课程
      • 事件
      • 展示
      • 问答

      支持

      • 文档
      • 论坛
      • 语言包
      • 发行状态

      推荐

      • iHub汉语代码托管
      • iLAB耗材管理
      • WooCommerce
      • 丁香园论坛

      weinformatics 即 恒诺新知。ICP备案号:粤ICP备19129767号

      • 关于我们
      • 博客
      • 联系
      • 成为一名讲师

      要成为一名讲师吗?

      加入数以千计的演讲者获得100%课时费!

      现在开始

      用你的站点账户登录

      忘记密码?

      还不是会员? 现在注册

      注册新帐户

      已经拥有注册账户? 现在登录

      close
      会员购买 你还没有登录,请先登录
      • ¥99 VIP-1个月
      • ¥199 VIP-半年
      • ¥299 VIP-1年
      在线支付 激活码

      立即支付
      支付宝
      微信支付
      请使用 支付宝 或 微信 扫码支付
      登录
      注册|忘记密码?