????轉(zhuǎn)錄組是指細(xì)胞在某一功能狀態(tài)下轉(zhuǎn)錄出來的所有RNA的總和恬汁。轉(zhuǎn)錄組測(cè)序（Transcriptome sequencing）是基于Illumina HiSeq測(cè)序平臺(tái)檢測(cè)細(xì)胞內(nèi)所有mRNA的一項(xiàng)技術(shù)伶椿，能夠快速獲得細(xì)胞在某一狀態(tài)下所有的轉(zhuǎn)錄本信息，因而被廣泛應(yīng)用于基礎(chǔ)研究蕊连、藥物研發(fā)和臨床診斷等多個(gè)領(lǐng)域悬垃。

? ? 在獲得原始測(cè)序數(shù)據(jù)（FASTQ文件格式）后游昼，該如何就FASTQ測(cè)序文件進(jìn)行后續(xù)分析處理呢甘苍？下面，讓我們一起來學(xué)習(xí)吧烘豌！

前期準(zhǔn)備：

硬件：

Linux系統(tǒng)载庭，>4G運(yùn)行內(nèi)存

軟件：

Miniconda、FastQC廊佩、Cutadapt囚聚、Hisat2、SAMtools标锄、subread（FeatureCounts）

PART1 測(cè)序數(shù)據(jù)質(zhì)量評(píng)估采用FastQC軟件進(jìn)行分析

? ? FastQC（http://www.bioinformatics.babraham.ac.uk/projects/fastqc/）是一款基于Java的軟件顽铸。無論是Windows/Linux系統(tǒng)運(yùn)行FastQC，首先都必須安裝java（注意Java的版本不能超過1.6料皇，否則FastQC會(huì)無法運(yùn)行）谓松。其作用是對(duì)測(cè)序數(shù)據(jù)進(jìn)行質(zhì)量評(píng)估。

FastQC安裝：

$ cd path/to/fastqc

$ chmod 755 fastqc

$ ./fastqc

$ sudo ln -s /path/to/FastQC/fastqc /usr/local/bin/fastqc

運(yùn)行格式：

$ fastqc -o outdir -t threads fastq1 fastq2 ...

主要參數(shù)：

-o --outdir FastQC 所生成的報(bào)告文件的儲(chǔ)存路徑

-t --threads 程序運(yùn)行的線程數(shù)

示例：

$ 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 ./

運(yùn)行結(jié)果：

運(yùn)行結(jié)束后践剂，生成兩個(gè)文件：.html網(wǎng)頁文件及.zip壓縮文件鬼譬。

PART2 測(cè)序數(shù)據(jù)過濾采用Cutadapt軟件進(jìn)行處理

? ? 測(cè)序過程中少量reads會(huì)測(cè)到接頭序列，或者測(cè)序長度過長時(shí)活導(dǎo)致3’端堿基質(zhì)量過低逊脯，因此需要對(duì)原始數(shù)據(jù)進(jìn)行預(yù)處理优质。采用Cutadapt（https://cutadapt.readthedocs.io/en/stable/installation.html）可以幫助我們：(1)去除接頭序列；(2)去除5’或3’末端質(zhì)量值較低或含N的堿基军洼；(3)去除平均質(zhì)量值低于30的序列巩螃；(4)去除trim后reads長度過短的序列。

Cutadapt安裝（需提前安裝miniconda）：

$ conda install -c bioconda cutadapt

運(yùn)行格式：

$ cutadapt -a ADAPTER_FWD -A ADAPTER_REV -o out.1.fq -p out.2.fq reads.1.fq reads.2.fq

主要參數(shù)：

-a 正向接頭序列（單端測(cè)序時(shí)僅有該參數(shù)）

-A 反向接頭序列（雙端測(cè)序時(shí)增加該參數(shù)）

-q 表示最低質(zhì)量值匕争，將低于此數(shù)值的堿基去除避乏，一般設(shè)為30

-m 去除去接頭后短于此數(shù)值的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

運(yùn)行結(jié)果：

獲得過濾后的Clean Data。此時(shí)可再次運(yùn)行一次FastQC軟件查看過濾后的數(shù)據(jù)質(zhì)量汗捡。

