生物信息笔记Bioinformatics转录组

RNA-seq数据分析---方法学文章的实战练习

2017-07-02  本文已影响9426人  面面的徐爷

前言

这次给大家带来的是16年发表在NATURE PROTOCOLS上面的一篇处理RNA-seq数据的文章:Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown
和它的12年发表于同一个杂志的兄弟文章:Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks都是NATURE PROTOCOLS上阅读量最大的文章。

NATURE PROTOCOLS

当然还有很多其它的介绍不止生信的方法文章,大家有时间可以去探究。
同时我这里就不再赘述RNA-seq的具体原理,有需要了解的请移步:一个简略的RNA-seq演示
至于软件的安装到官网下载,解压后将bin/添加进路径即可,这里不再做讲解。
#######所有操作皆在LINUX&R上完成,默认基本处理软件已经安装


本体介绍

大佬的文章

事实上作者团队一直致力于开发出更好的解决数据处理的软件,就目前12年推出的Tohat2和Cufflinks已经不是那么的令人满意了,所以他们又开发了 HISAT, StringTie and Ballgown三件套完成大家对于一个RNA-seq所有的幻想。#但实际上目前还存在着其他的性能很优秀软件也可以满足需求,例如STAR等等,但作为菜鸟来说只有先一步一步的学习了。

好了,我们先对这篇文章给出的处理方法有个大概的了解。

pipeline

  1. alignment of the reads to the genome
  2. assembly the alignments into full-length transcripts
  3. quantification of the differencesin expression levels of each gene and transcript
  4. calculation of the differences in expression for all genes among the different conditions
An overview of the 'new Tuxedo' protocol

这个protocol首先从原始RAN-seq数据入手,输出数据包括基因list,转录本,及每个样本的表达量,能够表现差异表达基因的表格并完成显著性的计算。

  1. 第一步是使用HISAT将读段匹配到参考基因组上,使用者可以提供注释文件,但HISAT依旧会检测注释文件没有列出来的剪切位点
  2. 第二步,比对上的reads将会被呈递给StringTie进行转录本组装,StringTie单独的对每个样本进行组装,在组装的过程中顺带估算每个基因及isoform的表达水平
  3. 第三步,所有的转录本都被呈递给StringTie的merge函数进行merge,这一步是必须的,因为有些样本的转录本可能仅仅被部分reads覆盖,无法被第二步的StringTie组装出来。merge步骤可以创建出所有样本里面都有的转录本,方便下一步的对比
  4. merge的数据再一次被呈递给StringTie,StringTie可以利用merge的数据重新估算转录本的丰度,还能额外的提供转录本reads数量的数据给下一步的ballgown
  5. 最后一步:Ballgown从上一步获得所有转录本及其丰度,根据实验条件进行分类统计

注意事项Warning:


数据处理设计

Read alignment with HISAT

总体上来说HISAT利用了BWA和Bowtie的算法,同时解决了mRNA中不存在内含子难以比对的问题,比对上代主流RNA-seq比对工具能快50倍,同时需求更少的内存<8G(这就意味着你可以在PC上跑数据),20个样本,每个样本一亿reads的估计,用一台电脑一天之内能够跑完。使用者可以提供精确的基因注释来提高在已知基因区域的准确性,但这是可选项。

事实上,针对RNA-seq数据的align目前有很多工具,16年12月12日出版的NATURE METHODS上的一篇文章,利用分别使用人和一种叫malaria的寄生虫的数据,模拟出三种复杂度的数据集(T1、T2和T3)。对于复杂度的衡量,定义了三个尺度:

Transcript assembly and quantification with StringTie

RNA-seq的分析依赖于精准的对于基因isoform的重建以及对于基因相对丰度的预测。继承于Cufflinks,StringTie相对于之前开发的工具更为准确,需求内存和耗时也更少。
使用者一样可以使用注释文件来帮助StringTie运行,对于低丰度的数据比较有帮助。此时StringTie依旧会对非注释区域进行转录本的组装,所以注释文件在这里是可选选项。

转录本组装完成后,使用者可以利用gffcompare(StringTie工具包含)工具来得知有多少转录本和注释文件相同,又有多少新的转录本:

#input: gff or gtf format
transcripts.gtf

command line example:
$ gffcompare -G -r chrX.gtf transcripts.gtf
#-r : reference
#-G :compare all transcripts
#-o :prefix

