前言

這次給大家?guī)淼氖?6年發(fā)表在NATURE PROTOCOLS上面的一篇處理RNA-seq數(shù)據(jù)的文章：Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown
和它的12年發(fā)表于同一個(gè)雜志的兄弟文章：Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks都是NATURE PROTOCOLS上閱讀量最大的文章。

NATURE PROTOCOLS

當(dāng)然還有很多其它的介紹不止生信的方法文章，大家有時(shí)間可以去探究。
同時(shí)我這里就不再贅述RNA-seq的具體原理，有需要了解的請(qǐng)移步：一個(gè)簡(jiǎn)略的RNA-seq演示
至于軟件的安裝到官網(wǎng)下載，解壓后將bin/添加進(jìn)路徑即可，這里不再做講解。
#######所有操作皆在LINUX&R上完成，默認(rèn)基本處理軟件已經(jīng)安裝

本體介紹

大佬的文章

事實(shí)上作者團(tuán)隊(duì)一直致力于開發(fā)出更好的解決數(shù)據(jù)處理的軟件，就目前12年推出的Tohat2和Cufflinks已經(jīng)不是那么的令人滿意了，所以他們又開發(fā)了 HISAT, StringTie and Ballgown三件套完成大家對(duì)于一個(gè)RNA-seq所有的幻想。#但實(shí)際上目前還存在著其他的性能很優(yōu)秀軟件也可以滿足需求，例如STAR等等，但作為菜鳥來說只有先一步一步的學(xué)習(xí)了。

好了，我們先對(duì)這篇文章給出的處理方法有個(gè)大概的了解。

pipeline：

1. alignment of the reads to the genome
2. assembly the alignments into full-length transcripts
3. quantification of the differencesin expression levels of each gene and transcript
4. calculation of the differences in expression for all genes among the different conditions

An overview of the 'new Tuxedo' protocol

這個(gè)protocol首先從原始RAN-seq數(shù)據(jù)入手，輸出數(shù)據(jù)包括基因list，轉(zhuǎn)錄本，及每個(gè)樣本的表達(dá)量，能夠表現(xiàn)差異表達(dá)基因的表格并完成顯著性的計(jì)算。

1. 第一步是使用HISAT將讀段匹配到參考基因組上,使用者可以提供注釋文件，但HISAT依舊會(huì)檢測(cè)注釋文件沒有列出來的剪切位點(diǎn)
2. 第二步，比對(duì)上的reads將會(huì)被呈遞給StringTie進(jìn)行轉(zhuǎn)錄本組裝，StringTie單獨(dú)的對(duì)每個(gè)樣本進(jìn)行組裝，在組裝的過程中順帶估算每個(gè)基因及isoform的表達(dá)水平
3. 第三步，所有的轉(zhuǎn)錄本都被呈遞給StringTie的merge函數(shù)進(jìn)行merge，這一步是必須的，因?yàn)橛行颖镜霓D(zhuǎn)錄本可能僅僅被部分reads覆蓋，無法被第二步的StringTie組裝出來。merge步驟可以創(chuàng)建出所有樣本里面都有的轉(zhuǎn)錄本，方便下一步的對(duì)比
4. merge的數(shù)據(jù)再一次被呈遞給StringTie，StringTie可以利用merge的數(shù)據(jù)重新估算轉(zhuǎn)錄本的豐度，還能額外的提供轉(zhuǎn)錄本reads數(shù)量的數(shù)據(jù)給下一步的ballgown
5. 最后一步：Ballgown從上一步獲得所有轉(zhuǎn)錄本及其豐度，根據(jù)實(shí)驗(yàn)條件進(jìn)行分類統(tǒng)計(jì)

注意事項(xiàng)Warning：

• HISAT, StringTie, Ballgown并不是唯一處理RNA-seq的方法，也并不能處理所有的RNA-seq數(shù)據(jù)。
• 例如質(zhì)量控制或者去除測(cè)序污染/去除接頭/去除低質(zhì)量數(shù)據(jù)（可以使用FastQC以及FASTXtoolkit做到）
• 這個(gè)pipeline無法處理第三代測(cè)序的數(shù)據(jù)
• Ballgwon可以使用stringtie/cufflinks/RSEM產(chǎn)出的數(shù)據(jù)，但是可能無法接受其它程序產(chǎn)出的結(jié)果
• 默認(rèn)參數(shù)適用于小至三個(gè)，大至小數(shù)百的樣本處理，對(duì)于特殊需求還需要參考manual

