分類: 未分類

Principle Component Anayasis練習：口腔癌、胃癌和大腸直腸癌基因表現資料（一）

1 則迴響

notes from the GSB Bioinformatics study group

Microarray資料來源

相對於Next Generation Sequence 資料動輒幾十gb，這次我們選擇用microarray的資料練習，
首先我們分別從GEO、EMBL和TCGA裡找尋用affymetrix 測expression的實驗資料，畢竟affymetrix是microarray年代裡最可靠的平台，另外我們也希望使用的chip 是一樣的，避免到probe annotation的問題，最後我們找到三組可行的資料來練習分析，都為使用Affymetrix Human Exon 1.0 ST array(GPL5175)：

第一組：口腔癌 GSE25099

這組資料來自於台灣長庚醫院頭頸癌團隊的研究論文，此組包含有79個樣本，其中22個為正常組織

第二組:胃癌 GSE33335 GSE13195

這組為來自中國大陸的研究團隊，其基因表現資料來自於27個paired胃癌組織樣本

第三組：大腸直腸癌 GSE29638

這組來自於挪威的研究團隊，主要有207組來自於挪威的大腸直腸癌病人檢體，其臨床分期有第一期到第五期，實驗主要是用來找尋第二期大腸直腸癌是否有可以用來預測預後的基因表現biomarker

匯入R並且簡易分析

參考此篇來做匯入

先安裝GEOquery R bioconductor package

source(“http://bioconductor.org/biocLite.R")
biocLite(“GEOquery“)

library(GEOquery)

接下來使用GEOquery來下載我們找到的檔案，先嘗試GSE25099

gse <- getGEO(“GSE25099“, GSEMatrix = TRUE)

screenshot.png

接下來我們只要這資料裡頭的gene expression matrix，於是我們把exprs subsetting出來，因為這是s4 class，先用＠在用$，交替尋出我們要的matrix

screenshot.png

接下來快速使用R內建函數來做簡單分析，可參考這篇，使用prcomp()

screenshot.pngscreenshot.png

RNAseq:Differential Gene Expression

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]

轉錄體（Transcriptome）代表著單一細胞或是一群細胞中所有的RNA，而RNAseq其中一個很重要的應用，便是藉由比較不同情況下，轉錄體的表現差異，來探討細胞在不同發育階段、特定刺激下的基因表現等等。比較基因表現量差異的關鍵便是要能辨別不同isoform transcripts的量，另一方面，也代表要對genomic functional element的劃分要更清楚，因為這些都會決定某一個transcripts的表現量高低，這種比較不同狀況下轉錄體表現量差異的分析就叫做DEG analysis，目前已有些統計的方式來計算這些轉錄體中不同transcripts的表現量，有部分方法是來自於之前microarrary資料處理的方法。

最常用來表達某一個transcripts的表現量為RPKM(reads per kilo base per million mapped reads)，這數值顧名思義就是將transcripts的長度以及總量對此transcripts的表現量做normalization後的數值，但其實還有許多因素會影響到評估transcripts的表現量，像是sequencing depth、gene length、isoform abundance。目前RPKM這統計數值被質疑的地方便是其每有考慮到同一個transcripts isoform的影響，而另一種新的數值RSEM(reads per expectation maximization)，其有考慮進去isoform的角色，

大多數用來分析的演算法，其只用最簡單的count-based 機率分布（Poisson distribution）以及Fisher’s exact test來做表現量差異的統計檢定。

這樣的檢定方法，並沒有考慮到biological variability，而要降低biological variability可以藉由分析重複實驗的資料，使用permutational-derived methods。

另一方面，假如資料有足夠多實驗重複資料時，可以使用extended Poisson distribution來做更進階的分析。

但RNAseq的研究非常的貴，要有多組重複實驗顯得不切實際，也有人在有限的檢體實驗中，直接模擬biological variability，並且使用pariwise或是multiple group comparisons來分析。有幾個工具已被開發出來進行表現量比較的分析，如Cuffdiff、DESeq、EdgeR。

