Principle Component Anayasis練習:口腔癌、胃癌和大腸直腸癌基因表現資料(一)

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

samtools 使用

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

常見的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中的表現差異分析呢?

論文閱讀: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次)

如何使用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: Normalization

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中的表現差異分析呢?

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

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