數(shù)據(jù)處理設(shè)計(jì)

Read alignment with HISAT

總體上來說HISAT利用了BWA和Bowtie的算法，同時(shí)解決了mRNA中不存在內(nèi)含子難以比對(duì)的問題，比對(duì)上代主流RNA-seq比對(duì)工具能快50倍，同時(shí)需求更少的內(nèi)存<8G（這就意味著你可以在PC上跑數(shù)據(jù)），20個(gè)樣本，每個(gè)樣本一億reads的估計(jì)，用一臺(tái)電腦一天之內(nèi)能夠跑完。使用者可以提供精確的基因注釋來提高在已知基因區(qū)域的準(zhǔn)確性，但這是可選項(xiàng)。

事實(shí)上，針對(duì)RNA-seq數(shù)據(jù)的align目前有很多工具，16年12月12日出版的NATURE METHODS上的一篇文章，利用分別使用人和一種叫malaria的寄生蟲的數(shù)據(jù)，模擬出三種復(fù)雜度的數(shù)據(jù)集（T1、T2和T3）。對(duì)于復(fù)雜度的衡量，定義了三個(gè)尺度：

• 突變率
• indel比例
• 錯(cuò)誤率
  T1復(fù)雜度最低，T3最高。
  隨后使用了目前主流的比對(duì)工具進(jìn)行了比對(duì)。
  因?yàn)镽NA-Seq測(cè)序的特性，天然的會(huì)有一部分?jǐn)?shù)據(jù)延伸到內(nèi)含子區(qū)，這部分跨越外顯子和內(nèi)含子的reads就稱為『junction reads』，所以RNA-Seq比對(duì)軟件需要針對(duì)此進(jìn)行優(yōu)化，而文章做benchmark也考慮到這點(diǎn)。分兩個(gè)水平做比較：
• base and read level，這點(diǎn)和一般的DNA aligner一樣
• junction-level
  度量軟件的結(jié)果時(shí)，用了兩個(gè)值：
• recall： 所有base中正確比對(duì)的比例。Tophat2在T1 base&read的表現(xiàn)，recall是85%-95%，T2降到80%，T3更是雪崩式降至20%以下。這部分表現(xiàn)好的是Novoalign和MapSplice2。

• precision：所有比對(duì)上的base中，正確比對(duì)的比例。對(duì)于T1和T2，大部分軟件這個(gè)值都較高。
  結(jié)果上圖：

  
  
  
  

  接下來看看調(diào)參（自行設(shè)置參數(shù)，而不是使用軟件默認(rèn)參數(shù)）比對(duì)，直接說結(jié)論吧，不管是否調(diào)參，表現(xiàn)robust的是CLC，GSNAP、Novoalign、STAR以及MapSplice2。針對(duì)復(fù)雜度高的T3數(shù)據(jù)，Tophat2調(diào)參和默認(rèn)參數(shù)得到的比對(duì)率相差67%還多。

  
  
  
  另外還有運(yùn)行時(shí)間的比較，這輪表現(xiàn)好的是HISAT/HISAT2，其實(shí)也能看出，數(shù)據(jù)復(fù)雜度越高，比對(duì)耗時(shí)越長(zhǎng)。

Transcript assembly and quantification with StringTie

RNA-seq的分析依賴于精準(zhǔn)的對(duì)于基因isoform的重建以及對(duì)于基因相對(duì)豐度的預(yù)測(cè)。繼承于Cufflinks，StringTie相對(duì)于之前開發(fā)的工具更為準(zhǔn)確，需求內(nèi)存和耗時(shí)也更少。
使用者一樣可以使用注釋文件來幫助StringTie運(yùn)行，對(duì)于低豐度的數(shù)據(jù)比較有幫助。此時(shí)StringTie依舊會(huì)對(duì)非注釋區(qū)域進(jìn)行轉(zhuǎn)錄本的組裝，所以注釋文件在這里是可選選項(xiàng)。

