分類: RNAseq

文獻閱讀：RNAseq在臨床應用的可能和挑戰(中)

發表留言

原始論文：Byron, S. A., Van Keuren-Jensen, K. R., Engelthaler, D. M., Carpten, J. D., & Craig, D. W. (2016). Translating RNA sequencing into clinical diagnostics: opportunities and challenges. Nature Reviews Genetics, 17(5), 257–271. doi:10.1038/nrg.2016.10

文獻閱讀：RNAseq在臨床應用的可能和挑戰(中)

承襲上篇閱讀的內容，這邊談到論文的第二部分，即如何應用RNA sequencing於感染性疾病:

 可將RNAsequencing用來作為RNA-based相關病原體的診斷使用

目前臨床上重要的RNA virus如愛滋病（HIV）、伊波拉病毒（Ebola）、登革熱（Dengue）、肝炎（Hepatitis） 、流行性感冒等在醫院裡面都是使用qRT-PCR assay作為檢驗工具，這些都有可能在未來可直接用RNAsequencing來取代，除此之外，也可以用來追蹤某個疫情的爆發，像2014年的伊波拉病毒大爆發，便是使用amplicon sequencing的技術來追蹤。

用來診斷傳統培養方式難以確診的細菌感染

一部份的感染疾病，其感染部位的不可性，再加上感染的細菌量少，用培養的方式通常要很久，且不一定正確，如臨床上腦膜炎的感染，確診的話會需要抽脊髓液且不一定能培養出東西，要是能從病人血液中檢驗到少許病原體的RNA便能解決，畢竟抽脊髓液本身在臨床處置上是很侵入式的，另一方面，如Mycobacterium tuberculosis的診斷用傳統方式會有問題，使用 RNAseq來提早確診外（這類細菌培養的時間都非常長，甚至幾週），也能從其transcriptome來判斷其抗藥性。

但目前來說，相對於傳統的 qt-PCR assay，使用RNAsequence在傳統的實驗中會有很大的分佈統計上的問題！

HTSeq:Counting Reads per Genes with DESeq

發表留言

在之前有討論過RNAseq的分析中，基因表現的分析有很多不同的策略，一種採用定序片段（reads）對應到基因區域的數量代表基因表現（counting reads per genes），另一種則是使用RPKM(Reads Per Kilobase per Million mapped reads) ，兩種方式都有相對應針對資料正規化（normalization）的想法，可以參考這篇。

HTSeq 是用來處理alignment好的bam檔，將其產生counting reads table。而之後便可以把這資料丟到DESeq做基因表現分析（gene expression analysis），DESeq是R語言的package，可以用來做基因表現分析，其中用的原理便是採用counting reads per genes來做分析。

首先，HTSeq是一款python的軟體，除了計算reads count per genes 也可以用來處理基本的bam file writing 或是genomic interval的處理，其必須搭配安裝python 的一些週邊package，如PyPI、pysam，在MacOS X 上則先安裝xcode、MacPort(套件管理程式，用來安裝SciPy的套件)，另外要注意其支援的為python2.7版本。

在使用htseq-count時最重要的是他用三種mode，關於他如何定義reads對到genes的關係，一般defaut為union。官網說明連結screenshot.png

 

使用counting reads功能的簡單語法如下：

screenshot.png

 

 

基本上，輸入的bam檔必須經過sorting，最好是by reads names的sort setting，可以節省記憶體的要求，另外，annotation上，使用ensembl的annotation其gene ID和transcripts ID有區別性，對於程式來說比較好。

輸出會呈現如下：

screenshot.png

RNA-SeQC :annotation quality control tool

1 則迴響

之前有討論過在拿到定序資料時，一開始的FASTQC處理，這步驟可以用來品管一些關於“序列品質“的問題，可能樣品放過久等等，而alignment後的bam檔，因為有了序列片段為應到reference上的資訊，此時就可以利用這樣的資訊來進一步看是否有一些實驗設計上所造成的bias，在此篇有稍微介紹一下annotation quality control的概念和工具，這邊介紹RNA-SeQC tools的使用。

