转录组是指细胞在某一功能状态下转录出来的所有RNA的总和。转录组测序(Transcriptome sequencing)是基于Illumina HiSeq测序平台检测细胞内所有mRNA的一项技术,能够快速获得细胞在某一状态下所有的转录本信息,因而被广泛应用于基础研究、药物研发和临床诊断等多个领域。
在获得原始测序数据(FASTQ文件格式)后,该如何就FASTQ测序文件进行后续分析处理呢?下面,让我们一起来学习吧!
前期准备:
硬件:
Linux系统,>4G运行内存
软件:
Miniconda、FastQC、Cutadapt、Hisat2、SAMtools、subread(FeatureCounts)
PART1 测序数据质量评估采用FastQC软件进行分析
FastQC(http://www.bioinformatics.babraham.ac.uk/projects/fastqc/)是一款基于Java的软件。无论是Windows/Linux系统运行FastQC,首先都必须安装java(注意Java的版本不能超过1.6,否则FastQC会无法运行)。其作用是对测序数据进行质量评估。
FastQC安装:
$ cd path/to/fastqc
$ chmod 755 fastqc
$ ./fastqc
$ sudo ln -s /path/to/FastQC/fastqc /usr/local/bin/fastqc
运行格式:
$ fastqc -o outdir -t threads fastq1 fastq2 ...
主要参数:
-o --outdir FastQC 所生成的报告文件的储存路径
-t --threads 程序运行的线程数
示例:
$ fastqc -o FASTQC/ -t 8 Ctrl-1_combined_R1.fastq.gz Ctrl-1_combined_R2.fastq.gz Ctrl-2_combined_R1.fastq.gz Ctrl-2_combined_R2.fastq.gz Ctrl-3_combined_R1.fastq.gz Ctrl-3_combined_R2.fastq.gz Treated-1_combined_R1.fastq.gz Treated-1_combined_R2.fastq.gz Treated-2_combined_R1.fastq.gz Treated-2_combined_R2.fastq.gz Treated-3_combined_R1.fastq.gz Treated-3_combined_R2.fastq.gz
$ multiqc ./
运行结果:
运行结束后,生成两个文件:.html网页文件及.zip压缩文件。
PART2 测序数据过滤采用Cutadapt软件进行处理
测序过程中少量reads会测到接头序列,或者测序长度过长时活导致3’端碱基质量过低,因此需要对原始数据进行预处理。采用Cutadapt(https://cutadapt.readthedocs.io/en/stable/installation.html)可以帮助我们:(1)去除接头序列;(2)去除5’或3’末端质量值较低或含N的碱基;(3)去除平均质量值低于30的序列;(4)去除trim后reads长度过短的序列。
Cutadapt安装(需提前安装miniconda):
$ conda install -c bioconda cutadapt
运行格式:
$ cutadapt -a ADAPTER_FWD -A ADAPTER_REV -o out.1.fq -p out.2.fq reads.1.fq reads.2.fq
主要参数:
-a 正向接头序列(单端测序时仅有该参数)
-A 反向接头序列(双端测序时增加该参数)
-q 表示最低质量值,将低于此数值的碱基去除,一般设为30
-m 去除去接头后短于此数值的reads
--trim-n 去除含N的碱基
-o reads.1.fq去接头后的输出文件
-p reads.2.fq去接头后的输出文件
示例:
$ cutadapt -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC -A AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT -q 30 -m 75 --trim-n --report=minimal -o Ctrl-1_out1_R1.fastq.gz -p Ctrl-1_out1_R2.fastq.gz Ctrl-1_combined_R1.fastq.gz Ctrl-1_combined_R2.fastq.gz
运行结果:
获得过滤后的Clean Data。此时可再次运行一次FastQC软件查看过滤后的数据质量。
PART3 采用Hisat2软件与参考基因组进行比对分析
获得Clean data后,即可与参考基因组进行比对分析,我们采用Hisat2软件(http://ccb.jhu.edu/software/hisat2/index.shtml)进行短reads的比对,以人类参考基因组为例。
Hisat2安装:
$ unzip hisat2-2.0.1-beta-OSX_x86_64.zip
$ cp hisat2 hisat2-align-s hisat2-align-l hisat2-build hisat2-build-s hisat2-build-l hisat2-inspect hisat2-inspect-s hisat2-inspect-l $HOME/bin
$ vi ~/.bashrc
$ export PATH=/lustre/home/lcn/chenwen/bin/hisat2-2.0.1-beta:$PATH
$ source ~/.bashrc
构建索引:
i)提取剪接位点及外显子信息
$ extract_splice_sites.py genes.gtf >genome.ss
$ extract_exons.py genes.gtf >genome.exon
ii)创建HISTA2 index
$ hisat2-build -p 20 --s genome.ss --exon genome.exon genome.fa genome_tran
创建人基因组的索引往往需要很大的内存(200G RAM),因此不建议大家自己创建,常用的物种均可以去下面网站下载现成的(http://ccb.jhu.edu/software/hisat2/index.shtml)。
运行格式:
$ hisat2 -p 8 --dta -x /path/to/file/hg19/genome -1 Ctrl-1_out_R1.fastq.gz -2 Ctrl-1_out_R2.fastq.gz -S Ctrl-1_out.sam
主要参数:
-p 线程数
-x 指定索引文件
-1 指定第一个FASTQ文件
-2 指定第二个FASTQ文件
-S 输出的SAM文件
运行结果:
获得比对后的SAM文件。
PS:若想将测序文件与自己的序列进行比对,可以参考小编之前的推送。
如何将转录组数据mapping到自己的序列并可视化?(HISAT2+Samtools+IGV)
PART4 SAMtools软件将SAM文件转化为BAM文件
采用SAMtools软件(http://samtools.sourceforge.net)进行sort及格式转化,sort之后文件占位更小,其转化为BAM二进制文件。
SAMtools安装:
$ tar jxvf samtools-0.1.19.tar.bz2
$ cd samtools-0.1.19
$ make
$ cd.
$ cp samtools-0.1.19/samtools $HOME/bin
运行格式:
$ samtools sort -@ 8 -o Ctrl-1_out.bam Ctrl-1_out.sam
运行结果:
获得BAM文件,该结果可采用IGV进行可视化。
PART5 采用FeatureCounts软件进行reads计数
FeatureCounts是subread软件包中的一个工具,尽管短序列比对工具subread没有bwa和hisat2流行,但其中FeatureCounts工具却应用广泛。其作用是将比对后的结果进行计数,即计算每个区域比对到的reads数。之后对reads数进行归一化处理,计算RPKM,FPKM,TPM等,然后进行差异比较分析。
Subread安装(直接找到官网下载解压即可):
$ wget https://jaist.dl.sourceforge.net/project/subread/subread-1.6.0/subread-1.6.0-Linux-x86_64.tar.gz
$ tar -zxvf subread-1.6.0-Linux-x86_64.tar.gz
运行格式:
$ featureCounts -T 5 -p -t exon -g gene_id -a /path/to/file/genes.gtf -o Ctrl-1.featureCounts.txt /path/to/file/Ctrl-1_out.bam
其中参考基因组需自行从UCSC下载。
主要参数:
-T 多线程数
-p 针对双端测序数据
-t 指定注释文件中的功能类型,默认为外显子exon
-g 指定注释文件中的属性信息,默认为gene_id
-a 注释文件,GTF或GFF格式文件,区分外显子区域
-o 输出结果文件,同时会输出以该名称结尾的.summary统计文件
运行结果:
得到rawcounts文件。之后可在R中进行后续差异比较分析。
参考文献:
1. Mihaela Pertea, Daehwan Kim, Geo M Pertea, Jeffrey T Leek, Steven L Salzberg. Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown. Nat Prot 2016. 11(9):1650-1667.
2. Yang Liao, Gordon K. Smyth, Wei Shi. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 2014. 30(7):923-930.