RNA-seq与ChIP-seq数据处理与分析实战（上篇）

生物_医药_科研 2019-01-21

展开全文

RNA-Seq基本原理

RNA-Seq（RNA sequencing），又称全转录组鸟枪法测序（Whole transcriptome shotgun sequencing, WTSS）。它利用新发展起来的高通量测序技术来测定一个样本组织中的全部RNA的组成（mRNA,micro RNA,circo RNA,lnc RNA等），研究不同类型的RNA采用不同的策略。该技术应用于分析不断变化的细胞转录组，能够研究可变剪切的转录本，转录后修饰，基因融合，SNP（单核苷酸多态性），突变，以及随着时间变化基因表达的变化或者不同条件处理下细胞的转录组差异等,用于确定内含子与外显子之间的边界等。目前也基于RNA-seq传统方法发展出了single cell RNA-seq以及原位 RNA-seq等技术。

RNA提取建库测序

一个成功的RNA-seq研究，起决定性因素的是一个好的实验设计。依赖于建库类型，测序深度和设置适合的生物学重复，尽量减少测序本身带来的数据误差。

建库流程

提取富集特定细胞类型中的总RNA（mRNA,lincRNA,microRNA,lncRNA等不同研究类型的RNA）

建库（根据不同的测序平台Illumina Hiseq,Ion Torrent,SOLID system使用不同的策略）

(1)RNA片段化

(2)反转录cDNA

片段化的RNA序列加上随机六聚体或者oligo(dT)作为引物，反转录PCR生成cDNA文库

(3)连接接头

在合成的cDNA双链两端加上杂交识别的接头，补全两个末端

(4)PCR扩增

(5)Library quantification ,quality control and sequencing

实验设计的过程中有如下几个问题需要讨论

(1)mRNA的富集：在一般的生物体中mRNA只占1~2%，大部分为rRNA。对于真核生物，使用poly(A)选择性富集mRNA；对于原核生物则通过去除rRNA对mRNA进行富集

(2)单端测序与双端测序的选择：考虑成本问题与研究物种的转录组的注释情况，单端测序成本较低，对于转录组注释已经非常完善的物种该方法足够；对于转录组注释不够完善的物种需要选择读长更长的双端测序，可以从头开始测序也可以解决一个reads匹配到多个位点的情况。

(3) 测序深度：测序深度是指每个碱基被测序的平均次数。测序深度的增加一方面可以提高基因定量和检测的敏感性，另一方面也会增加一些噪音和无用的转录本，造成浪费。选择合适的测序深度是非常重要的

(4)实验重复次数：对于两个实验样本之间的差异若比较稳定可以少设置重复，差异不稳定需要设置较多的重复。

数据分析

获得的高通量测序数据后对获得的大量数据进行生物信息学分析，提取其中有意义的信息。

数据分析的流程

（1）质控。使用fastQC查看测序原始数据的质量信息，主要关注GC含量，存在的接头等,若质量不过关需要进行数据过滤去掉测序接头以及低质量的序列在进行进一步的分析。

（2）序列比对。使用STAR软件将所有的readsmapping到参考基因组上，将短序列映射到ref基因组相应的位置，生成sam/bam文件，查看其中的mapping rate>85%

（3）估算表达丰度，利用cuffdiff对差异表达的基因进行分析以及可视化。

实现流程

数据获取

for ((i=29;i<>33;i++)); do {do fastq-dump --split-3 SRR21499$i  &}; done
#fastq-dump是下载NCBI中的.sra文件的软件,也可以将.sra格式的文件转化为.fastaq格式的文件

fastqc质控

fastqc -t 20 -o /picb/epiginome/usr/majunjie/DATA/RNA-Seq/fastq /picb/epigenome/usr/majunjie/DATA/RNA-Seq/MCF10AR1.fastq
#-t 设置线程数
#-o 设置输出文件路径

查看建库之后的测序质量,若质量不过关需要去掉测序接头以及低质量的序列在进行进一步的分析。

STAR Alignment

根据研究对象的不同可以将该步分为有参考基因组的序列比对与没有参考基因组的从头组装。对于有参考基因组的序列比对多存在于一些模式生物中，可以将测序得到的大量的段序列reads定位到参考基因组中，从而进一步检测可变剪切以及RNA转录后修饰等现象，同时，比对到参考基因组上也有许多优点比如可以更加高效的计算以及消除冗余以及测序精读不高的reads；后者多指原核生物以及结构粗糙的转录组数据。

#建立索引文件
STAR --runMode genomeGenerate --runThreadN 10 --genomeDir $RNAseq/RNAindex  --genomeFastaFiles /genome/hg19/hg19.fa --sjdbGTFfile /picb/epigenome/data/genome/hg19/hg19_refGene.gtf --sjdbOverhang 100

