干货 | 一文教会你如何采用Linux系统处理RNAseq测序数

    转录组是指细胞在某一功能状态下转录出来的所有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论文发表以及最新文献资讯等方面的推送~