output file 1: gffcmp.annotated.gtf
# 显示StringTie组装的转录本与注释文件内的转录本的差别(会给每个转录本加入一个class code,我理解为一个标识,释义如下图)
output dile 2: gffcmp.stats
# 根据注释文件显示组装结果的准确性和阳性预测率

class codes

由于在实验中,我们会处理多个RNA-seq数据,单个样本里面的基因和iosforms与其它样本的数据很少相同,但是为了进行比较就需要它们以一个相同的形式进行组装,作者通过StringTie的merge工具解决了这个问题,将所有样本中发现的基因进行merge。下图的例子可以帮助理解。


example

如图所示,四个样本中组装得到四个转录本,最后merge得到A/B两个转录本。1/2样本与参考基因组相同merge得到A转录本,3/4样本相互吻合但与参考基因组不一样,所以merge得到B转录本。
在merge步骤之后,需要StringTie再运行一遍去重新估算merge得到的转录本的丰度。

Chevy理解:这就是相当于重新定义了一个annotation file进行二次重新组装,相对于第一次组织可以提高准确度和敏感性。比已有注释文件的优势在于可以发现更多的isoforms。

Differential expression analysis with Ballgown

组装和定量完转录本之后自然需要进行基因表达差异分析,统计建模和可视化。R和Bioconductor提供了一揽子的工具处理这些任务,包括raw data的plot作图,数据的标准化,下游的统计建。此处使用Ballgown包作为承上启下的桥梁。

在R里面使用Ballgown处理需要得到:

而大多数的差异表达分析都会包括下面几个步骤:

Ballgown包可以完成以上的几个步骤并且可以联合R语言的其它操作进行分析

Ballgown使用:

实战示例

准备阶段

实际上本轮操作就是按照文章给的示例数据走了一遍流程:

文章中,为了减少计算时间,方便读者重复,作者截取了一批实验数据中能够匹配到chromosomeX上的数据作为示例数据,chromosomeX是一个基因相对较为丰富的染色体,可以占到基因组的5%左右。

首先是下载数据随后解压:

$nohup wget ftp://ftp.ccb.jhu.edu/pub/RNAseq_protocol/chrX_data.tar.gz 2>download.log &
$tar zxvf chrX_data.tar.gz

其中包括如下数据:


文件数据

sample文件夹下有12个fastq格式的paired-end RNA-seq文件,从两地人群中(英格兰岛住民和约鲁巴住民)各取六个样,六个样又分为男女性别各三个(最少的生物学重复)。

sample

随后介绍一个骚命令:

HISAT2可以直接从NCBI下载sra格式的文件并作为输入文件格式
下面以ERR188245测序数据为例
$ hisat2 -p 8 --dta -x chrX_data/indexes/chrX_tran --sra-acc ERR188245 –S ERR188245_chrX.sam
可以说是相当的帮人偷懒了

index文件夹下包含着HISAT2对于染色体X的index文件

index
当然HISAT2已经为各位老爷准备好了常见基因组的index文件和genes文件。
点我

genome文件夹下包含这染色体X的序列文件(GRCh38 build 81)

genome

genes文件夹下则包含着针对基因组的注释文件(来自于RefSeq数据库)


annotation

geuvadis_phenodata.csv和mergelist.txt则是用户着要自己创建表明数据关系的文件,这里作者提供出来方便使用。

如果需要建立官网上没有基因组的index,就需要需要使用Hisat2包里面的python脚本:
extract_splice_sites.py和extract_exons.py,从注释文件里面抽取出剪切位点和外显子信息

$ extract_splice_sites.py chrX_data/genes/chrX.gtf >chrX.ss
$ extract_exons.py chrX_data/genes/chrX.gtf >chrX.exon
先提取剪切位点和外显子数据
$ hisat2-build --ss chrX.ss --exon chrX.exon chrX_data/genome/chrX.fa chrX_tran
随后建立HISAT2 index

正式开始

1.将reads比对上genome

2017-07-24更新:我自己使用的是Ensembl发表的GRCh38版本的基因组进行比对,所以这里我附上这个版本的genome和注释文件的下载地址:

Ensembl版本全基因组的注释文件下载
HISAT2-index下载genome_tran

# 样本比对操作
nohup hisat2 -p 8 --dta -x ~/RNA-seq/chrx/chrX_data/indexes/chrX_tran -1 ~/RNA-seq/chrx/chrX_data/samples/ERRERR188044_chrX_1.fastq.gz -2 ERR188044_chrX_2.fastq.gz -S ERR188044_chrX.sam &

