文章題目：Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown
利用Hisat、Stringtie、Ballgown進行RNA測序中轉(zhuǎn)錄本定量
來源：Nature Protocols DOI：10.1038/nprot.2016.095
第一部分：詞句經(jīng)典
- High-throughput sequencing of mRNA (RNA-seq) has become the standard method for measuring and comparing the levels of gene expression in a wide variety of species and conditions.
【一句很通用的話介紹了一下NGS對于RNA-seq的意義】 - RNA-seq experiments capture the total mRNA from a collection of cells and then sequence that RNA in order to determine which genes were active, or expressed, in those cells.
【說明了RNA-seq的作用 p.s. 本文做的都是mRNA的研究】 - generate enormous numbers of raw sequencing reads, typically numbering in the tens of millions,
- three is the minimum number of replicates for valid statistical results
第二部分：文章結(jié)構(gòu)
摘要前言：
RNA測序能夠在一次試驗中捕獲成千上萬基因的表達(dá)量，能產(chǎn)生原始reads數(shù)十M，通過檢測每個基因測出的reads數(shù)目可以估算基因豐度，當(dāng)有兩個或多個重復(fù)時能夠判斷差異表達(dá)基因。另外利用RNA-seq可以檢測注釋文件中沒有的變異基因，也可以尋找每個基因不同的調(diào)控表達(dá)途徑。
分析RNA-seq一般需要四步:
- 將reads比對到參考基因組上（這篇文章是以有參轉(zhuǎn)錄組為例）；
- 將比對上的alignments組裝成全長轉(zhuǎn)錄本；
- 基因與轉(zhuǎn)錄本的表達(dá)定量；
- 計算不同實驗條件下，所有基因的表達(dá)差異
之前有一套公認(rèn)的好用又快的處理轉(zhuǎn)錄組的工具Tophat2和Cufflinks是由約翰霍普金斯大學(xué)團隊開發(fā)，后來被被創(chuàng)作團隊淘汰，當(dāng)大家不理解的時候，他們發(fā)聲：別用之前的啦，不好用！不夠快！～我們有了新的三件套！于是新版三件套工具就這么出來了，還命名為“Tuxedo ”～～真是應(yīng)了一句話無敵是多么寂寞
Hisat相對于Tophat2比對和發(fā)現(xiàn)可變剪切位點速度更快，消耗內(nèi)存更少；【可以提供gff基因注釋文件，當(dāng)然如果沒有程序也會自動檢測剪切位點】
StringTie負(fù)責(zé)拼接轉(zhuǎn)錄本，構(gòu)建isoforms用于估算基因表達(dá)量；【先接收Histat傳來的alignments進行轉(zhuǎn)錄本拼接，每組數(shù)據(jù)之間是獨立的，拼接的過程中估算每個gene和每個isoform的表達(dá)量；拼接完以后，所有的拼接好的序列進行merge，這一步是必須的，因為某些樣本中的轉(zhuǎn)錄本只是被reads部分覆蓋，導(dǎo)致的結(jié)果就是僅有那些被覆蓋的區(qū)域得到了拼接，利用merge可以降低這一誤差?！縨erge完成的轉(zhuǎn)錄本重新返回給StringTie，計算轉(zhuǎn)錄本豐度，它使用的算法和一開始拼接的一樣，但如果在merge過程中轉(zhuǎn)錄本的豐度雖然不變，但結(jié)構(gòu)發(fā)生了改變，reads也要因此進行調(diào)整。StringTie還提供了對每個轉(zhuǎn)錄本進行read-count，這個數(shù)據(jù)會傳給BallGown
BallGown利用StringTie拼接的結(jié)果，計算得出差異表達(dá)轉(zhuǎn)錄本【接受StringTie傳來的轉(zhuǎn)錄本和定量數(shù)據(jù)，根據(jù)設(shè)置的不同的實驗條件進行分組】它包含了某些R包可以直接可視化
這三件套流程支持帶有時間線的實驗(比如研究某些病的發(fā)育歷期)以及兩個以上實驗條件下的結(jié)果比較
流程很容易操作，但需要會使用命令行及R腳本運行
當(dāng)然，這三個工具并不是缺一不可：如果你是要做無參轉(zhuǎn)錄組，可以用其他工具如trinity拼接好以后傳給StringTie；也有一些項目只需要驗證已知的某些基因在樣本中的表達(dá)量；另外即使對于模式物種，基因注釋文件也不完善，因此可以結(jié)合許多其他工具使用。雖然這三種工具能夠檢測差異基因表達(dá)，但對于外顯子的差異不能探索，好在可以結(jié)合使用其他工具如：DEXseq 、rMATS 、MISO 。
這三件套不包括數(shù)據(jù)預(yù)處理（去接頭、去污染、去低質(zhì)量序列），但可以用第三方軟件FASTX (http://hannonlab.cshl.edu/fastx_toolkit) 、FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc) ；
適用于二代illumina測序，對于三代Pac Bio、Nanopore的比對環(huán)節(jié)需要用其他工具；
Ballgown目前支持三個組裝軟件的結(jié)果：StringTie、Cufflinks、RSEM
本文帶了實例數(shù)據(jù)：
人類的RNA-seq數(shù)據(jù)，選取了其中基因比較豐富的X染色體數(shù)據(jù)151M，相當(dāng)于全基因組的5%
技術(shù)路線如下：
具體方法：
-
HISAT（Hierarchical Indexing for Spliced Alignment of Transcripts）比對
比對的工具目前常用的有Bowtie2、BWA，他們運行時使用了BWT（Burrows-Wheeler transform）的數(shù)據(jù)格式，這是一種壓縮格式，能將參考基因組高度壓縮。然后再使用特殊的index機制——Ferragina–Manzini (FM) index，因此能快速在基因組上搜索匹配reads的區(qū)段，速度高達(dá)每小時比對幾百萬條reads。以上方法在DNA比對中最為常用，例如組裝基因組需要比對拼接。但是對于RNA來講，有一個DNA中不存在的問題需要注意：成熟的mRNA將introns切除，因此許多RNA測序的reads要跨越內(nèi)含子進行比對，因此一條短序列可能會被拆開比對到相隔1萬bp甚至更遠(yuǎn)的位置 【人類平均內(nèi)含子平均長度大于6000bp，長的有超過1M的】人類RNA測序采用的得到的是100bp reads，會有超過35%的reads能跨越式比對到多個外顯子。
要將一條跨越多個外顯子的reads比對到基因組上，其難度要大于 本就是在一個外顯子中的reads。Hisat使用了疊加的比對算法：一是全局比對，整個基因組建立index；輔助成千上萬個小的局部index；
這兩種算法也是基于BWT/FM體系建立起來的。因此這種特制的比對算法使用時比BWA、bowtie要快好幾倍，但是內(nèi)存用量僅是他們的兩倍。相比較于他的前任Tophat2，速度提升了50倍。如果比對人類參考基因組（3G大小）的話，內(nèi)存需要最少8G，比較人性化，一般電腦就能跑起來；
8G內(nèi)存，20個樣本，每個樣本100Mreads，一天能跑完比對環(huán)節(jié)。
用戶可以在比對的時候，可以添加基因注釋--描述已知基因的位置（包括它們外顯子、內(nèi)含子邊界），這樣比對會增加結(jié)果準(zhǔn)確性，可以更快找到新的剪切位點、轉(zhuǎn)錄起始和終止位點 -
StringTie拼接、融合
首先將比對完的reads劃分成不同的基因座組，然后把每個基因座拼接成isoform。運用了network flow algorithm ，從reads數(shù)最集中(也就是表達(dá)量最高)的轉(zhuǎn)錄本開始，然后移除已拼接完的reads，再次重復(fù)這個過程，最后拼接成眾多的isoform，結(jié)束的標(biāo)志是：全部的reads都被拼接完或者剩下的reads數(shù)目低于某個限制；
Network flow algorithm這個算法能夠同時計算轉(zhuǎn)錄本reads豐度和外顯子-內(nèi)含子結(jié)構(gòu)，因此它重構(gòu)各個轉(zhuǎn)錄本的速度比以前版本Cufflinks更快，消耗內(nèi)存更少，但最少要求內(nèi)存8G。如果用戶提供了注釋文件，stringtie會以它為指導(dǎo)進行拼接，另外還會將拼接的基因根據(jù)已注釋好的直接命名。注釋文件對低表達(dá)量的基因很有幫助，因為reads數(shù)量太少，不好拼接得到該基因。拼接完后，可以用
gffcompare檢測有多少拼接的轉(zhuǎn)錄本或全部或部分地map到了已知的基因，并且統(tǒng)計有多少是全新的基因。【此環(huán)節(jié)可以跳過】想要比較不同樣本基因表達(dá)量差異，如果通過一個一個基因?qū)ふ疫M行比較是很難的，但是如果將每個樣本中所有的基因融合到一起，再比較樣本的差異就方便多。StringTie提供了
merge的功能，將各個相似的樣本中的所有基因融合起來，再比較融合后的大序列不同。即使某一個樣本中缺少了某一段外顯子，通過和其他樣本的比較參考，也能將這段缺少外顯子序列找到和它相似的進行merge，**最根本的目的就是：相似的拼接alignments進行merge，再進行比較。merge這個圖就解釋了merge的用處：綠色的sample1-4是來自不同的四個樣本的部分拼接片段；黑色的是基因組注釋文件的部分；sample1和2都和參考基因組一致，所以把他們merge到一起形成了轉(zhuǎn)錄本A；sample3和4與參考基因組不一致，但彼此之間是一致的，所以把他們merge到一起形成轉(zhuǎn)錄本B
merge完成以后，會重新估算一遍merged transcripts表達(dá)量。
StringTie可以不依賴于注釋文件進行新的基因、轉(zhuǎn)錄本預(yù)測，注釋文件只是個佐證
-
Ballgown差異表達(dá)分析
這算是上游Hisat、StringTie與下游R/bioconductor的中間橋梁，他將數(shù)據(jù)進行標(biāo)準(zhǔn)化
基本上差異表達(dá)分析都會包括下面幾個步驟：(i) data visualization and inspection, (ii) statistical tests for differential expression, (iii) multiple test correction and (iv) downstream inspection and summarization of results.- 數(shù)據(jù)可視化和檢查
- 差異表達(dá)的統(tǒng)計分析
- 多重檢驗校正
- 下游檢查和數(shù)據(jù)summary
在R里面使用Ballgown處理需要得到：
- 表型數(shù)據(jù)：關(guān)于樣本的信息
- 表達(dá)數(shù)據(jù)：標(biāo)準(zhǔn)化和未標(biāo)準(zhǔn)化的外顯子，轉(zhuǎn)錄本，基因的表達(dá)數(shù)量
- 基因信息：有關(guān)外顯子，轉(zhuǎn)錄本，基因的坐標(biāo)以及注釋信息
使用Ballgown：
- 讀入數(shù)據(jù)：將上游Stringtie輸出的轉(zhuǎn)錄本表達(dá)量數(shù)據(jù)與表型數(shù)據(jù)結(jié)合起來【注意保證二者的ID號一致】
- 豐度預(yù)測：以FPKM（fragments per kilobase of transcript per million mapped reads）為單位的豐度預(yù)測根據(jù)library size進行標(biāo)準(zhǔn)化。
- 差異表達(dá)分析：FPKM的數(shù)據(jù)解讀需要取log轉(zhuǎn)化一下，再建立線性模型
- 可以在基因，轉(zhuǎn)錄本，外顯子水平上進行差異分析
- 結(jié)果以表格形式展現(xiàn)，樣本量大的話還有p值和q值
第三部分：數(shù)據(jù)測試
軟件準(zhǔn)備：
HISAT2 software (http://ccb.jhu.edu/software/hisat2 or http://github.com/infphilo/hisat2, version 2.0.1 or later) 【hisat2支持--sra-acc從NCBI SRA數(shù)據(jù)庫直接下載數(shù)據(jù)，然后自動轉(zhuǎn)為fa/fq格式】
StringTie software (http://ccb.jhu.edu/software/stringtie or https://github.com/gpertea/stringtie , version 1.2.2 or later)
SAMtools (http://samtools.sourceforge.net, version 0.1.19 or later)
數(shù)據(jù)準(zhǔn)備：
下載RNA-seq測試數(shù)據(jù)
其中包括人類基因組X染色體RNA-seq，基因注釋文件 ?！綳染色體參考基因組、測試RNA-seq數(shù)據(jù)都在其中】
-
解壓chrX_data.tar.gz會得到四個文件夾：samples、indexes、genome、genes
samples中包含2個種群(GBR (British from England) and YRI (Yoruba from Ibadan, Nigeria) )各6個sample（男女各3個重復(fù)），一共12個sample。然后使用雙端測序，結(jié)果共有24個測序fq文件
測序文件indexes 中是hisat2構(gòu)建好的chrX索引文件；
genome 中是參考基因組chrX.fa (GRCh38 build 81 );
genes 中是chrX.gtf , 是從RefSeq數(shù)據(jù)庫中下載的GRCh38的基因組注釋文件 -
hisat2構(gòu)建索引
我們其實下載好了indexes文件，但是如果從頭開始自己建立索引應(yīng)該怎么做呢？
新建一個文件夾index，用于存放一會生成的chrX_tran索引文件- 第一步利用hisat2的兩個python程序,將剪切位點、外顯子區(qū)域找出來：
extract_splice_sites.py chrX_data/genes/chrX.gtf >chrX.ss
extract_exons.py chrX_data/genes/chrX.gtf >chrX.exon - 第二步使用hisat2-build
hisat2-build --ss chrX.ss --exon chrX.exon chrX_data/genome/chrX.fa chrX_tran
【如果注釋文件不存在，--ss、--exon可以省略】
這一步僅構(gòu)建chrX的索引就需要內(nèi)存9G左右，如果針對全基因組構(gòu)建需要160G內(nèi)存以上。
沒有注釋文件就不用這么多內(nèi)存，另外對chrX構(gòu)建索引，使用一個CPU核心差不多10min；
構(gòu)建全基因組的使用8個核心，大概需要2h
我這里運行的時間是00:16:43，可能服務(wù)器有其他程序在跑
- 第一步利用hisat2的兩個python程序,將剪切位點、外顯子區(qū)域找出來：
-
比對環(huán)節(jié)
原來文章中是一個一個樣本運行的，這樣很麻煩，所以我寫了一個小腳本，提高效率-
將每個sample的reads比對到基因組上：
基本的程序使用是這樣的：hisat2 -p 24 --dta -x index/chrX_tran -1 chrX_data/samples/chrX_data/samples/ERR188273_chrX_1.fastq.gz -2 chrX_data/samples/ERR188273_chrX_2.fastq.gz -S ERR188273_chrX.samhisat2有幾個選項設(shè)置：
-p :線程數(shù)
--dta：輸出和拼接工具兼容的alignments(此處默認(rèn)Stringtie)
--dta-cufflinks:如果用cufflinks拼接，就要用這個選項
--known-splicesite-infile ：如果之前沒有建立剪切位點索引，直接用下載的gtf作為索引是不行的，可以用這個選項添加剪切位點文件【一般還是建議按上一步提前建立好索引文件】
-x : 索引文件【這里就有一個問題：因為有12組24個文件，所以我不可能一個一個去運行，這里寫一個循環(huán)腳本，讓程序自動運行】
【問題的關(guān)鍵是如何讓輸入的-1、-2雙端測序文件匹配到12個sample數(shù)據(jù)上？】觀察發(fā)現(xiàn)這些文件的命名都是有規(guī)律的，不同的只是前邊數(shù)字部分，如果我們可以把1.fastq.gz之前的部分（ERRxxxxxx_chrX_）提取出來，那么運行程序的時候只需要引用名字，再加上對應(yīng)的1.fastq.gz或者2.fastq.gz就可以啦！但是還要注意的是提取的時候只需要用1.fastq或者2.fastq其中一類就行，否則提取出來的會有重復(fù)
文件命名【腳本如下】 for i in *1.fastq.gz;do names=${i%_*} #上面??這一句的意思是：匹配提取出文件_前面的內(nèi)容（或許你會問：為什么數(shù)字編號后面也有_，卻不會匹配到前面的數(shù)字呢？這是因為存在一個貪婪匹配原則，匹配越多越好） hisat2 -p 24 --dta -x ../../index/chrX_tran -1 ${names}_1.fastq.gz -2 ${name s}_2.fastq.gz -S ../../align/${names}.sam done -
用samtools進行排序和轉(zhuǎn)換：
vi sam2bam.sh for i in *.sam;do names=${i%.*} samtools sort -@ 8 -o ${names}.bam ${names}.sam done time sh sam2bam.sh 【運行時間: 4m33.061s】 -
Stringtie拼接轉(zhuǎn)錄本：
vi assemble.sh for i in *.bam;do names=${i%.*} stringtie $i -p 8 -G ../chrX_data/genes/chrX.gtf -o ${names}.gtf done 【參數(shù)：-G：參考注釋文件；-p：線程數(shù)；-o：輸出的合并轉(zhuǎn)錄本GTF格式】 【運行時間：1m53.666s】 -
Stringtie轉(zhuǎn)錄本合并：
標(biāo)準(zhǔn)格式：stringtie --merge [Options] { gtf_list | strg1.gtf ...}先構(gòu)建merge_list: 其實就是把所有的輸出gtf文件放到一個文本文檔里，方便merge程序調(diào)取，
當(dāng)然如果merge的轉(zhuǎn)錄本量少的話，也可以一個一個輸入vi mergelist.sh for i in *.gtf;do echo $i > mergelist.txt done sh mergelist.sh > mergelist.txt運行merge： stringtie --merge -p 8 -G ../chrX_data/genes/chrX.gtf -o stringtie_merged.gtf mergelist.txt 【運行時間：0m2.597s 】 【查看結(jié)果有31027行】 -
檢查拼接的轉(zhuǎn)錄本有多少匹配到了參考基因組注釋文件【可選】
使用gffcompare
安裝：http://ccb.jhu.edu/software/stringtie/gff.shtml or http://github.com/gpertea/gffcompare標(biāo)準(zhǔn)使用：
gffcompare –r chrX.gtf –G transcripts.gtf
-r：參考注釋文件
-G：拼接的轉(zhuǎn)錄本，也是要被比較的對象
-o：指定輸出的結(jié)果前綴，否則默認(rèn)gffcmp這里使用的程序是：
gffcompare –r ../chrX_data/genes/chrX.gtf –G –o merged stringtie_merged.gtf
匹配的結(jié)果存放在merged.annotated.gtf中，其中在附加列中有一項內(nèi)容"class code"
可以幫助我們快速判斷轉(zhuǎn)錄本與參考基因組的關(guān)系
class code還有一個文件
gffcmp.stats,他計算了不同的gene features (e.g., exons, introns, transcripts and genes) 的Sensitivity和Precision分析數(shù)據(jù)
【其中，Sensitivity是指：參考注釋文件中的gene在拼接轉(zhuǎn)錄本中正確重構(gòu)的比例；
Precision：拼接轉(zhuǎn)錄本與參考注釋重疊的比例】
還計算了拼接轉(zhuǎn)錄本中新外顯子、內(nèi)含子、基因座的數(shù)量!其他信息
-
-
Stringtie估算轉(zhuǎn)錄本表達(dá)量并生成Ballgown能使用的統(tǒng)計表格
首先新建一個文件夾ballgown，并在其子目錄下新建sample名的文件夾，一會用來存放各個統(tǒng)計數(shù)據(jù),像這樣：
Ballgown文件夾
【如果你要一個個mkdir敲進去肯定很慢，而且還容易錯】
【批量新建文件夾簡單的辦法： 把文件夾名都放在一個文本文檔中，這里我們先用mergelist.txt復(fù)制過來就行，但是mergelist.txt不能直接用，因為它里面只有前半部分是我們需要的】
【本著能用一行命令絕不寫腳本的原則：
cat mergelist.txt | sed 's/_.*f//' |xargs -n1 mkdir】
直接按我們的需求新建好了文件夾
> 接下來運行stringtie
> vi estimate_abundance.sh
>
> for i in *.bam;do
> names=${i%_*}
> stringtie -e -B -p 8 -G stringtie_merged.gtf -o ../ballgown/$names/${names}_chrX.gtf $i
> done
>
> -e:只估算給定的參考轉(zhuǎn)錄本（就是合并之后的）【與-G搭配】
> -G：參考注釋文件【這里用stringtie_merged.gtf】
> -B：生成Ballgown格式的表格，和output文件放一起【與-o搭配】
> 運行結(jié)果如下：
stringtie運行結(jié)果
-
Ballgown差異表達(dá)分析
7.1首先下載并加載R包
下載：“ballgown” source("https://bioconductor.org/biocLite.R") biocLite("ballgown") “devtools” install.packages("devtools") "genefilter" source("https://bioconductor.org/biocLite.R") biocLite("genefilter") "RSkittleBrewer" devtools::install_github('RSkittleBrewer', 'alyssafrazee') 這里也能夠看出來R語言安裝包的三種途徑：CRAN的install.packages、bioconductor的biocLite 以及devtools::install_githublibrary(devtools) #包安裝器
library(RSkittleBrewer) #設(shè)置顏色
library(genefilter) #計算平均數(shù)和方差
library(dplyr) #排序整理結(jié)果7.2 第二步加載表型數(shù)據(jù)
pheno_data = read.csv("geuvadis_phenodata.csv") 這里使用測試文件自帶的csv文件。當(dāng)然實際運行時要自己用excel創(chuàng)建，其中包含了關(guān)于測序樣本的信息 csv文件(comma-separated values)就是逗號隔開的；表型數(shù)據(jù)7.3 第三步加載表達(dá)量數(shù)據(jù)
ballgown可以識別Stringtie、Cufflinks、RSEM的數(shù)據(jù)文件，讀取時需要設(shè)置三個地方
數(shù)據(jù)存放文件夾dataDir、樣本識別號samplePAttern，表型數(shù)據(jù)pData
bg_chrX = ballgown(dataDir = "ballgown", samplePattern = "ERR", pData=pheno_data)
將數(shù)據(jù)打包存放到bg_chrX這個對象中，它的核心就是一個矩陣，每一行是基因，每一列是樣本7.4 第四步過濾低表達(dá)基因
根據(jù)方差，過濾掉低于1的
bg_chrX_filt = subset(bg_chrX, "rowVars(texpr(bg_chrX)) > 1", genomesubset=TRUE)texpr是提取bg_chrX的表達(dá)量(reads_count)，rowVars對行(基因)求方差，意思就是：過濾掉只在一個樣本中表達(dá)，其他的樣本中都不表達(dá)的基因
7.5 第五步確定組間有差異的轉(zhuǎn)錄本/基因
例如，表型數(shù)據(jù)有兩個變量：性別與種群。要找不同性別(也就是這里的組)之間表達(dá)量差異的轉(zhuǎn)錄本，需要用種群的數(shù)據(jù)進行矯正。
7.5.1 使用stattest進行數(shù)據(jù)檢驗
#官方解釋：Test each transcript, gene, exon, or intron in a ballgown object for differential expression, using comparisons of linear modelsresults_transcripts = stattest(bg_chrX_filt, feature="transcript",covariate="sex",adjustvars = c("population"), getFC=TRUE, meas="FPKM") #featureL：在哪個水平進行計算，可選"gene", "transcript", "exon", or "intron" #covariate協(xié)變量：對因變量有影響的變量。它不是研究者研究的自變量,不為實驗者所操縱,但仍影響實驗結(jié)果。設(shè)置性別為協(xié)變量，因為本實驗主要是分析不同國家人群之間的差異。 #adjustvars：主變量，即不同的國家人群--> 字符串表示 #getFC：得到不同性別組之間利用無關(guān)變量矯正后的差異倍數(shù)(FC=Exp male/female) --> 向量表示 #meas：measurement采用哪種計算方式當(dāng)然也可以確定組間有差異的基因
results_genes = stattest(bg_chrX_filt, feature="gene",covariate="sex",adjustvars = c("population"), getFC=TRUE, meas="FPKM")7.5.2 在差異轉(zhuǎn)錄本中添加基因名稱與ID號：
results_transcripts = data.frame(geneNames=ballgown::geneNames(bg_chrX_filt), geneIDs=ballgown::geneIDs(bg_chrX_filt), results_transcripts)7.5.3 將轉(zhuǎn)錄本/基因按P值從小到大排序：
results_transcripts = arrange(results_transcripts, pval)
results_genes = arrange(results_genes, pval)
轉(zhuǎn)錄本的結(jié)果如下：排序后的轉(zhuǎn)錄本7.5.4 將轉(zhuǎn)錄本/基因結(jié)果寫入csv文件：
write.csv(results_transcripts, "chrX_transcripts_results.csv", row.names=FALSE)
write.csv(results_genes, "chrX_genes_results.csv", row.names = FALSE)7.5.5 找q value小于0.05的基因和轉(zhuǎn)錄本：
subset(results_transcripts, results_transcripts$qval < 0.05)
subset(results_genes, results_genes$qval < 0.05)差異表達(dá)轉(zhuǎn)錄本的結(jié)果如下：【匹配到了兩個已知基因的isoform】
表達(dá)差異轉(zhuǎn)錄本差異表達(dá)基因的結(jié)果如下：
表達(dá)差異基因
8. 數(shù)據(jù)可視化
8.1 【可選】圖層設(shè)置—使圖片更美觀
tropical= c('darkorange', 'dodgerblue', 'hotpink', 'limegreen', 'yellow')
palette(tropical)
8.2 展示樣本間基因豐度(FPKM值)的分布，根據(jù)顏色區(qū)分性別
fpkm = texpr(bg_chrX,meas="FPKM") #這里提取的是轉(zhuǎn)錄本的FPKM，當(dāng)然也可以提取其中的exon、gene等
fpkm = log2(fpkm+1)
boxplot(fpkm,col=as.numeric(pheno_data$sex),las=2,ylab='log2(FPKM+1)')
#取log的目的是為了讓數(shù)據(jù)能存在正值、負(fù)值，否則數(shù)據(jù)都為正不容易看出差別。這樣低表達(dá)小于0，高表達(dá)大于0，不表達(dá)等于0
8.3 展示單個轉(zhuǎn)錄本信息
比如要展示第12條：
ballgown::transcriptNames(bg_chrX)[12] --> 得到轉(zhuǎn)錄本名稱
ballgown::geneNames(bg_chrX)[12] --> 得到包含此轉(zhuǎn)錄本的基因名稱
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é)：最好的有參轉(zhuǎn)錄組分析軟件組合應(yīng)該是Hisat2+Stringtie+Deseq2
關(guān)于Deseq2以后會有介紹