PART3 采用Hisat2軟件與參考基因組進(jìn)行比對(duì)分析

? ? 獲得Clean data后淑际，即可與參考基因組進(jìn)行比對(duì)分析畏纲，我們采用Hisat2軟件（http://ccb.jhu.edu/software/hisat2/index.shtml）進(jìn)行短reads的比對(duì)，以人類參考基因組為例春缕。

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

構(gòu)建索引：

i)提取剪接位點(diǎn)及外顯子信息

$ extract_splice_sites.py genes.gtf >genome.ss

$ extract_exons.py genes.gtf >genome.exon

ii)創(chuàng)建HISTA2 index

$ hisat2-build -p 20 --s genome.ss --exon genome.exon genome.fa genome_tran

? ? 創(chuàng)建人基因組的索引往往需要很大的內(nèi)存（200G RAM）盗胀，因此不建議大家自己創(chuàng)建，常用的物種均可以去下面網(wǎng)站下載現(xiàn)成的（http://ccb.jhu.edu/software/hisat2/index.shtml）锄贼。

運(yùn)行格式：

$ 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

主要參數(shù)：

-p 線程數(shù)

-x 指定索引文件

-1 指定第一個(gè)FASTQ文件

-2 指定第二個(gè)FASTQ文件

-S 輸出的SAM文件

運(yùn)行結(jié)果：

獲得比對(duì)后的SAM文件票灰。

PS：若想將測(cè)序文件與自己的序列進(jìn)行比對(duì)，可以參考小編之前的推送宅荤。

如何將轉(zhuǎn)錄組數(shù)據(jù)mapping到自己的序列并可視化屑迂？(HISAT2+Samtools+IGV)

PART4 SAMtools軟件將SAM文件轉(zhuǎn)化為BAM文件

? ? 采用SAMtools軟件（http://samtools.sourceforge.net）進(jìn)行sort及格式轉(zhuǎn)化，sort之后文件占位更小冯键，其轉(zhuǎn)化為BAM二進(jìn)制文件惹盼。

SAMtools安裝：

$ tar jxvf samtools-0.1.19.tar.bz2

$ cd samtools-0.1.19

$ make

$ cd.

$ cp samtools-0.1.19/samtools $HOME/bin

運(yùn)行格式：

$ samtools sort -@ 8 -o Ctrl-1_out.bam Ctrl-1_out.sam

運(yùn)行結(jié)果：

獲得BAM文件，該結(jié)果可采用IGV進(jìn)行可視化惫确。

PART5 采用FeatureCounts軟件進(jìn)行reads計(jì)數(shù)

? ? FeatureCounts是subread軟件包中的一個(gè)工具手报，盡管短序列比對(duì)工具subread沒有bwa和hisat2流行，但其中FeatureCounts工具卻應(yīng)用廣泛改化。其作用是將比對(duì)后的結(jié)果進(jìn)行計(jì)數(shù)掩蛤，即計(jì)算每個(gè)區(qū)域比對(duì)到的reads數(shù)。之后對(duì)reads數(shù)進(jìn)行歸一化處理陈肛，計(jì)算RPKM揍鸟，F(xiàn)PKM，TPM等句旱，然后進(jìn)行差異比較分析阳藻。

Subread安裝（直接找到官網(wǎng)下載解壓即可）:

$ 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

運(yùn)行格式：

$ 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下載。

主要參數(shù)：

-T 多線程數(shù)

-p 針對(duì)雙端測(cè)序數(shù)據(jù)

-t 指定注釋文件中的功能類型前翎，默認(rèn)為外顯子exon

-g 指定注釋文件中的屬性信息稚配，默認(rèn)為gene_id

-a 注釋文件，GTF或GFF格式文件港华，區(qū)分外顯子區(qū)域

-o 輸出結(jié)果文件道川，同時(shí)會(huì)輸出以該名稱結(jié)尾的.summary統(tǒng)計(jì)文件

運(yùn)行結(jié)果：

得到rawcounts文件。之后可在R中進(jìn)行后續(xù)差異比較分析立宜。

參考文獻(xiàn)：

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.