转录组是指细胞在某一功能状态下转录出来的所有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.
怎么样，是不是很简单呢？快来试一试吧！
扫描下方的二维码并点击关注，获取更多关于Linux、生信、实验、SCI论文发表以及最新文献资讯等方面的推送~
网友评论