轉(zhuǎn)錄本組裝完成后，使用者可以利用gffcompare(StringTie工具包含)工具來得知有多少轉(zhuǎn)錄本和注釋文件相同，又有多少新的轉(zhuǎn)錄本:

#input: gff or gtf format
transcripts.gtf

command line example:
$ gffcompare -G -r chrX.gtf transcripts.gtf
#-r : reference
#-G :compare all transcripts
#-o :prefix

output file 1: gffcmp.annotated.gtf
# 顯示StringTie組裝的轉(zhuǎn)錄本與注釋文件內(nèi)的轉(zhuǎn)錄本的差別（會(huì)給每個(gè)轉(zhuǎn)錄本加入一個(gè)class code，我理解為一個(gè)標(biāo)識(shí)，釋義如下圖）
output dile 2: gffcmp.stats
# 根據(jù)注釋文件顯示組裝結(jié)果的準(zhǔn)確性和陽(yáng)性預(yù)測(cè)率

class codes

由于在實(shí)驗(yàn)中，我們會(huì)處理多個(gè)RNA-seq數(shù)據(jù)，單個(gè)樣本里面的基因和iosforms與其它樣本的數(shù)據(jù)很少相同，但是為了進(jìn)行比較就需要它們以一個(gè)相同的形式進(jìn)行組裝，作者通過StringTie的merge工具解決了這個(gè)問題，將所有樣本中發(fā)現(xiàn)的基因進(jìn)行merge。下圖的例子可以幫助理解。

example

如圖所示，四個(gè)樣本中組裝得到四個(gè)轉(zhuǎn)錄本，最后merge得到A/B兩個(gè)轉(zhuǎn)錄本。1/2樣本與參考基因組相同merge得到A轉(zhuǎn)錄本，3/4樣本相互吻合但與參考基因組不一樣，所以merge得到B轉(zhuǎn)錄本。
在merge步驟之后，需要StringTie再運(yùn)行一遍去重新估算merge得到的轉(zhuǎn)錄本的豐度。

Chevy理解：這就是相當(dāng)于重新定義了一個(gè)annotation file進(jìn)行二次重新組裝，相對(duì)于第一次組織可以提高準(zhǔn)確度和敏感性。比已有注釋文件的優(yōu)勢(shì)在于可以發(fā)現(xiàn)更多的isoforms。

Differential expression analysis with Ballgown

組裝和定量完轉(zhuǎn)錄本之后自然需要進(jìn)行基因表達(dá)差異分析，統(tǒng)計(jì)建模和可視化。R和Bioconductor提供了一攬子的工具處理這些任務(wù)，包括raw data的plot作圖，數(shù)據(jù)的標(biāo)準(zhǔn)化，下游的統(tǒng)計(jì)建。此處使用Ballgown包作為承上啟下的橋梁。

在R里面使用Ballgown處理需要得到：

• 表型數(shù)據(jù)。關(guān)于樣本的信息
• 表達(dá)數(shù)據(jù)。標(biāo)準(zhǔn)化和未標(biāo)準(zhǔn)化的關(guān)于外顯子，junction，轉(zhuǎn)錄本，基因的表達(dá)數(shù)量
• 基因信息。有關(guān)外顯子，junction，轉(zhuǎn)錄本，基因的坐標(biāo)以及注釋信息

而大多數(shù)的差異表達(dá)分析都會(huì)包括下面幾個(gè)步驟：

• 數(shù)據(jù)可視化和檢查
• 差異表達(dá)的統(tǒng)計(jì)分析
• 多重檢驗(yàn)校正
• 下游檢查和數(shù)據(jù)summary

Ballgown包可以完成以上的幾個(gè)步驟并且可以聯(lián)合R語(yǔ)言的其它操作進(jìn)行分析

Ballgown使用：