在read counts計算上，其分布常呈現skewed，需要使用演算法來把分佈轉換，才能繼續使用統計檢定，有研究者開發出一款工具PoissonSeq，便是使用Poisson log-linear model來做DEG分析。

除此之外，有科學家使用從microarrary分析而來的工具來做RNAseq的表現量分析，其使用limma package中的voom 函式來將資料轉換成Gaussian distribution，在進行統計檢定！

這邊有一篇研究來比較目前所有的DEG分析，可以參考此論文的意見再來決定要用什麼分析工具

Soneson C, Delorenzi M. A comparison of methods for differential expression analysis of RNA-seq data. BMC Bioinformatics. 2013;14:91.(被引用231次)
91. Kvam VM, Liu P, Si Y. A comparison of statistical methods for detecting dif- ferentially expressed genes from RNA-seq data. Am J Bot. 2012;99(2):248–56(被引用125次)

關於DEG分析還有許多挑戰需要解決，像是

How to uniform the reads coverage along the genome with the nucleotide composition variation;

How to detect the “within-sample” variations without simply assuming that the underlying conditions or treatments affect all individual gene equally;

How to improve current methods to detect differences in gene isoform preferences and abundance level in varying conditions;

How to account for the dif- ferent probability in read coverage in long genes versus short genes since we can gain great sequencing depth nowadays.

如何使用SRA toolkits下載NIH的NGS data

發表留言

screenshot.png

SRA database (Sequence Read Archive)是目前次世代定序實驗資料原始檔案的儲存地方，在越來越追求研究reproducible的時代，大多數基因體研究都會要求將資料上傳到雲端，讓其他科學家也能使用， 所以看到跟自己相關且興趣的基因體相關論文，想要自己下載其實驗定序出來的資料來分析一下，便可以先來NIH的網站瀏覽一番，看有無自已有興趣的資料，在用其原始資料accession number來下載，但次世代定序的資料動輒好幾GB，明智的做法是使用公用的linux server去下載，此使可以利用NIH 開發的SRA toolkit，這軟體主要是在command line環境下操作，只要知道自己想要下載資料的accession number，便能用一行指令下載好幾gb的檔案，甚至可以進行簡單的檔案轉換處理，而只要使用SRA toolkit便能簡單做到，這邊介紹一下如何用SRA toolkit來下載心儀的資料。

在此之前，我們先來了解一下SRA database其對檔案命名的眉角，其實是有一定的邏輯的，大致如下的命名編號架構（這邊指的是accession number）：

 

screenshot.png

簡單來說，每個研究都會有個計畫編號，開頭大多會是PRJN或是SRP,這種計畫通常都會不只做一組檢體，而每個送去做定序的檢體都會有專屬編號，以SAMN開頭，而此檢體使用的定序實驗設計，也會有個編號，以SRX開頭，裡面會記錄此檢體所使用的定序機器種類、library type、定序設計（paired-end, single ….)，而在同一次定序中其跑的run也會有編號，以SRR開頭。

接下來寫一下，SRA toolkit的下載和使用：

第一步：到此頁面下載符合自己系統的版本 link

screenshot.png

第二步：在系統下解壓縮
這邊以linux/mac系統的指令為例
＄tar -xzf sratoolkit.current-centos_linux64.tar.gz

第三步：設定至系統路徑
此tool解壓縮完，其實就可以使用了，其tools的執行擋放在/sratoolkit/bin/的目錄下，直接打想要的tool linux是無法直接啟動，所以要去修改bash_profile，參考這篇網誌 如何設定檔案室系統路徑變數中

第四步：常用tools
其中最常用的工具就是fastq-dump了，可以用其直接使用accession number來query NIH 的SRA database並且下載fastq檔案，這邊是其指令介紹：screenshot.png
主要指令依功能分成三大類：

  1. data formating：下載後的檔案希望為什麼樣的格式，可以為fasta、fastq、sam，可以存成一個檔案，或是每個reads分開存 
  2. filtering：指定要下載特定區間內的reads
  3. workflow and piping：檔案最後要直接顯示在螢幕、或是儲存到特定的目錄中

    screenshot.png

 

 

 

 

 

