一种非主流数据类型的分析之路(一)

-回复 -浏览
楼主 2019-08-12 16:58:53
举报 只看此人 收藏本贴 楼主


最近笔者接到一种较为特殊的数据类型分析,区别于传统的免疫沉淀--二代测序套路,这种数据比较类似于重亚硫酸盐处理的甲基化数据(RRBS,WGBS)处理。即:传统的实验是通过序列的有无来判定富集信息的,这种实验则要通过序列的变化来判断。

具体情况是:在正常的高等生物细胞内存在着RNA编辑(RNA-editing)这一生物学现象。这种所谓的编辑是发生在转录后的,通过将mRNA中的A替换成为I(对,你没看错,是I!当年本科生化怒拿及格的小编自然是搞不清楚I是什么奇葩碱基的...)来实现遗传信息的变异的。好在这个I化学性质类似G(这个我知道,叫腺嘌呤),因此三联密码子只是从??A等价变化成??G。由于二代测序实际上就是在做pcr,也服从碱基互补原理,因此测序结果上的变化与生物学现象也还是一致的。【有兴趣同学可以参看这篇文章:Nascent-Seq Indicates Widespread Cotranscriptional RNA Editing in Drosophila】

这种生物学现象很快被聪明的科学家拿来做人工改造(想想牛逼的CRISPR-Cas9),在2016年就有科学家将这样一个RNA-editing酶和RNA上的结合蛋白融合在一起。这样有结合蛋白的地方就很容易被顺带的RNA-editing酶所催化反应。这样改变的序列往往预示着被携带的结合蛋白所结合。【本项目首先要重复练手的也就是这篇文章:TRIBE : Hijacking an RNA-editing enzyme to identify cell-specific targets of RNA-binding proteins】

可能有同学会问,绕了一大圈,这不就是RIP(RNA-immunoprepitation)能解决的问题嘛。不错,理论上来说RIP也是可以解决的。不过免疫沉淀实验一大弊病就是看抗体的脸色,而且过程中会以大量核酸的损失为代价(这也是为什么近些年很多人着手开发low-input乃至single-cell的技术)。此外,免疫沉淀实验还有一个通病就是精度不够,这里作者上了一张与免疫沉淀的对比图,对比起来的确立竿见影。阳性就是阳性,没有call-peak这种剪不断理还乱的痛苦。

话多说一句,以前是有一种ChIP-exo的技术,可得到的sharp的免疫沉淀实验峰的。然而大家觉得普通的ChIP-seq够用了,于是这个技术就没市场了。不过从日常做数据分析的经验来看(区分不那么明显的结果占多数),这种实验其实潜力还很大。


(这是一条从闲谈进入分析模式的分界线。)

我们不妨试想一下,分析这样的数据会有哪些困难呢?

  1. 类似于RRBS数据,首先我们还是要带‘突变’进行比对。这要求我们要相比正常的比对要提高容错率。这种调参需要摸着石头过河,在比对率与真阳性之间做出平衡;

  2. 这种数据处理一般要排除SNP的影响。因此除了比对RNA外,我们还要同时比对品系本身的DNA(也就是筛出与标准基因组的差别部分);

  3. 暂时想不到更多了,开始干活吧....

数据下载

在接手真正数据前,我们需要重复TRIBE这篇文章的结果。作者已经给出了GSE号,只要跟着过去下就可以了。

不过最近NCBI的FTP服务器似乎是封闭了,传统的wget写循环下载的简单粗暴的方式似乎已经失效。不过官方的意思很明显了,我们规规矩矩用官方的SRA-toolkit吧。

先把SRR号全部提供好

接着用prefetch命令进行下载:

  1. nohup prefetch --option-file SraAccList.txt &

文件还挺大的,挂着回去睡觉吧...

重命名

一般的作者都会很良心的在目录页把SRR名与生物学意义名称联系起来,也就是这样的:

批量用mv命令那个痛快简直美滋滋。不过这次情况有点不一样,仔细看后缀为‘3’的样本似乎和‘5’的样本一模一样嘛。细胞系一样,蛋白一样,处理似乎也一样。这下尴尬了,这重命名之后也无法联系到完整的生物学意义了。于是我开始脑补,这作者是不是上传数据的时候才发现“Duplicated item illegal!”,最后干脆用excel表格拉了一个阶梯数列出来...

嘴上笑嘻嘻,心里mmp。我只能打开内容页接着看区别在哪里。里面写的是很清楚的

做了一个对照,有没有诱导蛋白表达。

这下就恶心了,我们当然是不可能(当然对有些人还真可能...)一个个戳开近40个样本去看详细内容嘛。于是没办法只能写爬虫解决,好在生物网站对于爬虫还是比较友好的:

注意到这种网站URL的规律性

接下来的的任务就交给python了:

  1. #!/usr/bin/env python3

  2. # -*- coding: utf-8 -*-

  3. """

  4. Created on Thu Jan  4 14:29:29 2018

  5. @author: ****

  6. """

  7. import urllib

  8. from urllib.request import urlopen

  9. import re

  10. import pandas as pd

  11. raw = pd.read_csv('/Users/*****/Documents/TRIBE.csv')

  12. raw = raw[['Run','SampleName','Name']]

  13. raw["Name"]=raw["Name"].str.replace(' ','_')

  14. dic = raw.to_dict('index')

  15. f = open('detail.txt','w')

  16. for i in dic:

  17.    name = dic[i]['SampleName']

  18.    url = 'https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc='+name

  19.    req = urllib.request.Request(url,headers = headers)

  20.    response = urlopen(req)

  21.    the_page = response.read()

  22.    the_page = the_page.decode('utf-8')

  23.    start = the_page.find('Characteristics')

  24.    end = the_page.find('Extracted molecule')    

  25.    seq = re.sub("<.+?>",'',the_page[start:end])

  26.    seq = seq.replace('celltribe','cell\ntribe')

  27.    print(dic[i]['Name']+'\n'+seq)

  28.    f.write(dic[i]['Name']+'\n'+seq)

  29. f.close()

结果如下:

看着信息刷刷刷的出来,那叫一个痛快~

今天先写到这里。小伙伴们如果想上车一起分析的,那就开始动手吧~~~不要掉队哦~

点击以下「关键词」,查看往期内容:

一个物种一个家

TCGA | 小工具 | 数据库 |组装注释 |   基因家族  |  Pvalue

基因预测  |bestorf |  sci NAR | 在线工具 | 生存分析 | 热图

 生信不死 | 初学者 | circRNA | 一箭画心| 十二生肖 | circos

 舞台|基因组 | 黄金测序 | 套路 | 杂谈组装 |  进化 | 测序简史


生信人


我要推荐
转发到