samtools的功能非常強大，也有很多用法，這篇有概略介紹，而這邊仔細講解一下比較細緻關於alignment檔案的extraction，就是檢視特定區域的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
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的方式串接起來
第二步，會使用這種graph產生出isoforms的可能組合
–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較多！
notes from the GSB Bioinformatics study group
相對於Next Generation Sequence 資料動輒幾十gb，這次我們選擇用microarray的資料練習，
首先我們分別從GEO、EMBL和TCGA裡找尋用affymetrix 測expression的實驗資料，畢竟affymetrix是microarray年代裡最可靠的平台，另外我們也希望使用的chip 是一樣的，避免到probe annotation的問題，最後我們找到三組可行的資料來練習分析，都為使用Affymetrix Human Exon 1.0 ST array(GPL5175)：
這組資料來自於台灣長庚醫院頭頸癌團隊的研究論文，此組包含有79個樣本，其中22個為正常組織
這組為來自中國大陸的研究團隊，其基因表現資料來自於27個paired胃癌組織樣本
這組來自於挪威的研究團隊，主要有207組來自於挪威的大腸直腸癌病人檢體，其臨床分期有第一期到第五期，實驗主要是用來找尋第二期大腸直腸癌是否有可以用來預測預後的基因表現biomarker
參考此篇來做匯入
先安裝GEOquery R bioconductor package
source(“http://bioconductor.org/biocLite.R")
biocLite(“GEOquery“)library(GEOquery)
接下來使用GEOquery來下載我們找到的檔案，先嘗試GSE25099
gse <- getGEO(“GSE25099“, GSEMatrix = TRUE)
接下來我們只要這資料裡頭的gene expression matrix，於是我們把exprs subsetting出來，因為這是s4 class，先用＠在用$，交替尋出我們要的matrix
接下來快速使用R內建函數來做簡單分析，可參考這篇，使用prcomp()
samtools是在inspect或是處理bam/sam黨非常重要的工具，這邊簡單介紹他的功能和基本語法使用。
基本上，samtools對於sam/bam檔案有三大功能：
第一 ，就是檔案格式的轉換（bam to sam/sam to bam）
第二，sam/bam reads的篩選，過濾、查看
第三，進階的variant calling分析（mpileup）
一 、檔案格式的轉換我們可以使用samtools將bam檔儲存成sam黨，反之亦然，雖然tools一開始的defaut是input bam，但都是可以調整的。
bam 轉 sam: $samtools view input.bam >output.sam
sam 轉 bam: $samtools view -b input.sam > output.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疊加數
＄samtools idxstates data.bam
黃色代表是左邊染色體上有被mapped上的reads數量，右邊代表的是unmapped量
這邊就需要調整比較多參數，其實這些功能蠻方便的，可以快速地瞭解整體bam檔的狀態！
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 壓縮的模式，如下圖所事是：
簡單來說，就是以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的檔案架構分兩部分，一部分是Header，另一部分是Alignment。外觀如下圖：
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的狀況：
依序是每個Reads獨一無二的編號，通常還有在定序機器上的位置，FLAG，則是代表此reads的mapping狀態，chr則很明顯的是對到的染色體，Start代表reads對應到reference上的開始位置，CIGAR則代表這個reads每一個位置跟reference mapping的狀態，是否100Match到，後面還有一些進階資訊，假如是paired-read則會顯示他的mate reads的位置。
因為sam檔沒有壓縮過，所以大多數會以bam檔形式儲存，必須使用samtools來轉換，以便閱讀上面那些資訊，這邊samtools在處理bam/sam檔案上就非常重要了！
Introduction
多變量統計分析相對於一般的(單一變量)統計分析更為進階，探討的層面也更為多元。學習正統的多變量統計分析，也許需要先修一些數學的工具，如線性代數、矩陣代數等等。然而並不是每個來自各個專業領域想要運用多變量統計分析的人都能熟悉那些生硬的數學工具，但它仍是社會科學、商學和生物學等領域重要的研究工具。 繼續閱讀 “Basic Multivariate Statistics and Principal Component Analysis (PCA)"
論文閱讀：Detecting differentially expressed genes by smoothing effect of gene length on variance estimation
表現差異分析(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次）
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 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）：
簡單來說，每個研究都會有個計畫編號，開頭大多會是PRJN或是SRP,這種計畫通常都會不只做一組檢體，而每個送去做定序的檢體都會有專屬編號，以SAMN開頭，而此檢體使用的定序實驗設計，也會有個編號，以SRX開頭，裡面會記錄此檢體所使用的定序機器種類、library type、定序設計（paired-end, single ….)，而在同一次定序中其跑的run也會有編號，以SRR開頭。
接下來寫一下，SRA toolkit的下載和使用：
第一步：到此頁面下載符合自己系統的版本 link
第二步：在系統下解壓縮
這邊以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檔案，這邊是其指令介紹：
主要指令依功能分成三大類：
閱讀自Reproducible, Scalable Fusion Gene Detection from RNA-Seq.Arsenijevic V1, Davis-Dusenbery BN2.
在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的機率上升！
我們的基因體時代 Our "Gene"ration 