RNAseq: Annotation Based Quality Control

當我們將定序出來的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較多!