Fusion genes是什麼？

發表留言

閱讀自Reproducible, Scalable Fusion Gene Detection from RNA-Seq.Arsenijevic V1, Davis-Dusenbery BN2.

screenshot.png

在1960年代，便在chronic myelogenous leukemia的病人身上發現BCR-ABL1 fusion genes的發生（第22對染色體 和 第 9對染色體），而最近因為sequencing的定序進步，幾乎所有分型的癌細胞上面都多少可以找到fusion genes的蹤跡，也懷疑這些fusiongene是否有driven 這些癌症的角色。

Fusion genes可由兩種方式來分類，第一種根據其嵌合基因的位置，兩個位於不同染色體的基因嵌合在一起為(interchromosomal)， BCR-ABL1 fusion gene的機制便是屬於這一類的，或是fusion genes來自同一個染色體上(interachromosomal)，像是造成攝護腺癌症的TMPRSS2-ERG fusion gene（發生在百分之五十的攝護腺癌），另一種功能上來看，則是嵌合的位置在regulatory region，影響其基因的表現量，不產生新蛋白質產物，或是發生在coding region，造成新的蛋白質產物發生。

但為何在癌細胞會產生fusion gene呢？是由什麼樣的機制造成fusion gene呢？目前一種觀察發現這些fusion gene所發生的chromosomal rearrangement的區域其有copy number transitions 和 高度表現量，是否是這樣的原因讓這樣的fusion能有selective advantage而繼續存在於細胞中呢？

另外有論文說可能是因在double strand 分開後DNA repair機制被改變，造成chromosomal translocation的機率上升！

RNAseq: Reads counting

發表留言

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]

在RNAseq資料裡面，常使用對應到某個基因序列的reads量來代表其基因表達量，當然要直接把reads數量對應成基因表現量，這邊一定要仔細考慮在定序前檢體處理的步驟都有相關，比如是否是strand-specific，這些都要列入考慮的因素！

這邊有一些工具可以計算reads counting的量，如Bedtools，給定其GFF和bam檔，其可以在特定的區段計算出counts，其內建是會計算兩個strands上有興趣區域的reads一起計算，但這能調整。

另一款工具HTseq則是有把strands-specific概念在預設計算中，其只會把reads對應到特定的strands上。


在R語言裡，easyRNAseq使用上簡易上手、summarizeOverlap為GenomicRanges裡頭的韓式、featureCount (Rsubread)裡頭則採用了非常有效率的chromosome hashing和feature blocking 的方式。

除了使用工具不同，計算出來的count值也會不同，所使用的gene model 同樣影響到後續的計算，尤其是在跨過junction的reads上，有論文發現在使用RefGene、Ensembl、UCSC三種版本對照序列，發現其中21958個常見基因只有16.3%的基因其算出的count相同，其中9.3%基因有大於50%的差異。

RNAseq: Reads mapping

發表留言

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]

在將RNAseq 定序出來的reads經過preprocessing後，接下來的步驟便是要將其組成的contigs對應到參考序列上，如何將短片段的序列對應到正確的位置一直以來都是生物資訊領域重要的問題！目前已經有很多不錯的軟體能處理這類問題，像是ELAND、SOAP、SOAP2、MAQ、Bowtie、BWA、ZOOM、STAR。這裡有一篇比較這些工具的論文可以參考看看。

Bao S, Jiang R, Kwan W, Wang B, Ma X, Song YQ. Evaluation of next-generation sequencing software in mapping and assembly. J Hum Genet. 2011;56(6):406–14

這些工具在處理一般序列區域的功能都相當不錯，但遇到像是poly(A) tail或是exon-intron splicing junction就必須要針對splicing處理的工具像是BLAT、TopHat、GEM、MapSplice。

