转录组那些事儿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介绍了