RNA-Seq数据使用方法

RNA-seq相关理论

概念
RNA-seq：比较不同条件或不同菌株间转录本转录差异。
获得转录本信息：

• 5’-TSS信息：dRNA-seq
• 全长转录本信息：SMART RNA-seq
• 转录本修饰信息：tag-seq
• 正在翻译的转录本信息：RNC-seq Ribo-seq
• 与蛋白相互作用转录本信息

文库构建
rRNAs depletion
探讨rRNA去除的必要性方法
分子数量tRNA较多，核糖体RNA(rRNA)总质量大，但其不活跃所以需要去除，mRNA不管是质量还是数量都是5%左右。

真核生物：55rRNAs 16rRNAs 23rRNAs
原核生物：55rRNAs 5.8rRNAs 18rRNAs 28 rRNAs

链特异性文库或非链特异性文库
链特异性转录组建库优点：

1. 基因定量更加准确
   通过reads定位到正确的链上(正义链)，可以分别得到mRNA和反义RNA的表达量，使基因表达量的计算更准确。只贴正义链，不贴反义链。
2. 适用于无参转录组组装
   链特异性文库通过区分正负链，避免有互补配对关系的编码与非编码转录本被组装成一条转录本，提升无参转录本组装的真实性。
3. 挖掘天然反义LncRNA
   很多LncRNA都来自反义RNA，如果使用普通转录组建库，较难以区分reads是来自mRNA，还是天然反义LncRNA。
   反义RNA是指与mRNA互补后，能抑制与疾病发生直接相关基因的表达的RNA。
4. 可变剪切事件检测更准确
   通过排除反义链上的转录本表达的干扰，链特异性文库可以有效降低可变剪接事件的假阳性率。比如由反向剪切形成的circRNA，在发枫新的circRNA时，链特异性的测序方式有助于提高CircRNA检测的淮确性。

背景
转录组：
广义上指某一生理条件下，细胞内所有转录产物的集合，包括信使mRNA、核糖体RNA、转运RNA及非编码RNA；
狭义上值所有mRNA的集合。

实验设计的注意事项

• 样本要有时空差异，至少设置两组样品；
• 每组样片止至少有三组生物学重复，统计学家认为重复数量越多越好，生物学来说3-5个就可；
• 测序深度(简单基因表达分析需要5M以上reads，小RNA至少30M)，常规3-5G测序基因质量；
• 文库构建(链特异性或非特异性)；
• 测序策略(单端或双端)；
• 测序平台(价格 读长 通量 准确率)；
• 确定分析流程(有参或无参 是否要分析选择性剪接 例如Htseq不能分析鉴定选择性剪接)

有参考组装基本流程与目标

1. 对照与处理，各取三个样品，液氮速冻，干冰保存送公司抽提RNA并测序，共6个RNA-seq测序样品；
2. 下机数据，fastqc检查质量，trimmomatics等去除低质量序列；
3. hista2/tophat2/bowtie2等软件，比对到参考基因组上，组装出转录本；
   NCBI下载该物种的参考基因组序列及注释，用bowtie2、hisat2等构建金银组索引。
4. htseq-count/cuffinks/stringtie等软件进行基因的表达量的统计，DESeq2/edger/RSEM/cuffdiff等进行基因差异的鉴定；
   RT-PCR，实时定量PCR等验证差异基因的表达。
5. 对差异基因进行聚类分析，包括GO、KEGG。(开展WGCNA共表达分析等，开展选择性剪接分析)

无参考基因组组装流程与目标

1. 对照与处理，各取三个样品，液氮速冻，干冰保存送公司抽提RNA并测序，共6个RNA-seq测序样品；
2. 下机数据，fastqc检查质量，trimmomatics等去除低质量序列；
3. trinity等软件组装出转录本，提取unigene;
4. RSEM等软件，进行差异表达基因的坚定；
5. 对差异基因进行聚类分析，包括GO、KEGG。RT-PCR验证部分基因，检验组装准确性。

文件格式

FASTQ文件