在對應回參考序列另一個問題就是序列的polymorphism，這會造成multiple-aligned reads的問題，在較短的reads於序列重複性較高的區域這問題會更嚴重，而在一些重複性較小的reads上，可以利用附近的reads來解決這類型的問題。

Cloonan N, Forrest AR, Kolle G, et al. Stem cell transcriptome profiling via massive-scale mRNA sequencing. Nat Methods. 2008;5(7):613–9

Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat Methods. 2008;5(7):621–8

另一種解決方式就是將paired-end reads 延升成200-500bp的長度，在進行mapping。

Holt RA, Jones SJ. The new paradigm of flow cell sequencing. Genome Res. 2008; 18(6):839–46.
Hillier LW, Marth GT, Quinlan AR, et al. Whole-genome sequencing and variant discovery in C. elegans. Nat Methods. 2008;5(2):183–8.
Campbell PJ, Stephens PJ, Pleasance ED, et al. Identification of somatically acquired rearrangements in cancer using genome-wide massively parallel paired- end sequencing. Nat Genet. 2008;40(6):722–9

給生物學家的計算分析（六）：shell下建立scripts

1 則迴響

閱讀自Practical Computing for biologists by Haddock Dunn

完整建立一個scripts是幫助自己在command line完成許多事情的力氣，就不用再做很多重複的工作。

基本概念：所有在command line底下的指令本就就是一組scripts

所有在command line下的使用指令，本身就是一個scripts，只是他一開始就安裝好了，如ls、cd、rm..等等，我們來看看他們放哪邊！screenshot.png

我們使用which這個指令可以用來找尋物件在系統中的路徑，會發現這些基本指令都放在 /bin 裡面，這資料夾放滿系統各式重要的資料，當我們輸入cd按下時，系統會去尋找這幾個資料夾內的scripts，看有沒有符合的名稱，有的話就會執行！那哪邊可以知道放置這些系統會去尋找的特定資料夾呢？

screenshot.png

使用echo $PATH，就能叫出目前系統當指令輸入後會去尋找的路徑！各個路徑用:隔開。所以我們只要把建好的檔案放進去這些資料夾中，便能像是cd一般的指令來使用，但把自己弄的“土炮”scripts弄髒系統裡面核心功能資料夾，似乎不是個很好的方式，所以另一方式就是建立一個資料夾，然後把其路徑放入＄PATH中，似乎是個不錯的方式！那要如何做呢？

以MacOS系統為例，每個帳號都有一個系統隱藏檔案.bash_profile，是來儲存各種個人的設定，（在linux中為.bashrc），可以在母目錄裡輸入 ls -a 看到此資料夾screenshot.png

在裡頭輸入

screenshot.png

那如何將command lines結合成scipts可以參考這篇文章：如何將command lines 組合成shell script

癌症基因體閱讀心得（一）

發表留言

閱讀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)

 

簡單來說，癌細胞就是突破了一般細胞生長的限制，正常細胞會有以下的機制去制衡其生長：

  1. contact inhibition
    正常細胞彼此有接觸後，就會停止，在cell line中長成monolayer，而癌細胞則不會受限制，會形成所謂的transformed loci
  2. apoptosis
    大部分的細胞在受損後，假如無法修復其基因的改變，就會執行自我凋亡，但癌細胞不會，即使其基因序列已經突變，無法修復，其依舊繼續存活者
  3. autocrine stimulation
    正常細胞在沒有接受外界刺激時，其不會proliferation，但癌細胞會分泌自我刺激的激素來活化自己分化或是複製

細胞的癌化

細胞累積了一定量的突變才轉化為癌細胞，對照流行病學，年齡增長和癌症發作比率同為正相關，這樣的解釋蠻合理的，分子層面來看，這些突變可能是來自於基因中特定的strands，慢慢累積成不同的clones，而癌細胞有兩個特色：prolifereation rate上升加上癌細胞genomic instability，整個為一個惡性循環。

 

RNAseq’s 資料前處理:Quality Control

發表留言