RNA-SeQC 是由Broad institute的研究團隊所開發的，主要是David S Deluca 專寫，為一款Java程式。其輸出會有html的介面和分析圖表，會有total and duplicate read count、mapping rate、read pairing stat、exon/intron/intergenic rate、strand specificity quantification、normalization by sownsampling、provide RPKM、variation of coverage、3’/5′ coverage rate、gaps in coverage report、stratification of gene expression、average transcript coverage、alignment mismatch rate per base、fragment size distribution等等，提供非常全面的資訊。

 

以下為使用的程式碼範例：

screenshot.png

 

使用這程式的前處理主要有fastq檔比需要先使用Picard 來產生sequence dictionary在fasta同一個資料夾中，bam黨必須要先使用samtools index過後，且bam檔的header必須要有@RG的註明，這部分也可以使用Picard的AddorreplaceReadGroup.jar來解決，另外事先bam檔案需要sorting過，根據chr1~chr22、chrX、 chrY、 chrM的順序處理，這部分可使用Picard ReorderSam.jar其可以搭配samtools將特定的reads除掉，像是chromosome patch等等。

螢幕快照 2016-02-23 上午12.49.54螢幕快照 2016-02-23 上午12.50.03螢幕快照 2016-02-23 上午12.50.12螢幕快照 2016-02-23 上午12.50.20

 

RNAseq: Annotation Based Quality Control

2 則迴響

當我們將定序出來的reads使用參考基因序列對應回位置後，我們便可以進一步計算每個基因上、exon、transcripts上對應的reads數量，來定量基因的表現，除此之外我們還可以用這些資訊來做進一步的實驗或是定序品質檢驗。

 

理論上來說，計算每個基因上所對應的reads數，可以很簡單的定量出基因表現，但在真核生物中，同樣的基因會有不同的isoforms，這樣便會出現一些問題，除此之外，在製備定序使用的library其protocol也會有影響，如rRNA有無去除乾淨等等，這些可以在alignment後對bam檔做annotation-based的quality control

 

雖然我們已經對fastq檔案做過初步的fastqc，但許多重要的資訊單純的fastq檔案無法提供，因為在fastq檔案裡我們只能對每個base-pair其定序出來的品質來檢查，如有沒有被rRNA汙染、或是製備protocol的檔案有沒有造成一些bias等等。

 

以下列出一些重要的品質指標(annotation-based quality metrics)

 

定序深度是否過深Saturation of sequencing depth

定序的深度其實會決定整個RNAseq實驗後續的基因表現分析、RNA剪接或是尋找特殊transcripts，所以到底實驗設計想要回答的問題需要定序多少深度會適合，而這決定於目標轉錄體的複雜性。

定序片段的在特定功能基因區分布是否有問題Read distribution between different genomic features

 可以看定序片段（reads）在不同的genomic features上的分佈，像是外顯子區域（exonic region）、內含子區域（intronic regions）、intergenic regions（基因間隔區），看分佈是否有特殊的地方，可能會發現新穎的基因或是isoforms，當然也可能會發現有無被污染，尤其是ribosomal RNA的比例，這可以知道是否在製備library時有清除乾淨。

transcripts上是否有均一地被定序到Coverage uniformity along transcripts

  不同實驗的步驟可能會造成transcripts上某一部份location被定序到的比例不同，比如說要是在製備library的時候，有poly-A capture的步驟，則其3’ ends就會有比較高被定序到的現象．

 