1. 文件用途：测序返回的一般数据格式，通常是压缩文件filename.fq.gz的格式。
2. 格式说明：
   fastq文件每四行代表一条序列
   第一行：记录序列测序时所用仪器以及在测序通道中坐标信息，以@开头；
   第二行：测序的序列信息，以ATCGN表示，由于荧光信号干扰无法判断是什么碱基是就用N表示；
   第三行：通常一个 号；
   第四行：与第二行碱基信息一个一个地对应，存储测序碱基的质量指(ASCⅡ字符表示)。
3. 查看方式
   zcat查看gzip压缩的文件；
   head -n 8显示前8行文件内容；

FASTA文件

1. 文件用途：fasta文件用于基因组或者基因的DNA或者蛋白质的序列信息存储；
2. 文件格式：
   以>符号开头，记录了该序列类型和所在基因组位置信息，也称"序列名字行"；
   序列行：一行或多行，为序列信息，soft-masked基因组会把所有重复区和低复杂区的序列用小写字母标出基因组，小写字母n标示未知碱基。
3. 不成文的小规范：
   第一部分是序列名字，与>相连；
   第二部分用空格与序列名字相连，表示注释信息，可以没有。

GFF/GTF格式

gff全称general feactureformat，主要用来注释基因组；
gtf全称general transfer format，主要用来对基因进行注释。
两种文件都有九列数据组成，前八列基本相同，最后一列有一些小差别：

• 根据所使用的软件不同，feature types是必须注明的；
• 第九列必须以gene_id以及transcript_id开头。标签与值以空格分开，键值有引号，且每个特征之后要有分号。

example：
gene_id “geneA”;transcript_id “geneA.1”

质量控制和预处理软件

质量问题通常来自测序本身或前面的文库制备。包括：

• 可信度低的碱基
• 序列特异性的偏差
• 3’/5’位置偏差
• 聚合酶链反应(PCR)假象
• 未修剪的接头
• 序列污染

通过过滤、修剪、纠错或偏差订正被矫正。

1FastQC

输入文件可以是FASTQ(未压缩或压缩的)或SAM/BAM文件。生成html的质量报告，包括：

• 读段的数目及质量编码
• 可视化有关碱基质量和内容
• 读取长度及k-mer内容
• 有含糊不清的碱基 过度代表的序列和重复的信息

1.2html结果文件解读

左侧会有Summary可以看到基本情况，绿色对号代表"PASS"，黄色感叹号代表"WARN"，红色叉号代表"FAIL"，当黄色出现的时候需要查看结果。

参考文章：FastQC结果解读
1.Basic Statistics

• Filename:文件名
• File type: 文件类型
• Encoding：测序平台的版本和相应的编码版本号 用于计算Phred反推error P时用
• Total Sequences: 输入文本的reads的数
• Sequence length: 测序长度
• %GC: GC含量 表示整体序列的GC含量 由于二代测序GC偏好性高 且深度越高 GC含量会越高

2.Per base sequence quality

横轴为read长度，纵轴为质量得分，柱状表示该位置所有序列的测序质量的统计，柱状(黄色)是25%-75%区间质量分布，error bar(触须)是10%-90%区间质量分布，蓝线代表平均数，红色代表中位数。

一般要求所有位置的10%小于20，即最多允许该位置10%的序列低于Q20，即90%的序列的碱基质量都大于Q20，即90%的序列碱基错误率不超过99%。当任何碱基质量低于10，或者任何中位数低于25时报WARN,需注意；当任何碱基质量低于5或者任何中位数低于20报FAIL。

碱基质量值Q越高表明碱基识别越可靠，准确度越高。Q20与Q30的含义：

• Q20为每100个碱基中会有一个识别错，即正确识别率为2个9，99%，当Phred = 20 时，碱基识别出错率为1/100，碱基识别正确率为99%，Q-score = -10 ✖ lg 10-2=20
• Q30为每1000个碱基中会有一个识别错，正确识别率为3个9，99.9%，当Phred = 30 时，碱基识别出错率为1/1000，碱基识别正确率为99.9%，Q-score = -10 ✖ lg 10-3=30

3.Per Sequence Quality Scores

• 每条reads的quality的均值的分布
• 横轴表示Q值，纵轴表示每个值对应的read数目，当测序结果主要集中在高分中，证明测序质量良好
• 当峰值小于27（错误率0.2%）时报"WARN"，当峰值小于20（错误率1%）时报"FAIL"

4.Per Base Sequence Content