# 数据太多我就写了个小脚本处理,in sample dir/
for i in *1.fastq.gz; 
do
i=${i%1.fastq.gz*}; 
nohup hisat2 -p 8 --dta -x ~/RNA-seq/chrx/chrX_data/indexes/chrX_tran -1 ${i}1.fastq.gz -2 ${i}2.fastq.gz -S ~/RNA-seq/chrx/chrX_data/result/${i}align.sam 2>~/RNA-seq/chrx/chrX_data/result/${i}align.log & 
done

运行结果如下:


比对结果

2.将sam文件转化为bam文件

# 转化操作
samtools sort @ 8 -o ERR188044_chrX.bam ERR188044_chrX.sam

# 批处理, in result dir/
for i in *.sam;
do 
i=${i%_align.sam*}; nohup samtools sort -@ 8 -o ${i}.bam ${i}_align.sam & 
done

结果如下:

sam to bam

3.组装转录本并定量表达基因

# 单文件操作
stringtie -p 8 -G ~/RNA-seq/chrx/chrX_data/genes/chrX.gtf -o ERR188044_chrX.gtf ERR188044_chrX.bam

# 批处理, in result dir/
for i in *.bam; 
do 
i=${i%.bam*}; nohup stringtie -p 8 -G ~/RNA-seq/chrx/chrX_data/genes/chrX.gtf -o ${i}.gtf ${i}.bam & 
done

结果如下:


strigtie

4.将所有转录本合并

warning: 此处的mergelist.txt是自己创建的,需要包含之前output.gtf文件的路径

cat ~/RNA-seq/chrx/chrX_data/mergelist.txt
ERR188044_chrX.gtf
ERR188104_chrX.gtf
ERR188234_chrX.gtf
ERR188245_chrX.gtf
ERR188257_chrX.gtf
ERR188273_chrX.gtf
ERR188337_chrX.gtf
ERR188383_chrX.gtf
ERR188401_chrX.gtf
ERR188428_chrX.gtf
ERR188454_chrX.gtf
ERR204916_chrX.gtf
#因为就在结果目录下面操作,所以直接显示文件名即可

stringtie --merge -p 8 -G  ~/RNA-seq/chrx/chrX_data/genes/chrX.gtf -o stringtie_merged.gtf  ~/RNA-seq/chrx/chrX_data/mergelist.txt

#output文件即为
stringtie_merged.gtf 

5.检测相对于注释基因组,转录本的组装情况

gffcompare -r ~/RNA-seq/chrx/chrX_data/genes/chrX.gtf  -G -o merged stringtie_merged.gtf
#输出文件信息可见上面的Transcript assembly and quantification with StringTie介绍

6.重新组装转录本并估算基因表达丰度,并为ballgown创建读入文件

# 单文件操作
stringtie -e -B -p 8 -G stringtie_merged.gtf -o ballgown/ERR188044_chrX/ERR188044_chrX.gtf ERR188044_chrX.bam

# 批处理, in result dir/
for i in *_chrX.bam; 
do 
i=${i%_chrX.bam*}; nohup stringtie -e -B -p 8 -G stringtie_merged.gtf -o ballgown/${i}/${i}_chrX.gtf  ${i}_chrX.bam & 
done

结果如下:

table

7.R包的安装与载入

R语言需要安装ballgown包和一些接下来处理会用到的包(RSkittleBrewer/genefilter/dplyr/devtools)
事实上我还没搞清楚LINUX里面的R包安装问题,老是报error,所以这里我就用Windows下的R进行操作。

# 首先是几个包的载入
>library(ballgown)
>library(RSkittleBrewer)
>library(genefilter)
>library(dplyr)
>library(devtools)

8.数据的读入

# 随后是数据的读入,CSV文件,我把所有文件放在桌面,上一步得到的ballgown文件夹直接拷到桌面
> setwd("C:/Users/DELL/Desktop")
> read.csv("geuvadis_phenodata.csv")
         ids    sex population
1  ERR188044   male        YRI
2  ERR188104   male        YRI
3  ERR188234 female        YRI
4  ERR188245 female        GBR
5  ERR188257   male        GBR
6  ERR188273 female        YRI
7  ERR188337 female        GBR
8  ERR188383   male        GBR
9  ERR188401   male        GBR
10 ERR188428 female        GBR
11 ERR188454   male        YRI
12 ERR204916 female        YRI
> pheno_data<-read.csv("geuvadis_phenodata.csv")
> bg_chrX=ballgown(dataDir = "C:/Users/DELL/Desktop/ballgown",samplePattern = "ERR",pData = pheno_data)
# dataDir指的是数据储存的地方,samplePattern则依据样本的名字来,pheno_data则指明了样本数据的关系,这个里面第一列样本名需要和ballgown下面的文件夹的样本名一样,不然会报错