• 數(shù)據(jù)的讀入：需要將StringTie輸出的數(shù)據(jù)結(jié)合表型數(shù)據(jù)，這里需要保證表型數(shù)據(jù)的identifier和StringTie輸出數(shù)據(jù)的一致，不然會(huì)報(bào)錯(cuò)
• 預(yù)測(cè)豐度的檢查：以FPKM（fragments per kilobase of transcript per million mapped reads）為單位的豐度預(yù)測(cè)將會(huì)根據(jù)library size進(jìn)行標(biāo)準(zhǔn)化。#極端差異此處需要引起注意
• 使用線性模型進(jìn)行差異表達(dá)分析，由于FPKM對(duì)于轉(zhuǎn)錄本解讀過于曲解，所以這里需要使用log轉(zhuǎn)化處理數(shù)據(jù)，隨后再使用線性模型進(jìn)行差異分析。
• Ballgown可以對(duì)time-course和fixed-condition數(shù)據(jù)進(jìn)行差異分析，但是無法避免批次效應(yīng)帶來的誤差。#使用stattest功能可以指定任何已知的cofounder
• Ballgown可以幫助你在基因，轉(zhuǎn)錄本，外顯子，junction水平上進(jìn)行差異分析
• 結(jié)果會(huì)以表格形式展現(xiàn)，如果樣本夠多還有p值和q值
• 結(jié)果數(shù)據(jù)可以用來進(jìn)行下一步的gene set分析（例如GSEA）等等

實(shí)戰(zhàn)示例

準(zhǔn)備階段

實(shí)際上本輪操作就是按照文章給的示例數(shù)據(jù)走了一遍流程：

文章中，為了減少計(jì)算時(shí)間，方便讀者重復(fù)，作者截取了一批實(shí)驗(yàn)數(shù)據(jù)中能夠匹配到chromosomeX上的數(shù)據(jù)作為示例數(shù)據(jù)，chromosomeX是一個(gè)基因相對(duì)較為豐富的染色體，可以占到基因組的5%左右。

首先是下載數(shù)據(jù)隨后解壓：

$nohup wget ftp://ftp.ccb.jhu.edu/pub/RNAseq_protocol/chrX_data.tar.gz 2>download.log &
$tar zxvf chrX_data.tar.gz

其中包括如下數(shù)據(jù)：

文件數(shù)據(jù)

sample文件夾下有12個(gè)fastq格式的paired-end RNA-seq文件，從兩地人群中（英格蘭島住民和約魯巴住民）各取六個(gè)樣，六個(gè)樣又分為男女性別各三個(gè)（最少的生物學(xué)重復(fù)）。

sample

隨后介紹一個(gè)騷命令：

HISAT2可以直接從NCBI下載sra格式的文件并作為輸入文件格式
下面以ERR188245測(cè)序數(shù)據(jù)為例
$ hisat2 -p 8 --dta -x chrX_data/indexes/chrX_tran --sra-acc ERR188245 –S ERR188245_chrX.sam
可以說是相當(dāng)?shù)膸腿送祽辛?/p>

index文件夾下包含著HISAT2對(duì)于染色體X的index文件

index

genome文件夾下包含這染色體X的序列文件（GRCh38 build 81）

genome

genes文件夾下則包含著針對(duì)基因組的注釋文件（來自于RefSeq數(shù)據(jù)庫(kù)）

annotation

geuvadis_phenodata.csv和mergelist.txt則是用戶著要自己創(chuàng)建表明數(shù)據(jù)關(guān)系的文件，這里作者提供出來方便使用。

如果需要建立官網(wǎng)上沒有基因組的index，就需要需要使用Hisat2包里面的python腳本:
extract_splice_sites.py和extract_exons.py，從注釋文件里面抽取出剪切位點(diǎn)和外顯子信息

$ extract_splice_sites.py chrX_data/genes/chrX.gtf >chrX.ss
$ extract_exons.py chrX_data/genes/chrX.gtf >chrX.exon
先提取剪切位點(diǎn)和外顯子數(shù)據(jù)
$ hisat2-build --ss chrX.ss --exon chrX.exon chrX_data/genome/chrX.fa chrX_tran
隨后建立HISAT2 index

正式開始

1.將reads比對(duì)上genome

2017-07-24更新：我自己使用的是Ensembl發(fā)表的GRCh38版本的基因組進(jìn)行比對(duì)，所以這里我附上這個(gè)版本的genome和注釋文件的下載地址：

