• 主页
  • 课程

    关于课程

    • 课程归档
    • 成为一名讲师
    • 讲师信息
    alphafold2培训

    alphafold2培训

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

      关于课程

      • 课程归档
      • 成为一名讲师
      • 讲师信息
      alphafold2培训

      alphafold2培训

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

      老俊俊的生信笔记

      • 首页
      • 博客
      • 老俊俊的生信笔记
      • 跟着 Cell Reports 学做图-CLIP-seq 数据可视化

      跟着 Cell Reports 学做图-CLIP-seq 数据可视化

      • 发布者 weinfoadmin
      • 分类 老俊俊的生信笔记
      • 日期 2022年3月24日
      测试开头

      跟着 Cell Reports 学做图-CLIP-seq 数据可视化
      点击上方,关注老俊俊!

      1引言

      承接上一篇推文,使用写的函数来复现一篇文章的结果。

      m6A enrichment on peak center

      文章标题:

      HITS-CLIP and Integrative Modeling Define the Rbfox Splicing-Regulatory Network Linked to Brain Development and Autism

      一篇关于 HITS-CLIP-seq 的文章,揭示了 Rbfox2 这个 RBP(RNA Binding Protein) 调控剪接的机制。

      原理:

      跟着 Cell Reports 学做图-CLIP-seq 数据可视化
      • 通过 紫外交联 RBP 于 RNA 序列的结合,然后通过 蛋白酶的消化 获取结合的 RNA 序列,这样处理过的序列在后面反转 cDNA 序列时会引起序列的 突变 或者 截断 。
      • 后面数据分析则是 寻找这些突变或者截断的位点 来准确识别 RBP 所结合的基序及位置,从而达到单碱基分辨率的效果。

      复现图如下:

      跟着 Cell Reports 学做图-CLIP-seq 数据可视化

      图注:

      跟着 Cell Reports 学做图-CLIP-seq 数据可视化

      2计算相对距离

      计算 RBFOX2 的 motif 的 depletion/truncation 相对距离,详见上篇推文:

      import re
      import pandas as pd
      from collections import Counter
      import seaborn as sns
      import matplotlib.pyplot as plt
      from pyfaidx import Fasta

      定义函数:

      def GetPeaksFa(peak,genomefie,outfasta,type='bed',a=0,b=0,minlen=5):
          # load geneome
          genome = Fasta(genomefie)
          outfa = open(outfasta,'w')
          # process peaks file
          with open(peak,'r') as bed:
              for line in bed:
                  fields = line.replace('n','').split()
                  chr = fields[0]
                  strand = fields[5]
                  start = int(fields[1])
                  end = int(fields[2])
                  if type == 'bed':
                      end = end
                  elif type == 'narrowpeak':
                      end = end + 1
                  else:
                      print('please mark your peaks file type(bed/narrowpeak)')
                      break
                  # extend upstream and downstram peaks
                  if strand == '+':
                      seq = genome[chr][(start - a):(end + b)].seq
                  elif strand == '-':
                      seq = genome[chr][(start + a):(end -b)].complement.reverse.seq
                  else:
                      seq = genome[chr][(start - a):(end + b)].seq
                  # mimimum seq length
                  if len(seq) >= minlen:
                      # fa name
                      name = fields[3] + '|' + '|'.join([chr,str(start),str(end)])
                      # output seq
                      outfa.write('>' + name + '|' + strand + 'n')
                      outfa.write(seq + 'n')


          outfa.close()

      ###############################################################################################

      def CalculatePeaksRelativeDist(peaksfa,motif='[G|A][G|A]AC[A|C|T]|[T|G|A]GT[C|T][C|T]',mid=3):
          # save in list
          relpos = []
          # open fa
          with open(peaksfa,'r') as seq:
              for line in seq:
                  line = line.replace('n','')
                  if line.startswith('>'):
                      next
                  peakmid = round((len(line) + 1)/2,ndigits=0)
                  # search all motif
                  pattern = re.compile(motif,flags=re.IGNORECASE)
                  pos_res = pattern.finditer(line)
                  # shift to A site position
                  allpos = [pos.start() + mid for pos in pos_res]
                  # calculate relative to peak center bases
                  relposition = [pos - peakmid for pos in allpos]
                  # relposition = [peakmid - pos for pos in allpos]
                  # save in list
                  relpos.extend(relposition)
          return relpos

      提取序列

      # 提取序列
      GetPeaksFa(peak='RBFOX2-CIMS.bed',
                 genomefie='mm10.fa',type='bed',
                 a=50,b=50,
                 outfasta='RBFOX2-CIMS-peaks_sequnce.fa')

      GetPeaksFa(peak='RBFOX2-CITS.bed',
                 genomefie='mm10.fa',type='bed',
                 a=50,b=50,
                 outfasta='RBFOX2-CITS-peaks_sequnce.fa')

      计算距离

      # 计算距离
      CIMS = CalculatePeaksRelativeDist(peaksfa='RBFOX2-CIMS-peaks_sequnce.fa',
                                          motif='TGCATG',mid=1)

      CITS = CalculatePeaksRelativeDist(peaksfa='RBFOX2-CITS-peaks_sequnce.fa',
                                          motif='TGCATG',mid=1)

      保存数据

      这里保存数据可以导入 R 里面进行可视化。

      先定义一个保存函数:

      # save data function
      def saveData(listres,outputname):
          tmp_freq = Counter(listres).most_common()
          # sort ascending
          tmp_freq.sort(key= lambda item:item[0])
          # to data frame and save
          tmp_freq_df = pd.DataFrame(tmp_freq,columns=['reldist','frequnecy'])
          tmp_freq_df.to_csv(outputname,header=True,index=False)

      保存:

      saveData(listres=CIMS,outputname='CIMS_freq_df.csv')
      saveData(listres=CITS,outputname='CITS_freq_df.csv')

      3python 可视化

      接下来对分析的数据进行绘图:

      plt.figure(figsize=[16,5])
      sns.set(style='ticks')

      plt.subplot(1,2,1)
      # label size
      plt.tick_params(labelsize=18)
      # density plot
      sns.kdeplot(data=CIMS,shade=False,color='#398AB9',alpha=1,bw_adjust=0.09,linewidth=4)
      # labels
      plt.xlabel("Distance to deletion site",fontsize=20)
      plt.ylabel("Density of UGCAUG",fontsize=20)
      # x lines
      plt.axvline(x=0,ls="--",lw=2,c='black')
      plt.axvline(x=-5,ls="--",lw=2,c='grey')
      plt.xlim(-50,50)
      plt.xticks(ticks=[-50,-30,-10,0,10,30,50],labels=[-50,-30,-10,0,10,30,50])
      plt.legend(['UGCAUG'],frameon=False,fontsize = 16,loc='upper left')

      plt.subplot(1,2,2)
      # label size
      plt.tick_params(labelsize=18)
      # density plot
      sns.kdeplot(data=CIMS,shade=False,color='#398AB9',alpha=1,bw_adjust=0.09,linewidth=4)
      # labels
      plt.xlabel("Distance to deletion site",fontsize=20)
      plt.ylabel("Density of UGCAUG",fontsize=20)
      # x lines
      plt.axvline(x=0,ls="--",lw=2,c='black')
      plt.axvline(x=-5,ls="--",lw=2,c='grey')
      plt.xlim(-6,6)
      # plt.xticks(ticks=[-6,-5,0,5,6],labels=[-6,-5,0,5,6])
      plt.legend(['UGCAUG'],frameon=False,fontsize = 16,loc='upper right')

      # save
      plt.savefig('Distance to deletion site.pdf',format='pdf')
      跟着 Cell Reports 学做图-CLIP-seq 数据可视化

      左边是主图,右边是放大的小图。


      绘制 truncation 的图:

      plt.figure(figsize=[16,5])
      sns.set(style='ticks')

      plt.subplot(1,2,1)
      # label size
      plt.tick_params(labelsize=18)
      # density plot
      sns.kdeplot(data=CITS,shade=False,color='#398AB9',alpha=1,bw_adjust=0.09,linewidth=4)
      # labels
      plt.xlabel("Distance to truncation site",fontsize=20)
      plt.ylabel("Density of UGCAUG",fontsize=20)
      # x lines
      plt.axvline(x=0,ls="--",lw=2,c='black')
      plt.axvline(x=-5,ls="--",lw=2,c='grey')
      plt.xlim(-50,50)
      plt.xticks(ticks=[-50,-30,-10,0,10,30,50],labels=[-50,-30,-10,0,10,30,50])
      plt.legend(['UGCAUG'],frameon=False,fontsize = 16,loc='upper left')

      plt.subplot(1,2,2)
      # label size
      plt.tick_params(labelsize=18)
      # density plot
      sns.kdeplot(data=CITS,shade=False,color='#398AB9',alpha=1,bw_adjust=0.09,linewidth=4)
      # labels
      plt.xlabel("Distance to truncation site",fontsize=20)
      plt.ylabel("Density of UGCAUG",fontsize=20)
      # x lines
      plt.axvline(x=0,ls="--",lw=2,c='black')
      plt.axvline(x=-5,ls="--",lw=2,c='grey')
      plt.xlim(-6,6)
      # plt.xticks(ticks=[-6,-5,0,5,6],labels=[-6,-5,0,5,6])
      plt.legend(['UGCAUG'],frameon=False,fontsize = 16,loc='upper right')

      # save
      plt.savefig('Distance to truncation site.pdf',format='pdf')
      跟着 Cell Reports 学做图-CLIP-seq 数据可视化

      4计算离 VGCAUG 的距离画图

      V 为 A/C/G 三个碱基:

      CIMSv = CalculatePeaksRelativeDist(peaksfa='RBFOX2-CIMS-peaks_sequnce.fa',
                                          motif='[A|C|G]GCATG',mid=1)

      CITSv = CalculatePeaksRelativeDist(peaksfa='RBFOX2-CITS-peaks_sequnce.fa',
                                          motif='[A|C|G]GCATG',mid=1)

      saveData(listres=CIMSv,outputname='CIMSv_freq_df.csv')
      saveData(listres=CITSv,outputname='CITSv_freq_df.csv')

      画图:

      plt.figure(figsize=[15,5])
      sns.set(style='ticks')

      plt.subplot(1,2,1)
      # label size
      plt.tick_params(labelsize=18)
      # density plot
      sns.kdeplot(data=CIMSv,shade=False,color='#398AB9',alpha=1,bw_adjust=0.09,linewidth=4)
      # labels
      plt.xlabel("Distance to deletion site",fontsize=20)
      plt.ylabel("Density of VGCAUG",fontsize=20)
      # x lines
      plt.axvline(x=0,ls="--",lw=2,c='black')
      plt.axvline(x=-5,ls="--",lw=2,c='grey')
      plt.xlim(-50,50)
      plt.xticks(ticks=[-50,-30,-10,0,10,30,50],labels=[-50,-30,-10,0,10,30,50])
      plt.legend(['VGCAUG'],frameon=False,fontsize = 16,loc='upper left')

      plt.subplot(1,2,2)
      # label size
      plt.tick_params(labelsize=18)
      # density plot
      sns.kdeplot(data=CIMSv,shade=False,color='#398AB9',alpha=1,bw_adjust=0.09,linewidth=4)
      # labels
      plt.xlabel("Distance to deletion site",fontsize=20)
      plt.ylabel("Density of VGCAUG",fontsize=20)
      # x lines
      plt.axvline(x=0,ls="--",lw=2,c='black')
      plt.axvline(x=-5,ls="--",lw=2,c='grey')
      plt.xlim(-6,6)
      # plt.xticks(ticks=[-6,-5,0,5,6],labels=[-6,-5,0,5,6])
      plt.legend(['VGCAUG'],frameon=False,fontsize = 16,loc='upper right')
      跟着 Cell Reports 学做图-CLIP-seq 数据可视化

      后面代码类似,这里只放图:

      跟着 Cell Reports 学做图-CLIP-seq 数据可视化

      5导入 R 进行可视化

      测试结尾

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

      • 分享:
      作者头像
      weinfoadmin

      上一篇文章

      m6A enrichment on peak center
      2022年3月24日

      下一篇文章

      关于序列提取的问题思考
      2022年3月26日

      你可能也喜欢

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

      搜索

      分类

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

      投稿培训

      免费

      alphafold2培训

      免费

      群晖配置培训

      免费

      最新博文

      白介素-17受体信号的自主激活,维持炎症并促进疾病进展
      048月2023
      MCT4依赖的乳酸分泌抑制LKB1缺陷肺腺癌的抗肿瘤免疫
      187月2023
      Nature | 单细胞技术揭示衰老细胞与肌肉再生
      301月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年
      在线支付 激活码

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