RNA-SeQC :annotation quality control tool

之前有討論過在拿到定序資料時,一開始的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

當我們將定序出來的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的程度

R Shiny 快速使用R搭建app應用

Rstudio 開發的shiny基本架構

Shiny application是一款以R語言為框架,可以快速且方便搭健application,當幫忙使用R寫完清理或是視覺化資料後,最重要的是如何讓不懂R語言的人也可以享受,此時,用RStudio開發的Shiny框架便能快速有這效果,且其跟Rstudio整合得非常好,貼心宜人,即使非寫網頁應用程式出生的生物學家,只要略懂R語言,便能完成!

screenshot.png

shiny語法有非常詳細的解說,且搭配實例,只要到官方網頁就可以瀏覽,大概一兩個小時就可以了解他的架構,最棒的事其可以跟html / javascripts等等串接在一起,做出更多的效果!

做好的application可以直接使用Rstudio裡頭的快捷鍵發佈到shinyapp的網路平台,其提供後台完整的服務,付費的話可以獲得更佳的效能,但免費版已經可以玩得很開心了!

screenshot.png

 

點擊寫好的code,便能直接deploy到shinyapps的平台,一條龍服務!

screenshot.png

基本上,shiny的概念,就是直接用R code搭建網站的前端(ui.R)和資料處理的後端(server.R),server的code會處理由ui端或是使用者輸入的資料,處理後再丟回給ui!

Kyruus 醫師資訊的彭博社

kyruus 公司網站

這是由本來為內科醫師之後轉往創投發展的Graham Gardner和亞裔工程背景的Julie Yoo在2010年創辦,其目標在於整合臨床醫師的資訊,提供給醫院、保險公司和藥廠,目前有搭建一個平台ProviderMatch,想要解決病人和醫療提供者之間的mismatch(因為在美國還要考慮到不同保險和藥廠的關係),所以建立這個平台來解決這個問題,將醫師是否availabe、位置、其隸屬醫院的保險公司等等全部整合在一起。

他會從多種來源來收集資料,像是人力資源資料庫、醫師認證系統、電子病歷、供應鍊資料庫等等超過一千個公共與營利性的資料庫。針對企業端,他們其中一個想要降低的是病人因轉診而流失在同一個醫療系統下的轉診流失率,這對很多醫療企業是非常在乎的。其公司在資料處理上主要有三大類:收集資料並整理、分析資料、

關於Tumor-Suppressor Genes

Tumor-suppressor genes其所轉譯出的蛋白質主要有三種功能存在:

  1. 其產生的蛋白質本身屬於細胞proliferation的機制中,但原本功能可能是減慢這個proliferation的速度
  2. 其本身為cell-cycle的檢查點調控的一部份
  3. 其本身跟細胞DNA損傷的修復或是細胞正常凋亡有關

那要如何研究此類的基因呢?相對於oncogen通常是顯性遺傳,tumor-suppressor genes則為隱性的alleles。所以兩個對偶基因都要有缺損才會造成疾病,而傳統上找尋tumor-suppressor genes可藉由分析特定的家族遺傳癌症來研究。

Retinoblastoma的研究歷史其實就是一個很好的例子。通常這類病人可以在13q14發現缺失片段,RB則被用來代表這個tumor-suppressor基因,比較有趣的地方是在細胞層次來說,RB-是隱性的,但在organism level卻表現出顯性的遺傳模式,目前解釋是某種rare event造成原本已經為RB-/RB+的另一個原本沒有缺損的RB gene被knock out,造成loss 0f heterozygosity。這種Retinoblastoma通常可分為兩大類,一種是有家族病史,另一種沒有家族病史,兩種的分子機制可能不同。

其他tumor-suppressor genes

  • Rb
    • 其為cell-cycle中的E2F Transcription factors inhibitor
  • p53
    • 為G1-S的檢查點之重要角色
  • HNPCC(Hereditary nonpolyposis colon cancer)
    • 此類大腸癌主要為chr2上的MSH2或是chr3上的MLH1有缺損造成,這些基因都跟DNA mismatch後的修復有關

如何找尋oncogenic mutation

閱讀Genetics: From Genes to Genomes, 5th edition 5th Edition by Leland H. Hartwell (Author), Michael L. Goldberg (Author), Janice A. Fischer (Author), Leroy Hood (Author), Charles F. Aquadro (Author)

如同上一篇所言,細胞癌化其實有兩類型的突變:oncogenes和tumor-suppressor genes,剛好是從相反的角度影響細胞的功能,大部分的細胞都是累積不同類的突變,慢慢轉變成癌細胞的,這篇來討論一下,科學家是如何找尋oncogenic mutation的呢?是哪種突變造成proto-oncogenes轉變成oncogenes的呢?oncogenes的突變是如何影響細胞癌化的?

