月份: 一月 2016

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的使用！

Basic Multivariate Statistics and Principal Component Analysis (PCA)

5 則迴響

Introduction

  • 多變量統計分析 | Multivariate Statistics

多變量統計分析相對於一般的(單一變量)統計分析更為進階，探討的層面也更為多元。學習正統的多變量統計分析，也許需要先修一些數學的工具，如線性代數、矩陣代數等等。然而並不是每個來自各個專業領域想要運用多變量統計分析的人都能熟悉那些生硬的數學工具，但它仍是社會科學、商學和生物學等領域重要的研究工具。 繼續閱讀 “Basic Multivariate Statistics and Principal Component Analysis (PCA)" →

基因的長度是否會影響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: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: 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中的表現差異分析呢？

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