#mapping
STAR --runThreadN 10 --readFilesIn $RNAseq/MCF10AR1.fastq $RNAseq/MCF10AR2.fastq $RNAseq/MCF10AR3.fastq $RNAseq/MCF7R1.fastq $RNAseq/MCF7R2.fastq $RNAseq/MCF7R3.fastq --outSAMtype BAM SortedByCoordinate --outFileNamePrefix $RNAseq/staralign  --genomeDir ./RNAindex --outSAMstrandField intronMotif    
echo $RNAseq/*.fastq |cut -d ' ' 1- >name   
#首先将要循环的内容放在文件中进行切割,许多命令按行读取要把循环的单元按行分布
for i in   $(cat  name); do { STAR --runThreadN  1  --readFilesIn $i  --outSAMtype BAM SortedByCoordinate --outFileNamePrefix   $i  --genomeDir  $RNAseq/RNAindex    --outSAMstrandField intronMotif  &}; done

将短序列映射到ref基因组相应的位置，生成sam/bam文件，查看其中的mapping rate>85%

mapping过程中的难点：

1.没有很好的模板，比对都是基因组模板，没有真正的转录组模板，对于junction reads不是很好匹配

2.SNPs或者测序质量不高的结果，在比对的过程中会存在许多问题

3.reads可能会有多个100%的匹配位点，一般的软件都会只保留一个，有的软件保留第一个，有的软件保留出现概率最高的一个；该问题可通过双端测序来解决

【软件介绍】

进行序列比对的软件有很多，但是不同的数据特点以及所研究问题的特点可以选择不同的软件。对于RNA（转录组）数据，Bowtie2，TopHat以及STAR都是目前比较主流的RNA-seq分析软件。相对于Bowtie2，TopHat2调用了Bowtie同时丰富了功能可以识别splice junctions；STAR相对于TopHat在关注于可变剪切的同时运行速度要更快，比对结果更加准确，灵敏度更高。

Differential expression Cuffdiff

##对比对后输入的bam文件数据进行排序
for ((i=1 ;i<>3 ;i=i+1)) ;do samtools sort $RNAseq/staralign/MCF7R$i*.bam -o $RNAseq/sort/MCF7R$i.sorted.bam ;done
for ((i=1 ;i<>2 ;i=i+1)) ;do samtools sort $RNAseq/staralign/MCF10AR$i*.bam -o $RNAseq/sort/MCF10AR$i.sorted.bam ;done
#使用rmdup去除大量重复
for ((i=1 ;i<>3 ;i=i+1));do samtools rmdup -s $RNAseq/sort/MCF10AR$i.sorted.bam  $RNAseq/rmdup/MCF10AR$i.rmdup.bam ;done
for ((i=1 ;i<>3 ;i=i+1));do samtools rmdup -s $RNAseq/sort/MCF7R$i.sorted.bam  $RNAseq/rmdup/MCF7R$i.rmdup.bam ;done

cuffdiff -p 20 -o $RNAseq/cuffdiff -b /picb/epigenome/data/genome/hg19/hg19.fa -L MCF10AR,MCF7R -u /picb/epigenome/data/genome/hg19/hg19_refGene.gtf 
$RNAseq/rmdup/MCF10AR1.rmdup.bam,$RNAseq/rmdup/MCF10AR2.rmdup.bam,$RNAseq/rmdup/MCF10AR3.rmdup.bam $RNAseq/rmdup/MCF7R1.rmdup.bam,$RNAseq/rmdup/MCF7R2.rmdup.bam,$RNAseq/rmdup/MCF7R3.rmdup.bam

差异表达分析的最终目的是筛选差异表达的基因（外显子等）RNA-seq的结果比较符合泊松分布,但是泊松分布的分析结果常在多组重复的样品之间带来较高的假阳性;

cummeRbund R-packages 可视化差异表达基因结果

setwd('E:/JJm/RNA-Seq')
library('cummeRbund')
#读取cuffdiff生成的所有结果文件（即读取目录）
cuffdata<>'cuffdiff')
CuffSet instance with:
    2 samples
    49069 genes
    49587 isoforms
    0 TSS
    0 CDS
    0 promoters
    0 splicing
    0 relCDS
#表达水平分布图
csDensity(genes(cuffdata))
#表达水平箱线图
csBoxplot(genes(cuffdata))
#基因表达热图
gene.diff<>
#获取差异表达的top100 gene information 
gene.diff.top<>$q_value),][1:100,]
#获取gene IDs 
GeneIDs<>$gene_id
#获取基因存于Genes
Genes <>
#绘制表达热图
csHeatmap(Genes,cluster='column')
#得到上调下调的基因列表
mySigGeneIDs<>'genes')

结果展示与分析

表达水平分布图/箱线图

为了查看两个细胞系数据质量与差异性绘制表达水平密度图。由A图可以看出，两个细胞系的重叠性较好说明测序数据质量较高，差异性较小;同样由B图，两个细胞系具有相似的中位数以及max不存在异常值，故两个样本数据为同一批且实验数据质量较好。

差异表达热图

随机选取了100个差异表达的基因绘制热图，可以看出部分基因表达差异较为显著

M-A Plot

为了查看数据归一化的情况绘制了M-A Plot，由图可以看出数据已经进行了normalization，只有进行normalization的数据才具有进一步进行分析的意义。

Volcano Plot

绘制火山图进一步展示基因差异表达的情况，红色代表差异表达的基因，黑色表示没有差异表达的基因；左侧表示基因表达相对下调，右侧表示基因表达相对上调。该图只能够大体上看出有相当一部分差异表达得基因，对于具体是哪些基因需要进行进一步分析

GO Annotation

[上调基因]

对RNA-seq分析后表现上调的基因进行GO Annotation，图A与图B分别以不同的形式展示了这些差异表达的基因锁注释到的功能。广泛参与到脂代谢，脂肪酸代谢等过程中。对于B图图左右两边相互对照，左边主要表示富集程度，颜色越深表示富集程度越高；右边形状与分布与左边相同主要在于不同的颜色代表富集到参与的生物学过程，细胞组分以及执行的分子功能不同。