目前找尋到oncogenes的實驗設計主要有兩種:

第一種是利用tumor virus,利用retrovirus將特定的cDNA嵌入目標的基因體中,再來研究其對細胞的生理影響,比如將可能的proto-oncogenes基因序列送入宿主基因體中,或是將想要表達的proto-oncogen或oncogen序列放入病毒基因體中enhancer或promoter的下游區段,使其在宿主中表達,另一種方式則是將病毒這段enhancer或是promoter的區段嵌入宿主可能的proto-oncogen或是oncogene的序列上游,觀察這些變化對細胞的影響,藉此來辨認oncogenes的正確序列位置。

另一種方式則是使用transformation assay,將癌細胞的DNA分離出來加到正常細胞株中,使正常細胞吸收那些癌細胞的DNA序列,看哪些正常細胞會因此影響,在看其是被哪些癌細胞的基因變段所影響藉此來尋找oncogenes,Alu sequence就是因此被找到的!

以下為三種目前被發現的oncogens:

Ras:其點突變造成其一直停留在GTP-actived form

c-Abl:在chronic myelogenous leukemia中被發現染色體9和染色體22上的c-Abl和bcr兩個基因互相交換,因此產生fusion gene,此fusion gene轉譯出的蛋白質為一種protein tyrosin kinase,但其不受調控

Her2:其本為human epidermal growth factor receptor 2,突變後變成over-expression 的狀態,不需要有growth factor的存在就能被激活

兩類主要癌化突變:oncogenes和tumor-suppressor genes

閱讀Genetics: From Genes to Genomes, 5th edition 5th Edition by Leland H. Hartwell (Author), Michael L. Goldberg (Author), Janice A. Fischer (Author), Leroy Hood (Author), Charles F. Aquadro (Author)

突變是如何造成細胞癌化的(一)

oncogenes和tumor-suppressor genes為兩類最主要造成細胞癌化的兩種突變基因類型

簡單來說,癌症基因在這先定義為造成細胞不受控制演變成癌細胞之突變的alleles,而這類的可以分成兩類:oncogenes、tumor-suppressor genes。oncogenes為突變的alleles,其主要是“顯性”地刺激癌細胞生長,在diploid的狀態下,單一alleles突變其實就可以造成細胞展現cancer-related phenotypes,其為gain-of-function mutation.而tumor-suppressor genes則是“隱性地”會使細胞癌化,為loss-of-function mutation,以汽車的踏板來譬喻,oncogenes是汽車的加油踏板出現問題,另一方面,tumor-suppressor genes則是汽車的煞車踏板出現問題。

古老的程式語言awk

程式語言不斷地隨者電腦科學的發展在演進,awk可以說是最早的scripting language,衍生自貝爾實驗室,在1977年出現,其為Perl(1987)和Python(1989)的前身,有一陣不太有人使用,但隨者有結構的資料大增,這種古老語言的處理速度,使其再度被使用。

就實用性來說,可以用其來取代cut,其在將資料分成field有非常大的彈性。

基本awk架構

awk ` pattern {action } `

  • awk處理邏輯是以行為單位的
  • 當pattern存在於line中,執行action
  • 將每行中的空白格當作是分隔
  • pattern的辨識以最大化為主

特殊變數special variables

在pattern或是action指令中有變數功能,如前述每個space中間即可被當作是一個單位的資訊

  • $0 整行
  • $1  第一個field(行中由空白格分隔裡,第一組文字)
  • $2  第二個field
  • NF 第幾個field
  • NR 當下的line數

指令operators

  • +-*/ 加減乘除
  • >< ==, !=比較
  • % 取mode
  • ~,!~ 有無match到

有用的awk教學資源

使用samtools來看特定區域的alignments狀況

samtools的功能非常強大,也有很多用法這篇有概略介紹而這邊仔細講解一下比較細緻關於alignment檔案的extraction,就是檢視特定區域的alignment狀況。

看特定區域的alignment狀況

在position-sorted和indexed過的bam file,我們可以擷取特定區域的alignment狀況使用samtools view這指令。

首先,先將想要使用的bam file用samtools index過

$samtools index sample.bam

接下來可以看指定的區域

$samtools view sample.bam 12:7788426-7796061

是將這段區域儲存成另一個bam 檔

$samtools view sample.bam 12:7788426-7796061 > Nanog_region.bam

或是有使用BED file儲存多的region,也可直接用samtools來看這些區域

$samtools view -L cancer_exons.bed sample.bam

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較多!