转录组那些事儿Part II

今天是生信星球陪你的第103天

你想找辆共享单车,发现满街都是别家车,没有一辆你能骑。

你想学点生信,搜了“初学者教程”,满眼尽是高大上,没有一句能看懂。

终于你跨越茫茫宇宙,来到生信星球,发现了初学者的新大陆

.阅读文献

选取的文献是哈佛大学Lee L. Rubin团队于2015年发表在Cell Stem Cell中的Genome-wide RNA-Seq of Human Motor Neurons Implicates Selective ER Stress Activation in Spinal Muscular Atrophy,doi是:http://dx.doi.org/10.1016/j.stem.2015.08.003

文章梗概

运动神经元存活蛋白(Survival Motor Neuron,SMN)是人体内的泛表达蛋白,它的会导致脊髓性肌萎速症(Spinal muscular atrophy,SMA),SMA是由SMN1基因突变引起的,但是这个基因是广泛表达的基因,为什么运动神经元(Motor Neurons,MNs)偏偏是最易受影响的细胞类型之一?实验通过对照组和SMA患者诱导多能干细胞(IPSC)从而产生纯化的运动神经元,通过固定、抗体标记,然后进行了RNA测序研究。在患者运动神经元中发现了SMA特异性的变化,其中包括内质网(ER)应激通路的过度激活。功能研究表明,抑制内质网的应激反应能提高运动神经元的存活率。在患病小鼠中,使用ER应激抑制剂穿过血脑屏障可以保护脊髓运动神经元。因此,选择性激活内质网应激反应,导致了SMA患病个体的运动神经元死亡

实验流程

RNA-seq数据存储在GSE69175中

GEO数据

.开始实战

万事之源,软件为先

配置好工作路径

下载测试数据

GEO数据库中的编号是:GSE69175

关于数据:

NGS测序数据一般会上传到SRA数据库里面,但是怎么从GEO数据库中找到SRA原始数据的下载地址?【GEO数据库地址:https://www.ncbi.nlm.nih.gov/geo/】

具体的层级关系是:SRP(项目)—>SRS(样本)—>SRX(数据产生)—>SRR(数据本身)

SRR数据库地址:https://ftp-private.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/

质控过滤

这里使用的是fastp,同时融合了质控、过滤的模式;

同时也可以使用fastqc + trimmomatic/ trim_galore进行

下载参考基因组和注释文件

人类参考基因组选择UCSC数据库;注释文件选择GENCODE,https://www.gencodegenes.org/

下载好基因组,需要构建基因组索引

如果自己的项目做非模式物种,可以用二代三代测序数据组装成参考转录组,例如trinity组装。如果做的物种有参考基因组,就方便一些,可以直接从相关数据库中下载参考基因组

比对

使用了samtools的三件套:转换(view)、排序(sort)、建索引(index)

转换: -S指输入文件格式(不加-S默认输入是bam),-b指定输出文件(默认输出sam)【如果要bam转sam,-h设置输出sam时带上头注释信息】

排序:对bam排序,-@ 线程数, -o输出文件

索引:结果会产生.bai文件【必须排序后才能建索引~就像体育课必须从高到矮排好以后再报数】

基本信息统计

reads计数

基于基因组水平,可以用Htseq-count和featureCounts

Htseq-count:它是用python写的一个脚本,安装好以后可以直接拿来用

需要比对文件sam/bam格式;需要基因组注释文件gff/gtf格式。这两个文件的染色体名称必须相同,因为需要将比对位置和特征信息相匹配就要借助染色体名。

它的工作原理是:先通过bam文件找到reads比对上的外显子,然后去gtf文件中将外显子的基因ID进行统计,得到的结果就是未经标准化的表达量

htseq-count的是strand意思,默认使用链特异性建库方式,也就是计算过程中要求read1只能和双端测序的其中一条比较,read2只能和另一条比较;我们一般设置no即可,表示每条read都和正负链比较一次;

两个选择 ,数据根据位置或者read名称排序,默认name;

输入文件格式bam/sam;

一般也就是默认union模式

结果每个样本输出一个count文件,其中包含了基因名和reads计数;另外,如果你看到文件倒数几行,会发现还有几行带文字的

no_feature:比对区域与任何基因都没有重叠

ambiguous:比对区域与多个基因都发生重叠

too_low_aQual:比对质量低于设定阈值(默认是10)

not_aligned:无法比对上

alignment_not_unique:比对位置不唯一

featureCounts:【相比于htseq更快】

设置线程,表示针对双端测序数据,默认将外显子认定为基因组的一个feature(搜索特征), 默认按照基因名来计数, 输入注释文件, 输出文件

结果有两个,一个count文件,一个summary文件

对两个软件的结果进行合并

统计一下两个软件的计数之和

本来想两部分,可是自己写脚本、优化需要大量的时间,关于R处理的部分

只能放在Part III介绍了

  • 发表于:
  • 原文链接https://kuaibao.qq.com/s/20180821G1UNAD00?refer=cp_1026
  • 腾讯「云+社区」是腾讯内容开放平台帐号(企鹅号)传播渠道之一,根据《腾讯内容开放平台服务协议》转载发布内容。

扫码关注云+社区

领取腾讯云代金券

年度创作总结 领取年终奖励