RNAseq資料處理的步驟非常繁複，而第一步就是要看定序出來的品質如何，品質不好的話是否要去掉差的定序reads，特意設計疊加的reads是否要先處理完後再繼續下一步的aligment、是否有被污染、過多的adaptors(比較短的library中容易碰到的問題)。

最常見用來評估定序品質的工具為FastQC

FASTQC 是一款用來分析RNA-seq FASTQ 檔案的軟體,可以在command line環境下執行.

分析後會取得的資料:

  • 檔案名稱:最原始分析前之檔案名稱
  • 檔案類型
  • 解碼(encoding):其品質數值所表示的ASCII
  • 全部定序數量(Total Sequence):總共有多少定序片段
  • 過濾定序(filtered sequence):在Casava模式下,在分析過程中片段序列被標誌且移除
  • 定序長度(Sequence Length):最短和最長的定序序列
  • GC比例(GC content)

FastQC 執行完,會給出11張圖:

Per Base Sequence Quality 每個base平均品質
screenshot.png

    1. 分析內容
      1. 可以用來看是否有特別低品質的定序序列,通常都是因為照向品質不好所造成的,像是位在定序晶片邊緣的點
      2. 達不同品質指標,圖標會顯示之圖例
        • 警告(Warming):當平均品質小於27/0.2% error rate
        • 失敗(Failure):當平均品質小於20/1% error rate
      3. 備註:經過quality trimming後可能會有偏差,或是看不出error batch?

Per Sequence Quality Scores 每條序列中不同位點的品質
screenshot.png

    1. 分析內容
      1. 顯示序列裡不同位置平均品質範圍(只要是fastq中之資料)
      2. BoxWhisker type Plot圖中標示的解釋:
        • 中間紅色線:中間值(median)
        • 黃色長條:四分位距內的值(25-75%)
        • 上下whiskers: 各代表10%和90%的點
        • 藍色線:代表平均值(mean)
        • 背景色
          1. 綠色very good quality
          2. 橘色reasonable quality
          3. 紅色poor quality
        • 標題上是FastQC猜測此定序所使用的方式,在虛列品質好時比較容易猜測錯誤.
      3. 達不同品質指標,圖標會顯示之圖例
        • 警告(Warming):當25%之任何序列位點<10或是median值<25
        • 失敗(Failure):當25%之任何序列位點<5或是median值<20
      4. 問題導因:
        • 原因一: 跑定序的時間過長,造成後面定序出來的品質下降
          1. 處理方式: Quality Trimming(依據平均品質來決定裁掉的長度)
          2. 通常在adapter read-through也會有問題
        • 原因二:在一開始的時候就產生問題,可能有泡泡跑過定序的夾層,這時候就不適合trim,會因此拋棄過多好品質的序列
        • 原因三:太低的coverage可能會造成false negative

Per Base Sequence Content序列中不同位置其ATCG所佔的比例

screenshot.png

    1. 分析內容:
      1. 在一個隨機的library中,其序列中不同base應該不會有太大的差別.
      2. 通常在開頭的地方容易有不平均的ATCG比例,通常在Library製備過程中priming使用到random hexamers(幾乎所有RNA-seq的製備都會使用),或是在fragmentation過程中使用transposases,通常會在reads開頭有偏差出現,使得在有enrichment在5’端的different k-mers處,這種技術性造成的偏差無法藉由trimming消掉,但在downstream analysis中也不會造成太大的問題
      3. 達不同品質指標,圖標會顯示之圖例
        • 警告(Warming):當任何一處的ATCG其差異超過10%
        • 失敗(Failure):當任何一處的ATCG其差異超過20%
      4. 問題導因:
        • 原因一:Overrepresented sequence
          1. 可能是做定序所使用的adapter dimers或是rRNA
        • 原因二:Biases Fragmentation
          1. 因製備Library中所使用的random hexamers在前12bases會有selected biase,但這在後續的process可以處理掉,不太會影響後續的分析
        • 原因三:Biases Composition Libraries
          1. 大多數的library本身就會有些composition的偏差,尤其是使用sodium bisulphite製備的,其會造成cytosine變成thymines,所以會使得C的比值偏高
        • 原因四:Extreme Trimming
          1. 在尾端會出現,因為很積極地將跟adapter match的序列去除掉