• 对所有reads的每一个位置，统计ATCG四种碱基(正常情况)的分布，用于检查是否有AT,GC分离现象
• 横轴为碱基长度分布，纵轴表示百分比，图中4条线分别代表A，C，T，G在每个位置上的平均含量
  由于测序平台及测序长度不同，以及测序仪开始状态不稳定经常出现前后波动情况
  好的样本中四条线应该平行且接近。当部分位置碱基的比例出现bias时，即四条线在某些位置纷乱交织，往往提示我们有overrepresented sequence的污染。当所有位置的碱基比例一致的表现出bias时，即四条线平行但分开，往往代表文库有bias(建库过程或本身特点)，或者是测序中的系统误差。
  在碱基含量分布图，前几个碱基可能会出现较大波动，这是由于随机引物扩增偏差原因造成的
  当任一位置的A/T比例与G/C比例相差超过10%，报"WARN"；当任一位置的A/T比例与G/C比例相差超过20%，报"FAIL"

5.Per Sequence GC Content

• 统计reads的平均GC含量的分布
• 红线是实际情况，蓝线是理论分布(正态分布，均值不一定在50%，而是由平均GC含量推断的)
• 曲线形状的偏差往往是由于文库的污染或是部分reads构成的子集有偏差(overrepresented reads)。形状接近正态但偏离理论分布的情况提示我们可能有系统偏差
• 偏离理论分布的reads超过15%时，报"WARN"；偏离理论分布的reads超过30%时，报"FAIL"

6.Per base N content

• 出现测序仪不能分辨的碱基时会产生N，横轴为碱基分布，纵轴为N比率
• 当任一位置N的比率超过5%报WARN，超过20%报FAIL

7. Sequence Length Distribution

• reads长度的分布
• 理论上每次测序仪测出的read长度时一致的，但是由于建库等因素通常会导致一些小片段，reads长度不一致时报"WARN"；当有长度为0的read时报“FAIL

8. Sequence Duplication Levels

• 统计序列完全一致的reads频率，横坐标是duplication的次数，纵坐标是duplicated reads的数目，以unique reads的总数作为100%
• 一般测序深度越高，越容易产生一定程度的重复序列
• 当非unique的reads占总数的比例大于20%时，报"WARN"；当非unique的reads占总数的比例大于50%时，报"FAIL"

9.Overrepresented sequences

• 如果有某个序列大量出现，就叫做over-represented
• fastqc的标准是占全部reads的0.1%以上。和上面的duplicate analysis一样，为了计算方便，只取了fq数据的前200,000条reads进行统计，所以有可能over-represented reads不在里面。而且大于75bp的reads也是只取50bp。如果命令行中加入了-c contaminant file，出现的over-represented sequence会从contaminant_file里面找匹配的hit(至少20bp且最多一个mismatch)，可以给我们一些线索
• 当发现超过总reads数0.1%的reads时报"WARN"，当发现超过总reads数1%的reads时报"FAIL"

10.Adapter Content

• 横轴表示碱基位置，纵轴表示百分比
• 当fastqc分析时没有选择参数-a adapter list时，默认使用图例中的4种通用adapter序列进行统计
• 若有adapter残留，后续必须去接头

Bowtie2-HTseq

下载miniccoda
针对不考虑可变剪切的情况：原核生物转录组RNA-seq

1. 直接由Bowtie2比对到参考基因组，得到sam文件
2. santools将sam转为bam并排序
3. Htseq-count根据bam文件统计个基因counts

Htseq-Count运行

htseq-count -f bam/sam -r bam文件排序方式 -s 链特异性文库 -a 质量控制默认10 -n 线程数量 -t gene -i Name -m 计数模式 bam 文件路径 注释GFF/GTF文件>转录本counts文件

htseq-count -f bam -r name -s no -n 64 -t gene -i Name -m intersection-strict
SRR5176531.bam ../GCA 000009725.1_ASM972v1_genomic.gff>
SRR5176531.counts.txt

对结果进行合并操作

awk vim 去除最后5行 加表头

awk 'NR==FNR{a[$1]=$2}NR>FNR&&a[b=$1]{print $0,a[b]}'
SRR5176532.counts.txt SRR5176531.counts.txt >merged2.counts.txt

基因表达数据标准化方法

这篇好文章是转载于：学新通技术网