10.滤掉低丰度的基因

# 这里滤掉了样本间差异少于一个转录本的数据
> bg_chrX_filt=subset(bg_chrX,"rowVars(texpr(bg_chrX))>1",genomesubset=TRUE)

11.确认组间有差异的转录本

2017-07-21更新:实际上为了确定基因表达差异是由于某个变量造成的,我们这边需要使用剩下的变量进行修正。Example:我们为了比较male and female之间的基因表达差异,这里就需要指定分析参数为"transcripts",主变量为"sex",修正变量为"population",getFC可以指定输出结果显示组间表达量的foldchange。

并且Ballgown的统计算法基于标准的线性模型,对于组内数据较少(<4)时较为准确。这里也可以使用其它的一些软件例如limma/DESeq/edgeR等。

> result_transcripts=stattest(bg_chrX_filt,feature = "transcript",covariate = "sex",adjustvars = c("population"), getFC=TRUE,meas="FPKM")
> result_transcripts
        feature   id           fc         pval         qval
1    transcript    1  0.941367186 7.353184e-01 9.468599e-01
2    transcript    2  1.207949354 8.668638e-01 9.744055e-01
3    transcript    3  1.013100019 9.920824e-01 9.986659e-01
4    transcript    4  0.387671042 5.233364e-01 9.189906e-01
5    transcript    5  0.615797877 3.324164e-01 9.117015e-01
......

12.确认组间有差异的基因

这里指定feature为gene

> result_genes=stattest(bg_chrX_filt,feature = "gene",covariate = "sex",adjustvars = c("population"), getFC=TRUE,meas="FPKM")
> result_genes
     feature           id          fc         pval         qval
1       gene      MSTRG.1 1.114999374 7.305975e-01 9.218157e-01
2       gene     MSTRG.10 1.749747142 2.767538e-01 7.922755e-01
3       gene   MSTRG.1000 1.399358968 6.039065e-01 8.945639e-01
4       gene   MSTRG.1002 0.937668743 8.825216e-01 9.816878e-01
5       gene   MSTRG.1003 1.044840944 6.155345e-01 8.977043e-01
6       gene   MSTRG.1004 1.220626116 4.160214e-01 8.388688e-01
.....

13.对结果 result_transcripts增加基因名

> result_transcripts=data.frame(geneNames=ballgown::geneNames(bg_chrX_filt),geneIDs=ballgown::geneI Ds(bg_chrX_filt),result_transcripts)
> result_transcripts
          geneNames      geneIDs    feature   id           fc         pval         qval
1                 .      MSTRG.4 transcript    1  0.941367186 7.353184e-01 9.468599e-01
2            PLCXD1      MSTRG.4 transcript    2  1.207949354 8.668638e-01 9.744055e-01
3                 .      MSTRG.4 transcript    3  1.013100019 9.920824e-01 9.986659e-01
4                 .      MSTRG.4 transcript    4  0.387671042 5.233364e-01 9.189906e-01
5                 .      MSTRG.5 transcript    5  0.615797877 3.324164e-01 9.117015e-01
6            PLCXD1      MSTRG.4 transcript    6  0.630018786 2.938396e-01 9.002815e-01
......

14.按照P值排序(从小到大)

> result_transcripts=arrange(result_transcripts,pval)
> result_genes=arrange(result_genes,pval)

15.把结果写入csv文件

> write.csv(result_transcripts, "chrX_transcript_results.csv",row.names=FALSE)
> write.csv(result_genes, "chrX_gene_results.csv",row.names=FALSE)

16.筛选出q值小于0.05的transcripts和genes,也就是在性别间表达有差异的基因了。

注:这里无需过分关注geneIDs以及transcript name,实际上那个是stringtie在比对过程中赋予的一个符号。

> subset(result_transcripts,result_transcripts$qval<0.05)
   geneNames   geneIDs    feature   id          fc         pval         qval