可以用來做annotation-based quality control的工具

 

  • RNA-SeQC
    • 為Java撰寫的程式，採用GTF格式的annotation來分析，另外，其需要bam所使用的reference fasta且要有.fai index和sequence dictionary file(,dict)，其會輸出html格式的檔案，裡面包含mean coverage、針對不同location transcripts coverage的計算（3’或是5’ bias)，會將基因的表現和GC content分成三個程度：low、medium、high。這篇文章有進一步介紹其功能
  • Qualimap
    • 為Java撰寫的程式，裡頭會使用到R/bioconductor的package，除了使用GTF格式，也採用BED格式作為annotation，其對於saturation和biotype 分佈的分析還不錯。其對於測序深度的分析會包含不同genomic features的深度分析、且會提供每增加一百萬的定序深度會提高發現多少新的features的機會，在biotype分佈分析，會將定序片段（reads）分佈在蛋白質轉譯區域(protein coding region)、偽基因（pseudogens）、rRNA、miRNA等地方顯示出來。
  • RseQC
    • 為python撰寫的程式，其特色是在計算不同genomic features上定序片段分佈時，會把上游和下游的transcripts都一起展現出來，另外一點特色則是其對splice junction的分析會標示出是否為novel
  • Picard’s CollectRNASeqMetrics
    • Picard工具裡其中一個Java 程式，其特色為會分析出strand specific的程度

Cufflinks: Transcripts Assembly Tools

發表留言

reference:
Trapnell, C., Williams, B. a, Pertea, G., Mortazavi, A., Kwan, G., van Baren, M. J., … Pachter, L. (2010).,Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nature Biotechnology, 28(5), 511–515. http://doi.org/10.1038/nbt.1621

cufflinks是一款用c++撰寫的程式，可以用來處理RNAseq已經align玩的資料，將其alignmenta組合成transcipts，並且找尋不同isoforms的可能，操作上在command line下使用。可以在這邊的網站下載以及使用說明書。

原理

第一步：在paired-end的設計下，其會把一對的reads一起看，將其排列成overlapped的bundles，程式會先找出incompatible的fragment，因為這些特別的點代表者一種spliced isoforms的可能，如下圖所示，黃色藍色和紅色分別都是被標示的incompatible fragments，整個資訊會用node和edge的方式串接起來

screenshot.png

第二步，會使用這種graph產生出isoforms的可能組合

screenshot

軟體使用的注意事項：

  1. 假如是paired-end reads實驗設計，要把reads identified的數值要改成一樣的。
  2. samtools必須安裝好，且路徑設計在PATH變數下

語法：
screenshot

連結

–overlap–radius:default 50 base-pair
 這是用來找出比較長距離的gene models

-g  reference.gtf
可以給定已知的gene model reference

 

輸出：

總共會輸出四個檔案：transcripts.gtf   isoforms.fpkm_tracking genes.fpkm_tracking

小建議：

假如要有多組sample要做cufflinks的話，建議先把bam檔合成單一個檔案，因為當alignment數越多，使用其演算法能找出的獨特isoforms較多！

 

samtools 使用

1 則迴響

samtools是在inspect或是處理bam/sam黨非常重要的工具，這邊簡單介紹他的功能和基本語法使用。

基本上，samtools對於sam/bam檔案有三大功能：

    第一 ，就是檔案格式的轉換（bam to sam/sam to bam）

    第二，sam/bam reads的篩選，過濾、查看

    第三，進階的variant calling分析（mpileup）

screenshot.png一 、檔案格式的轉換

我們可以使用samtools將bam檔儲存成sam黨，反之亦然，雖然tools一開始的defaut是input bam，但都是可以調整的。

bam 轉 sam:   $samtools view input.bam >output.sam

sam 轉 bam:   $samtools view -b input.sam > output.bam

二、篩選，過濾、查看sam/bam

排序bam檔案： ＄samtools sort -f input.bam output.bam

索引bam檔案： ＄samtools index input.bam (會跑出 *.bai 檔案 )

通常再匯入IGV或是某些軟體前會需要進行這兩個步驟

過濾特定flag的reads

留下特定flag的reads $samtools view -f [FLAG value]

濾掉特定flag的reads $samtools view -F [FLAG value]
(關於samformat的flag可參考此)

可以看alignment的depth和reads mapping status，這功能蠻好用的！

看每個reference序列定序深度
＄ samtools depth data.bam
綠色框框代表每個序列對應到的reads疊加數

螢幕快照 2016-01-27 下午11.46.06

