文章題目：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測序能夠在一次試驗中捕獲成千上萬基因的表達量症副，能產(chǎn)生原始reads數(shù)十M发皿，通過檢測每個基因測出的reads數(shù)目可以估算基因豐度粘姜，當(dāng)有兩個或多個重復(fù)時能夠判斷差異表達基因形入。另外利用RNA-seq可以檢測注釋文件中沒有的變異基因，也可以尋找每個基因不同的調(diào)控表達途徑藕赞。
分析RNA-seq一般需要四步:
- 將reads比對到參考基因組上（這篇文章是以有參轉(zhuǎn)錄組為例）成肘；
- 將比對上的alignments組裝成全長轉(zhuǎn)錄本；
- 基因與轉(zhuǎ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負責(zé)拼接轉(zhuǎn)錄本囤萤，構(gòu)建isoforms用于估算基因表達量；【先接收Histat傳來的alignments進行轉(zhuǎn)錄本拼接是趴，每組數(shù)據(jù)之間是獨立的涛舍，拼接的過程中估算每個gene和每個isoform的表達量；拼接完以后唆途，所有的拼接好的序列進行merge富雅，這一步是必須的，因為某些樣本中的轉(zhuǎn)錄本只是被reads部分覆蓋肛搬，導(dǎo)致的結(jié)果就是僅有那些被覆蓋的區(qū)域得到了拼接没佑，利用merge可以降低這一誤差∥屡猓】merge完成的轉(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é)果章郁，計算得出差異表達轉(zhuǎn)錄本【接受StringTie傳來的轉(zhuǎn)錄本和定量數(shù)據(jù)，根據(jù)設(shè)置的不同的實驗條件進行分組】它包含了某些R包可以直接可視化
這三件套流程支持帶有時間線的實驗(比如研究某些病的發(fā)育歷期)以及兩個以上實驗條件下的結(jié)果比較
流程很容易操作痢掠，但需要會使用命令行及R腳本運行
當(dāng)然驱犹，這三個工具并不是缺一不可：如果你是要做無參轉(zhuǎn)錄組，可以用其他工具如trinity拼接好以后傳給StringTie足画；也有一些項目只需要驗證已知的某些基因在樣本中的表達量雄驹；另外即使對于模式物種，基因注釋文件也不完善淹辞，因此可以結(jié)合許多其他工具使用医舆。雖然這三種工具能夠檢測差異基因表達，但對于外顯子的差異不能探索象缀，好在可以結(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ū)段，速度高達每小時比對幾百萬條reads。以上方法在DNA比對中最為常用窟她，例如組裝基因組需要比對拼接陈症。但是對于RNA來講，有一個DNA中不存在的問題需要注意：成熟的mRNA將introns切除震糖，因此許多RNA測序的reads要跨越內(nèi)含子進行比對录肯，因此一條短序列可能會被拆開比對到相隔1萬bp甚至更遠的位置 【人類平均內(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é)果準確性蜡镶，可以更快找到新的剪切位點、轉(zhuǎn)錄起始和終止位點 -
StringTie拼接恤筛、融合
首先將比對完的reads劃分成不同的基因座組官还，然后把每個基因座拼接成isoform。運用了network flow algorithm 毒坛，從reads數(shù)最集中(也就是表達量最高)的轉(zhuǎn)錄本開始望伦，然后移除已拼接完的reads林说，再次重復(fù)這個過程，最后拼接成眾多的isoform屯伞，結(jié)束的標志是：全部的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ù)已注釋好的直接命名。注釋文件對低表達量的基因很有幫助勾习，因為reads數(shù)量太少浓瞪，不好拼接得到該基因。拼接完后语卤，可以用
gffcompare檢測有多少拼接的轉(zhuǎn)錄本或全部或部分地map到了已知的基因追逮，并且統(tǒng)計有多少是全新的基因〈舛妫【此環(huán)節(jié)可以跳過】想要比較不同樣本基因表達量差異钮孵，如果通過一個一個基因?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表達量。
StringTie可以不依賴于注釋文件進行新的基因贼急、轉(zhuǎn)錄本預(yù)測茅茂，注釋文件只是個佐證
-
Ballgown差異表達分析
這算是上游Hisat捏萍、StringTie與下游R/bioconductor的中間橋梁，他將數(shù)據(jù)進行標準化
基本上差異表達分析都會包括下面幾個步驟：(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ù)可視化和檢查
- 差異表達的統(tǒng)計分析
- 多重檢驗校正
- 下游檢查和數(shù)據(jù)summary
在R里面使用Ballgown處理需要得到：
- 表型數(shù)據(jù)：關(guān)于樣本的信息
- 表達數(shù)據(jù)：標準化和未標準化的外顯子玉吁，轉(zhuǎn)錄本照弥，基因的表達數(shù)量
- 基因信息：有關(guān)外顯子腻异，轉(zhuǎn)錄本进副，基因的坐標以及注釋信息
使用Ballgown：
- 讀入數(shù)據(jù)：將上游Stringtie輸出的轉(zhuǎn)錄本表達量數(shù)據(jù)與表型數(shù)據(jù)結(jié)合起來【注意保證二者的ID號一致】
- 豐度預(yù)測：以FPKM（fragments per kilobase of transcript per million mapped reads）為單位的豐度預(yù)測根據(jù)library size進行標準化。
- 差異表達分析：FPKM的數(shù)據(jù)解讀需要取log轉(zhuǎn)化一下悔常，再建立線性模型
- 可以在基因影斑，轉(zhuǎn)錄本，外顯子水平上進行差異分析
- 結(jié)果以表格形式展現(xiàn)机打，樣本量大的話還有p值和q值
第三部分：數(shù)據(jù)測試
軟件準備：
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ù)準備：
下載RNA-seq測試數(shù)據(jù)
其中包括人類基因組X染色體RNA-seq，基因注釋文件 残邀〗粤桑【X染色體參考基因組、測試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(此處默認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)錄本合并：
標準格式：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標準使用：
gffcompare –r chrX.gtf –G transcripts.gtf
-r：參考注釋文件
-G：拼接的轉(zhuǎn)錄本，也是要被比較的對象
-o：指定輸出的結(jié)果前綴碗旅，否則默認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)錄本表達量并生成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差異表達分析
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 第三步加載表達量數(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 第四步過濾低表達基因
根據(jù)方差蹋凝，過濾掉低于1的
bg_chrX_filt = subset(bg_chrX, "rowVars(texpr(bg_chrX)) > 1", genomesubset=TRUE)texpr是提取bg_chrX的表達量(reads_count)鲁纠，rowVars對行(基因)求方差，意思就是：過濾掉只在一個樣本中表達鳍寂，其他的樣本中都不表達的基因
7.5 第五步確定組間有差異的轉(zhuǎn)錄本/基因
例如，表型數(shù)據(jù)有兩個變量：性別與種群情龄。要找不同性別(也就是這里的組)之間表達量差異的轉(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)差異表達轉(zhuǎn)錄本的結(jié)果如下：【匹配到了兩個已知基因的isoform】
表達差異轉(zhuǎn)錄本差異表達基因的結(jié)果如下：
表達差異基因
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ù)能存在正值佑笋、負值翼闹，否則數(shù)據(jù)都為正不容易看出差別。這樣低表達小于0蒋纬，高表達大于0猎荠，不表達等于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以后會有介紹
