RNA-seq技术是指通过现有的测序方法技术手段获取某个物种或者特定细胞类型产生的所有转录本的集合。 转录组研究能够从整体水平研究基因功能以及基因结构,揭示特定生物学过程以及疾病发生过程中的分子机理,已被广泛应用于基础研究、临床诊断和药物研发等领域。
现在基本上都是利用的illumina平台进行转录组测序,illumina的测序文件中,一般采用双端测序(paired-end),一个样本得到的是seq_1.fastq.gz和seq_2.fastq.gz两个文件,每个文件存放一段测序文件。在illumina公司测得的序列文件经过处理以fastq文件协议存储为*.fastq、*.fq、*.fq.gz等fastq格式文件,分析转录组数据主要包含以下几个步骤:
- 原始数据质控
- 去除adapter和低质量reads
- 和参考序列比对
- 转录本组装、获取每个基因/转录本count数
- 基于基因的表达count数进行下游数据分析
- 基因富集、相互作用分析(网络关系)
注:下面提到的软件推荐利用conda安装,代码现在各个搜索引擎都有很多,本文不在添加具体代码,感兴趣的小伙伴可以打开官方文档或者借助搜索引擎进行探索。
1)原始数据质控
可以利用像fastqc(http://www.bioinformatics.babraham.ac.uk/projects/fastqc/)、multiqc(https://multiqc.info/)等来做原始数据的质控,软件分别能给出最终质控情况和完整的报告。
2)去除adapter和低质量reads
在去除adapter和低质量可以选择使用:fastqc、trimmomatic(http://www.usadellab.org/cms/?page=trimmomatic)、cutadapt(https://cutadapt.readthedocs.io/en/stable/)、trim-galore(https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/)等软件(注意参数的使用,不match的参数可能会导致异常的结果)
3)和参考序列比对
这个可以选用的软件有很多,如:star(推荐,https://github.com/alexdobin/STAR)、hisat2(推荐,http://daehwankimlab.github.io/hisat2/) 、 tophat2(https://ccb.jhu.edu/software/tophat/index.shtml)、bowtie2 (http://bowtie-bio.sourceforge.net/bowtie2/index.shtml)、 bwa(http://bio-bwa.sourceforge.net/)等,当然这根据你所需要比对的样本或者序列相关,某些软件擅长长序列或者针对某个平台下数据的灵敏度较高,某些软件可能针对短序列比对处理较好,或者某些软件处理速度较快等。需要根据具体的情况选择,此处不在赘述。
4)转录本组装、获取每个基因/转录本count数
序列比对完成后可以利用诸如samtools(http://samtools.sourceforge.net/)等软件将sam文件转换为bam文件(bam,二进制文件,方便排序)。然后利用htseq-count(推荐,现在支持多线程,https://htseq.readthedocs.io/en/master/index.html)、stringtie(推荐,https://ccb.jhu.edu/software/stringtie/)、bedtools(https://bedtools.readthedocs.io/en/latest/)、salmon(http://combine-lab.github.io/salmon/)、Subreads(http://subread.sourceforge.net/)将统计得到的每个基因/转录本的count数。
5)基于基因的表达count数进行下游数据分析
对于得到的count数,我们通常使用的方法有基于序列的CPM(counts per million)、log-CPM、FPKM(fragments per kilobase of transcript per million),和基于转录本数目的RPKM(reads per kilobase of transcript per million)等方法。此处我们可以参考RNAseq123
教程内容,用CPM进行下游分析,假设基因长度在比较中始终不变,且任何观测到的差异是来自于条件的变化而不是基因长度的变化。
对于RNA-seq数据,来自Bioconductor(https://bioconductor.org/)的 edgeR、DESeq2和limma包提供了一套用于处理此问题的完善的统计学方法(具体代码可以查看参考资料1,推荐用R分析)。
- count数据标准化(如CPM)
- 删除低表达基因
- 归一化数据,无监督聚类查看样本情况(如PCA、t-SNE、UMAP等)
- 差异表达分析(可以用像edgeR、DESeq2和limma等包)
6)基因富集、相互作用分析(网络关系)
对得到的差异基因进行下游的数据分析,包括利用GO/KEGG/Drug等数据获取差异基因所关联的已知的生物意义、功能,对经典研究进行深入探索和佐证或者推翻。
可以利用诸如GSEA、String等对通路进行深入探究。
7)高级分析
现在比较流行的利用ML(或者DL)来对大数据进行深度挖掘,寻找内在的关联模式,个人经验现在高影响力的文章主要用到的方法诸如:逻辑回归、RF、SVM、CNN等方法,如果是深入的数据挖掘可能会用到多阶机器学习方法。一般来说,对于常规的数据分析到第6步基本上算是可以close掉,准备下游的验证和投稿事宜了,如果数据涉及到队列可以利用第七步建议进行深入挖掘或许能发现一些更有意思的规律。
参考资料:
1.http://master.bioconductor.org/packages/release/workflows/vignettes/RNAseq123/inst/doc/limmaWorkflow_CHN.html
2.https://baike.baidu.com/item/RNA-seq