1          . MSTRG.531 transcript 1657 0.030820591 1.427864e-10 2.082377e-07
2       XIST MSTRG.531 transcript 1656 0.003014576 1.860927e-10 2.082377e-07
3          . MSTRG.531 transcript 1655 0.016144997 3.762791e-10 2.807042e-07
4          . MSTRG.531 transcript 1658 0.028308029 6.599039e-08 3.692163e-05
5       TSIX MSTRG.530 transcript 1654 0.078461818 1.690716e-06 7.567643e-04
6          . MSTRG.613 transcript 1848 7.391342987 1.249275e-05 4.659796e-03
7          . MSTRG.141 transcript  421 3.204058932 2.729898e-05 8.727874e-03
8          . MSTRG.618 transcript 1852 9.136857610 4.804244e-05 1.343987e-02
9          . MSTRG.777 transcript 2338 0.127674288 1.000751e-04 2.488533e-02
10     KDM6A MSTRG.258 transcript  736 0.054212867 1.173824e-04 2.627018e-02
11    PNPLA4  MSTRG.62 transcript  186 0.592778584 2.083496e-04 4.238966e-02
12     RPS4X MSTRG.519 transcript 1611 0.597532121 2.493976e-04 4.651265e-02

> subset(result_genes,result_genes$qval<0.05)
   feature        id          fc         pval         qval
1     gene MSTRG.531 0.002396124 2.469125e-11 2.523446e-08
2     gene MSTRG.141 3.165966412 1.666206e-06 8.514310e-04
3     gene MSTRG.530 0.082749314 7.018439e-06 2.390948e-03
4     gene MSTRG.613 7.314423295 1.128589e-05 2.883544e-03
5     gene MSTRG.618 9.050399157 5.022017e-05 1.026500e-02
6     gene MSTRG.519 0.637382385 6.953432e-05 1.184401e-02
7     gene MSTRG.376 0.620693674 1.134561e-04 1.656458e-02
8     gene  MSTRG.62 0.600686117 1.638615e-04 2.093330e-02
9     gene MSTRG.777 0.159178288 2.177206e-04 2.472338e-02
10    gene MSTRG.622 7.870447732 3.737270e-04 3.527295e-02
11    gene MSTRG.778 0.127712244 3.796501e-04 3.527295e-02
12    gene MSTRG.229 1.412379185 4.200674e-04 3.577574e-02
13    gene  MSTRG.48 4.158987103 5.279898e-04 4.150812e-02

17.数据可视化之颜色设定

# 赋予调色板五个指定颜色
> tropical= c('darkorange', 'dodgerblue','hotpink', 'limegreen', 'yellow')
> palette(tropical)

# 当然rainbow()函数也可以完成这个任务
> palette(rainbow(5))

18.以FPKM为参考值作图,以性别作为区分条件

# 提取FPKM值
> fpkm = texpr(bg_chrX,meas="FPKM")

#方便作图将其log转换,+1是为了避免出现log2(0)的情况
> fpkm = log2(fpkm+1)

# 作图
> boxplot(fpkm,col=as.numeric(pheno_data$sex),las=2,ylab='log2(FPKM+1)')

结果如下:

FPKM, male in blue, females in orange

19.就单个转录本的查看在样品中的分布

> ballgown::transcriptNames(bg_chrX)[12]
         12 
"NR_027232" 
> ballgown::geneNames(bg_chrX)[12]
         12 
"LINC00685" 

# 绘制箱线图
> plot(fpkm[12,] ~ pheno_data$sex, border=c(1,2),
+      main=paste(ballgown::geneNames(bg_chrX)[12],' : ',
+                 ballgown::transcriptNames(bg_chrX)[12]),pch=19, xlab="Sex",
+      ylab='log2(FPKM+1)')
> points(fpkm[12,] ~ jitter(as.numeric(pheno_data$sex)),
+        col=as.numeric(pheno_data$sex))

结果如下:

transcript

20.查看某一基因位置上所有的转录本

# plotTranscripts函数可以根据指定基因的id画出在特定区段的转录本
# 可以通过sample函数指定看在某个样本中的表达情况,这里选用id=1750, sample=ERR188234
> plotTranscripts(ballgown::geneIDs(bg_chrX)[1750], bg_chrX, main=c('Gene MSTRG.575 in sample ERR188234'), sample=c('ERR188234'))

结果如下:


21.以性别为区分查看基因表达情况

# 这里以id=575的基因为例(对应上一步作图)
> plotMeans('MSTRG.575', bg_chrX_filt,groupvar="sex",legend=FALSE)

结果如下:

MSTRG.575

临终讨论


参考文献

Pertea M, Kim D, Pertea G M, et al. Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown[J]. Nature protocols, 2016, 11(9): 1650-1667.


日常Bob镇楼
上一篇下一篇

猜你喜欢

热点阅读