Per Sequence GC Content 測量GC在序列中的比例
screenshot.png

    1. 分析內容:
      1. 其主要是分析不同序列的GC之分布圖,應該要跟理論分布差不多
      2. 比較陡或是尖:可能是被adapter汙染
      3. 比較寬:可能是被不同物種所污染
    2. 達不同品質指標,圖標會顯示之圖例:
      1. 警告:偏離理論值15%
      2. 失敗:偏離理論值30%

Per Base N Content 平均每個位點壞品質(N)

screenshot.png

    1. 分析內容:
      1. 當某個序列即位點定序品質很差,則會用N代表
      2. 通常當某個部分其N所佔的比例上升的話,其會使後面的分析無法獲得有效的base cal.
    2. 達不同品質指標,圖標會顯示之圖例:
      1. 警告:任一個位置其N所佔比例大於5%
      2. 失敗:任一個位置其N所佔比例大於20%
    3. 備註原因:
      1. 原因一:
        • 通常是整體定序品質不佳,可以檢查一下特定bin的coverage,可能有某一個bin其定到的序列非常少
      2. 原因二:
        • 有嚴重偏差的LIBRARY,在Per Base Sequence Content也可以看出來

 

   Sequence Length Distribution 定序長度分布圖
screenshot.png

    1. 分析內容:
      1. 顯示平均定序的長度,這部分就要先了解是使用何種定序方式,有的會有單一長度,有的則是不同
      2. 達不同品質指標,圖標會顯示之圖例:
        • 警告:假如定序的長度有不一樣的話
        • 失敗:任一個定序長度是零

Duplicate Sequences 重複序列
screenshot.png

    1. 分析內容:
      1. 在一個均勻的library中序列應該只會出現一次,很低程度的重複率可能暗示高coverage,高度序列重複則可能是library中有某種程度的enrichment bias
      2. 為了使分析的計算效率提高,這個分析模組只會分析在前100000的序列,另外在超過10個重複的序列,其會被放在一起呈現
      3. 比較長的定序讀數,通常會造成一定程度的低估重複率
      4. 藍色線:代表完整序列其重複率, 紅色線:代表去掉重複後不同序列的比例
  1. Overrepresented Sequences
  2. Adapter Content
  3. Kmer Content
  4. Per Tile Sequence Quality

去掉品質差Reads的工具：Trimmomatic

可以用來處理illumina機器裡的adaptor，其也可以將low-quality的base或是N base去掉，其可以處理single-end和paired-end的資料。此程式是java base的軟體，在command line下操作。

 

處理特殊設計的paired-end 實驗設計：BBMerge 和 FLASH

有些實驗設計時，會將paired-end read設計成疊在一起的樣子，而這兩款程式可以將重疊的paired-end reads合成單一條較長的reads，且能計算出innert-size的分布圖。

偷偷查了一下FLASH的使用趨勢，似乎暴增人數使用！
screenshot.png

RNAseq可用來研究的主題

發表留言

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]

RNAseq的資料為high resolution到single base，可應用範圍寬廣

RNAseq的技術可以用來區分許多不同種類的transcripts像是long non-coding RNA、miRNA、siRNA、snRNA、piRNA，以及他們如何調控蛋白質的translation或是調控chromatin states

RNAseq也可以用來研究enhancer region的transcripts表現量，這可能跟epigenetic gene regulation相關，可以稱這類型的transcripts為enhancer RNA。

從RNAseq的資料可以取得transciption start sites、alternative promoter usage、mRNA isoforms (alternative splicing)、mRNA 穩定性（premature transciption termination at the 3′ ）

Locating regulatory element

Allele-specific expression、disease-associated single nucleotide polymorphisms(SNP)、fusion gene for causal variants in cancer

Transcription of endogenous retrotransposons and parasitic repeat elements that influence the transcription of neighboring genes or result in somatic mosaicism in the brain