Ensembl版本全基因組的注釋文件下載
HISAT2-index下載genome_tran

# 樣本比對(duì)操作
nohup hisat2 -p 8 --dta -x ~/RNA-seq/chrx/chrX_data/indexes/chrX_tran -1 ~/RNA-seq/chrx/chrX_data/samples/ERRERR188044_chrX_1.fastq.gz -2 ERR188044_chrX_2.fastq.gz -S ERR188044_chrX.sam &

# 數(shù)據(jù)太多我就寫了個(gè)小腳本處理,in sample dir/
for i in *1.fastq.gz; 
do
i=${i%1.fastq.gz*}; 
nohup hisat2 -p 8 --dta -x ~/RNA-seq/chrx/chrX_data/indexes/chrX_tran -1 ${i}1.fastq.gz -2 ${i}2.fastq.gz -S ~/RNA-seq/chrx/chrX_data/result/${i}align.sam 2>~/RNA-seq/chrx/chrX_data/result/${i}align.log & 
done

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

比對(duì)結(jié)果

2.將sam文件轉(zhuǎn)化為bam文件

# 轉(zhuǎn)化操作
samtools sort @ 8 -o ERR188044_chrX.bam ERR188044_chrX.sam

# 批處理, in result dir/
for i in *.sam;
do 
i=${i%_align.sam*}; nohup samtools sort -@ 8 -o ${i}.bam ${i}_align.sam & 
done

結(jié)果如下：

sam to bam

3.組裝轉(zhuǎn)錄本并定量表達(dá)基因

# 單文件操作
stringtie -p 8 -G ~/RNA-seq/chrx/chrX_data/genes/chrX.gtf -o ERR188044_chrX.gtf ERR188044_chrX.bam

# 批處理, in result dir/
for i in *.bam; 
do 
i=${i%.bam*}; nohup stringtie -p 8 -G ~/RNA-seq/chrx/chrX_data/genes/chrX.gtf -o ${i}.gtf ${i}.bam & 
done

結(jié)果如下:

strigtie

4.將所有轉(zhuǎn)錄本合并

warning: 此處的mergelist.txt是自己創(chuàng)建的，需要包含之前output.gtf文件的路徑

cat ~/RNA-seq/chrx/chrX_data/mergelist.txt
ERR188044_chrX.gtf
ERR188104_chrX.gtf
ERR188234_chrX.gtf
ERR188245_chrX.gtf
ERR188257_chrX.gtf
ERR188273_chrX.gtf
ERR188337_chrX.gtf
ERR188383_chrX.gtf
ERR188401_chrX.gtf
ERR188428_chrX.gtf
ERR188454_chrX.gtf
ERR204916_chrX.gtf
#因?yàn)榫驮诮Y(jié)果目錄下面操作，所以直接顯示文件名即可

stringtie --merge -p 8 -G  ~/RNA-seq/chrx/chrX_data/genes/chrX.gtf -o stringtie_merged.gtf  ~/RNA-seq/chrx/chrX_data/mergelist.txt

#output文件即為
stringtie_merged.gtf 

5.檢測(cè)相對(duì)于注釋基因組，轉(zhuǎn)錄本的組裝情況

gffcompare -r ~/RNA-seq/chrx/chrX_data/genes/chrX.gtf  -G -o merged stringtie_merged.gtf
#輸出文件信息可見上面的Transcript assembly and quantification with StringTie介紹

6.重新組裝轉(zhuǎn)錄本并估算基因表達(dá)豐度，并為ballgown創(chuàng)建讀入文件

# 單文件操作
stringtie -e -B -p 8 -G stringtie_merged.gtf -o ballgown/ERR188044_chrX/ERR188044_chrX.gtf ERR188044_chrX.bam

# 批處理, in result dir/
for i in *_chrX.bam; 
do 
i=${i%_chrX.bam*}; nohup stringtie -e -B -p 8 -G stringtie_merged.gtf -o ballgown/${i}/${i}_chrX.gtf  ${i}_chrX.bam & 
done

結(jié)果如下：

table

7.R包的安裝與載入

