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

 

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