＄samtools idxstates data.bam
黃色代表是左邊染色體上有被mapped上的reads數量，右邊代表的是unmapped量

螢幕快照 2016-01-27 下午11.46.12

第三種功能則是calling variant

螢幕快照 2016-01-27 下午11.51.16

這邊就需要調整比較多參數，其實這些功能蠻方便的，可以快速地瞭解整體bam檔的狀態！

關於sam/bam 檔案格式可以參考這篇

SAM, BAM and CRAM

1 則迴響

常見的Alignment Genomic Data Archive Format

alignment後的基因定序資料通常是任何分析的起頭，如廚師把配菜弄好擺在流理台後，要大顯身手，此時一個好的工具來處理就是非常重要，不然會有種看者滿山資料卻不知從何下手。alignment後的定序資料通常就是已有配對到reference genome上的reads，目前此類檔案有三種format: SAM, BAM, CRAM，檔案大小依序也是從左到右，而SAM(Sequence Alignment/Map)和BAM(Binary Alignment/Map)兩個檔案基本上是同樣的內容，差別在SAM是人看得懂的，BAM是壓縮成binary的，前兩種檔案格式出現的時候，定序資料量還沒現在這麼大，最近連BAM格式還是被嫌太大，CRAM format是會逐漸取代前者的壓縮方式，EMBL已經改用CRAM的格式來儲存定序檔案了，其概念是所謂的reference 壓縮的模式，如下圖所事是：screenshot.png

簡單來說，就是以reference上相對位置和改變的資訊來儲存，而非直接儲存reads序列，
整個降低檔案儲存需要的大小。使用cram格式相對於bam可以節省約10-30%的檔案大小，但因為cram的格式是基於reference來減低檔案格式，需要固定的reference版本，才不會再還原時候有問題，所以每個CRAM檔案會有MD5 checksum，來確保記錄所使用的reference版本。

關於CRAM的相關reference:

Hsi-Yang Fritz, et al. (2011). Genome Res. 21:734-740

Cochrane G. et al. (2012). GigaScience 1:2

SAM Format的形式

SAM的檔案架構分兩部分，一部分是Header，另一部分是Alignment。外觀如下圖：

screenshot

Header的部分位在整個檔案的一開頭，裡面會有這份檔案的基本資料，如有無sorted過、裡頭使用的reference及其長度、使用過什麼工具處理和alignment時所下的指令。Header裡頭會有幾個關鍵組成：＠SQ/ @RG/ @PG，分別代表以下的意義

＠SQ 這個開頭含有做alignment所使用的reference序列資訊，SN代表的是sequence name，而緊跟者為LN，代表此參考序列的長度

@RG代表者read group的資訊和sample的基本資料，部分軟體會根據裡頭的ID，去辨識有無batch effect，有的會有PL指標去代表其所使用的定序平台

@PG裡頭含有此次alignment所使用的程式資訊，CL指標後面有所下的指令，VL後面有使用的軟體版本資訊

Alignment的部分就是實際一個個read的alignment資料，一行為一個reads，其用以下的資訊來說明每個reads的狀況：
screenshot.png

依序是每個Reads獨一無二的編號，通常還有在定序機器上的位置，FLAG，則是代表此reads的mapping狀態，chr則很明顯的是對到的染色體，Start代表reads對應到reference上的開始位置，CIGAR則代表這個reads每一個位置跟reference mapping的狀態，是否100Match到，後面還有一些進階資訊，假如是paired-read則會顯示他的mate reads的位置。

因為sam檔沒有壓縮過，所以大多數會以bam檔形式儲存，必須使用samtools來轉換，以便閱讀上面那些資訊，這邊samtools在處理bam/sam檔案上就非常重要了！

可參考這篇介紹samtools的使用！

基因的長度是否會影響RNAseq中的表現差異分析呢？

1 則迴響

論文閱讀：Detecting differentially expressed genes by smoothing effect of gene length on variance estimation

關於RNAseq Normaization更進一步閱讀