R語(yǔ)言需要安裝ballgown包和一些接下來處理會(huì)用到的包（RSkittleBrewer/genefilter/dplyr/devtools）
事實(shí)上我還沒搞清楚LINUX里面的R包安裝問題，老是報(bào)error，所以這里我就用Windows下的R進(jìn)行操作。

# 首先是幾個(gè)包的載入
>library(ballgown)
>library(RSkittleBrewer)
>library(genefilter)
>library(dplyr)
>library(devtools)

8.數(shù)據(jù)的讀入

# 隨后是數(shù)據(jù)的讀入，CSV文件，我把所有文件放在桌面，上一步得到的ballgown文件夾直接拷到桌面
> setwd("C:/Users/DELL/Desktop")
> read.csv("geuvadis_phenodata.csv")
         ids    sex population
1  ERR188044   male        YRI
2  ERR188104   male        YRI
3  ERR188234 female        YRI
4  ERR188245 female        GBR
5  ERR188257   male        GBR
6  ERR188273 female        YRI
7  ERR188337 female        GBR
8  ERR188383   male        GBR
9  ERR188401   male        GBR
10 ERR188428 female        GBR
11 ERR188454   male        YRI
12 ERR204916 female        YRI
> pheno_data<-read.csv("geuvadis_phenodata.csv")
> bg_chrX=ballgown(dataDir = "C:/Users/DELL/Desktop/ballgown",samplePattern = "ERR",pData = pheno_data)
# dataDir指的是數(shù)據(jù)儲(chǔ)存的地方，samplePattern則依據(jù)樣本的名字來，pheno_data則指明了樣本數(shù)據(jù)的關(guān)系，這個(gè)里面第一列樣本名需要和ballgown下面的文件夾的樣本名一樣，不然會(huì)報(bào)錯(cuò)

10.濾掉低豐度的基因

# 這里濾掉了樣本間差異少于一個(gè)轉(zhuǎn)錄本的數(shù)據(jù)
> bg_chrX_filt=subset(bg_chrX,"rowVars(texpr(bg_chrX))>1",genomesubset=TRUE)

11.確認(rèn)組間有差異的轉(zhuǎn)錄本

2017-07-21更新：實(shí)際上為了確定基因表達(dá)差異是由于某個(gè)變量造成的，我們這邊需要使用剩下的變量進(jìn)行修正。Example：我們?yōu)榱吮容^male and female之間的基因表達(dá)差異，這里就需要指定分析參數(shù)為"transcripts"，主變量為"sex"，修正變量為"population"，getFC可以指定輸出結(jié)果顯示組間表達(dá)量的foldchange。

并且Ballgown的統(tǒng)計(jì)算法基于標(biāo)準(zhǔn)的線性模型，對(duì)于組內(nèi)數(shù)據(jù)較少(<4)時(shí)較為準(zhǔn)確。這里也可以使用其它的一些軟件例如limma/DESeq/edgeR等。

> result_transcripts=stattest(bg_chrX_filt,feature = "transcript",covariate = "sex",adjustvars = c("population"), getFC=TRUE,meas="FPKM")
> result_transcripts
        feature   id           fc         pval         qval
1    transcript    1  0.941367186 7.353184e-01 9.468599e-01
2    transcript    2  1.207949354 8.668638e-01 9.744055e-01
3    transcript    3  1.013100019 9.920824e-01 9.986659e-01
4    transcript    4  0.387671042 5.233364e-01 9.189906e-01
5    transcript    5  0.615797877 3.324164e-01 9.117015e-01
......

12.確認(rèn)組間有差異的基因

這里指定feature為gene

> result_genes=stattest(bg_chrX_filt,feature = "gene",covariate = "sex",adjustvars = c("population"), getFC=TRUE,meas="FPKM")
> result_genes
     feature           id          fc         pval         qval
1       gene      MSTRG.1 1.114999374 7.305975e-01 9.218157e-01
2       gene     MSTRG.10 1.749747142 2.767538e-01 7.922755e-01
3       gene   MSTRG.1000 1.399358968 6.039065e-01 8.945639e-01
4       gene   MSTRG.1002 0.937668743 8.825216e-01 9.816878e-01
5       gene   MSTRG.1003 1.044840944 6.155345e-01 8.977043e-01
6       gene   MSTRG.1004 1.220626116 4.160214e-01 8.388688e-01
.....

13.對(duì)結(jié)果 result_transcripts增加基因名

> result_transcripts=data.frame(geneNames=ballgown::geneNames(bg_chrX_filt),geneIDs=ballgown::geneI Ds(bg_chrX_filt),result_transcripts)
> result_transcripts
          geneNames      geneIDs    feature   id           fc         pval         qval
1                 .      MSTRG.4 transcript    1  0.941367186 7.353184e-01 9.468599e-01
2            PLCXD1      MSTRG.4 transcript    2  1.207949354 8.668638e-01 9.744055e-01
3                 .      MSTRG.4 transcript    3  1.013100019 9.920824e-01 9.986659e-01
4                 .      MSTRG.4 transcript    4  0.387671042 5.233364e-01 9.189906e-01
5                 .      MSTRG.5 transcript    5  0.615797877 3.324164e-01 9.117015e-01
6            PLCXD1      MSTRG.4 transcript    6  0.630018786 2.938396e-01 9.002815e-01
......

14.按照P值排序（從小到大）

> result_transcripts=arrange(result_transcripts,pval)
> result_genes=arrange(result_genes,pval)

15.把結(jié)果寫入csv文件

> write.csv(result_transcripts, "chrX_transcript_results.csv",row.names=FALSE)
> write.csv(result_genes, "chrX_gene_results.csv",row.names=FALSE)

16.篩選出q值小于0.05的transcripts和genes，也就是在性別間表達(dá)有差異的基因了。

注：這里無需過分關(guān)注geneIDs以及transcript name，實(shí)際上那個(gè)是stringtie在比對(duì)過程中賦予的一個(gè)符號(hào)。

> subset(result_transcripts,result_transcripts$qval<0.05)
   geneNames   geneIDs    feature   id          fc         pval         qval
1          . MSTRG.531 transcript 1657 0.030820591 1.427864e-10 2.082377e-07
2       XIST MSTRG.531 transcript 1656 0.003014576 1.860927e-10 2.082377e-07
3          . MSTRG.531 transcript 1655 0.016144997 3.762791e-10 2.807042e-07
4          . MSTRG.531 transcript 1658 0.028308029 6.599039e-08 3.692163e-05
5       TSIX MSTRG.530 transcript 1654 0.078461818 1.690716e-06 7.567643e-04
6          . MSTRG.613 transcript 1848 7.391342987 1.249275e-05 4.659796e-03
7          . MSTRG.141 transcript  421 3.204058932 2.729898e-05 8.727874e-03
8          . MSTRG.618 transcript 1852 9.136857610 4.804244e-05 1.343987e-02
9          . MSTRG.777 transcript 2338 0.127674288 1.000751e-04 2.488533e-02
10     KDM6A MSTRG.258 transcript  736 0.054212867 1.173824e-04 2.627018e-02
11    PNPLA4  MSTRG.62 transcript  186 0.592778584 2.083496e-04 4.238966e-02
12     RPS4X MSTRG.519 transcript 1611 0.597532121 2.493976e-04 4.651265e-02

> subset(result_genes,result_genes$qval<0.05)
   feature        id          fc         pval         qval
1     gene MSTRG.531 0.002396124 2.469125e-11 2.523446e-08
2     gene MSTRG.141 3.165966412 1.666206e-06 8.514310e-04
3     gene MSTRG.530 0.082749314 7.018439e-06 2.390948e-03
4     gene MSTRG.613 7.314423295 1.128589e-05 2.883544e-03
5     gene MSTRG.618 9.050399157 5.022017e-05 1.026500e-02
6     gene MSTRG.519 0.637382385 6.953432e-05 1.184401e-02
7     gene MSTRG.376 0.620693674 1.134561e-04 1.656458e-02
8     gene  MSTRG.62 0.600686117 1.638615e-04 2.093330e-02
9     gene MSTRG.777 0.159178288 2.177206e-04 2.472338e-02
10    gene MSTRG.622 7.870447732 3.737270e-04 3.527295e-02
11    gene MSTRG.778 0.127712244 3.796501e-04 3.527295e-02
12    gene MSTRG.229 1.412379185 4.200674e-04 3.577574e-02
13    gene  MSTRG.48 4.158987103 5.279898e-04 4.150812e-02