表現差異分析(Differentially Expression Analysis)是RNAseq實驗很重要的目的之一，因為我們想了解各種狀態間基因表現的差異！但目前大部分用來做DE analysis的演算法有一個問題，便是傾向於將較長的基因辨識為有表現差異的基因，而比較不容易將長度短的基因辨識為有表現差異的。

這篇研究認為這演算法的偏差，是來自於其計算過程中沒有將gene length對於reads count的variance estimation考慮進去，所以想提出一種改進過的方法LenSeq，將基因的長度smoothing，以達到讓有表現的短片段基因能被辨識出來！

原本科學家都認為RNAseq可以unbiased的用來比較基因的表現量，但實際上定序出來的reads是正比於此transcripts的長度，且定序深度也跟reads的總處相關，

換句話說，在資料前處理的normalization需要考慮進去這個bias:reads count 正比transcripts length，有人會借用microarrary data所常使用的quantile normalization，但其實這個問題在microarrary的設計中不會遇到，所以這標準化的分析並沒有考慮進去這偏差。

目前第一步用來表準化transcripts表現量的數值為RPKM(reads per kilo base per million mapped reads)=(read counts  X 10^9)/( gene length X library size)或是後來的RPKM。可以用一個簡單的例子來思考：
Reads count for gene
Total Library Size
Library 1
800, 800, 800, 800, 800, 800, 800, 800, 800, 7800
15000
Library 2
1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 1000, 6000
15000

當這兩個library size一樣的時候，會發現上面十個基因都是DE genes, 但實際上是不合理的，所以單純的使用library size表準化是不夠的。

後來的DESeq和edgeR，則有更進一步的方式來進行標準化，他們的假設建立在大多數的基因表現是沒有Deffierential Expression的，edgeR有工具可以trimmed mean of M-value，其會將read counts較高的數值去掉，剩下的再來計算兩個要比較的library之weighted mean of log ratio ，而DESeq則是使用scaling factor來標準化，先計算每個基因在不同lanes間的geometric means，取得每個reads count對應到其geometric means的ratio，在使用這群數值的mean 來代表scaling factor，最後在應用這數值來標準化兩個library。本質上這樣的演算法是將read counts過大或是過小的transcripts濾掉，從新調整資料的centrality。

小結論，仔細思考RNAseq的normalization步驟是很重要的，其影響到downstream analysis。

Bullard JH, Purdom E, Hansen KD, Dudoit S, Evaluation of statistical methods for normalization and differential expression in MRNA-Seq experiments, BMC Bioinfor- matics 11(1):94, 2010（被引用673次）

留一些關鍵字待查：

RNAseq followed Poisson distribution

Robinson MD, Smyth GK, Moderated statistical tests for assessing diFFerences in tag abundance, Bioinformatics 23(21):2881–2887, 2007.（被引用356次）

Negative binomial distribution

Whitaker L, On the Poisson law of small numbers, Biometrika 10(4):36–71, 1914（被引用65次）

RNAseq: Normalization

1 則迴響

Han Y, Gao S, Muegge K, Zhang W, Zhou B. (2015) Advanced Applications of RNA Sequencing and Challenges. Bioinform Biol Insights 9(Suppl 1):29-46. [article]

計算好基因的表達量，之後的data normalization其實是很重要的一環，這一步驟要考慮的細節很多，像是transcripts size、GC content、sequencing depth、sequencing error rate、insert size。

在評估各種標準化的方式，可以使用measurement error model來評估其標準化的效力。

目前有很多標準化的方式，像是quantile normalization，這方法在處理microarray資料上很有效，在RNAseq中能提高data quality，即使是low amount的RNA。

在R語言中的EDASeq，其先使用within-lane normalization再來between-lane normalization，可以有效降低GC-content所造成的問題

Lowess normalization的方法在microRNA 的資料處理蠻有效的。

目前在處理RNAseq資料上的normalization依舊需要有更好的方法被發展出來。

進一步閱讀：基因的長度是否會影響RNAseq中的表現差異分析呢？