17.數(shù)據(jù)可視化之顏色設(shè)定

# 賦予調(diào)色板五個(gè)指定顏色
> tropical= c('darkorange', 'dodgerblue','hotpink', 'limegreen', 'yellow')
> palette(tropical)

# 當(dāng)然rainbow()函數(shù)也可以完成這個(gè)任務(wù)
> palette(rainbow(5))

18.以FPKM為參考值作圖，以性別作為區(qū)分條件

# 提取FPKM值
> fpkm = texpr(bg_chrX,meas="FPKM")

#方便作圖將其log轉(zhuǎn)換，+1是為了避免出現(xiàn)log2(0)的情況
> fpkm = log2(fpkm+1)

# 作圖
> boxplot(fpkm,col=as.numeric(pheno_data$sex),las=2,ylab='log2(FPKM+1)')

結(jié)果如下：

FPKM, male in blue, females in orange

19.就單個(gè)轉(zhuǎn)錄本的查看在樣品中的分布

> ballgown::transcriptNames(bg_chrX)[12]
         12 
"NR_027232" 
> ballgown::geneNames(bg_chrX)[12]
         12 
"LINC00685" 

# 繪制箱線圖
> plot(fpkm[12,] ~ pheno_data$sex, border=c(1,2),
+      main=paste(ballgown::geneNames(bg_chrX)[12],' : ',
+                 ballgown::transcriptNames(bg_chrX)[12]),pch=19, xlab="Sex",
+      ylab='log2(FPKM+1)')
> points(fpkm[12,] ~ jitter(as.numeric(pheno_data$sex)),
+        col=as.numeric(pheno_data$sex))

結(jié)果如下：

transcript

20.查看某一基因位置上所有的轉(zhuǎn)錄本

# plotTranscripts函數(shù)可以根據(jù)指定基因的id畫出在特定區(qū)段的轉(zhuǎn)錄本
# 可以通過sample函數(shù)指定看在某個(gè)樣本中的表達(dá)情況，這里選用id=1750, sample=ERR188234
> plotTranscripts(ballgown::geneIDs(bg_chrX)[1750], bg_chrX, main=c('Gene MSTRG.575 in sample ERR188234'), sample=c('ERR188234'))

結(jié)果如下：

21.以性別為區(qū)分查看基因表達(dá)情況

# 這里以id=575的基因?yàn)槔▽?duì)應(yīng)上一步作圖）
> plotMeans('MSTRG.575', bg_chrX_filt,groupvar="sex",legend=FALSE)

結(jié)果如下：

MSTRG.575

臨終討論

• 關(guān)于時(shí)間的使用：針對(duì)chrX的數(shù)據(jù)量的比對(duì)和組裝，在PC上可以40min內(nèi)搞定，可以說是比較快的了。作者論述如果將12個(gè)樣本的全數(shù)據(jù)下下來做分析的話，8核&8G內(nèi)存的配置花費(fèi)12.5h可以搞定比對(duì)，花費(fèi)8h可以搞定組裝和表達(dá)分析。

• RNA_seq的比對(duì)情況一覽

  
  

• 轉(zhuǎn)錄組組裝情況一覽

  
  

• 其實(shí)這篇文章說到底還是在講方法論，它只負(fù)責(zé)上游的數(shù)據(jù)處理到可以分析的一步，下游的分析和可視化還需要依賴自身的實(shí)驗(yàn)設(shè)計(jì)和研究目的。

• 官方給出的方案是：stringtie處理得到的數(shù)據(jù)是直接呈遞給ballgown進(jìn)行差異分析和可視化作圖，但是我還不是很清楚輸出文件是否可以無差別的被其它軟件使用。（我也是菜鳥）

• RNA-seq的文章其實(shí)很多，也不是非要走這個(gè)pipeline。另外我覺得ballgown的作圖其實(shí)也沒有很elegant，org

參考文獻(xiàn)

Pertea M, Kim D, Pertea G M, et al. Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown[J]. Nature protocols, 2016, 11(9): 1650-1667.

日常Bob鎮(zhèn)樓