給生物學家的計算分析(六):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

癌症基因體閱讀心得(一)

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

給生物學家的計算分析(五):shell下的文本處理

閱讀自Practical Computing for biologists by Haddock Dunn

shell環境下的文本處理效率相當高,但需要combine許多shell command line tool,且同樣的處理其實方法會有蠻多種的,可能隨者經驗增加或是使用的工具變多而調整,這邊就稍微介紹一些基本的工具,已可做非常多的應用。

shell環境裡的文本編輯器nano

在command line底下編寫文字,需要一個方便的工具,雖然之前介紹過圖形介面的TextWrangler,但直接在command line裡面編輯有時候比較方便,其實在unix核心下的系統(linux, OS X)有蠻多內建選擇的,如vi、 vim 、emac、pico,每種操作細節上都有一點不同,而這邊介紹nano,一款比較基礎的,可以直接在command line下呼叫 nano
螢幕快照 2016-01-16 上午7.35.37螢幕快照 2016-01-16 上午7.35.48

編輯器下排其實就有非常清楚的解說,關於如何離開、寫入等等,如輸入:

^X   離開
^O   寫入
^Y    向上翻頁
^V    向下翻頁

 

 

給生物學家的計算分析(四):進入shell

閱讀自Practical Computing for biologists by Haddock Dunn

command line的工作環境是一開始生物背景跨進生資分析的障礙之一,相對於已經被圖形介面(Graphical User Line)的電腦環境所嬌生慣養,要說服我們從舒適的圖形介面跳進黑壓壓的指令框框,想必要有強而有力的說詞,這邊就要說說圖形介面的壞處:
  • 生資分析常常要做複雜的步驟來處理資料,或是同樣的步驟操作多次,這些工作要是再圖形介面下不知道要用滑鼠點到什麼時候。
  • 一般在GUI的環境下,大部分的分析都不會留下log黨,讓我們去回顧可能是哪一個步驟做錯了,而在command line下則很方便能把不同分析步驟發生的事情紀錄在log中,等出問題的時候方便除錯
  • 叢集電腦系統不支持圖形化介面的操作,當我們需要使用更加計算能力時,使用圖形介面會是個問題
  • 對於軟體開發者來說,圖形介面的軟體開發複雜且不同平台使用會需要調整,所以在軟體更新上,command line操作環境下的軟體更新的較快!
以上的問題,其實很明顯地想要處理大型資料就必須要克服這個恐懼感。
基礎的unix command 知識:
  • 相對路徑和絕對路徑
  • 相對路徑的使用:
    • ~   代表home目錄
    • ..  代表前一個目錄
    • .   代表目前目錄
  • cd 可以用來在系統路徑中移動
  • 基礎指令:
指令
功能
ls
當下目錄的檔案
pwd
目前目錄之路徑
mkdir
建立資料夾
rmdir
移除資料夾
rm
移除檔案
less
看檔案
man
查看指令說明
cp
複製檔案
mv
移動檔案或是重新命名

給生物學家的計算分析(三):Regular Expression使用

思考如何寫Regular expression的thinking process
以下是資料
螢幕快照 2016-01-16 上午8.01.22
首先思考要如何整理,整理成什麼樣子
螢幕快照 2016-01-16 上午8.04.27
再來看目前資料每行的狀態
螢幕快照 2016-01-16 上午8.06.46
先看哪些是要的資訊,用(     ) 標示
(13) (Jan)uary, (1752) at (13):(53)    (-1.414)    ( 5.781)etc
先看括號裡的資料怎麼用Regular Expression表示
螢幕快照 2016-01-16 上午8.18.53
再將括號外的符號放進去
(13) (Jan)uary, (1752) at (13):(53)    (-1.414)    ( 5.781)etc
(\d+)\s+(\w{3})[\w\,\.]*\s+ (\d+)\sat\s (\d+):(\d+)\s+    ([-\d\.]+)\s+    ( [\d\.]+).*
    1               2                             3               4          5               6                         7
在調整取代裡面的括號資料順序
(13) (Jan)uary, (1752) at (13):(53)    (-1.414)    ( 5.781)etc
(\d+)\s+(\w{3})[\w\,\.]*\s+ (\d+)\sat\s (\d+):(\d+)\s+    ([-\d\.]+)\s+    ( [\d\.]+).*

    \3               \2 \.              \1               \4                 \5              \ 6                         \7

 

*:(quantifiers)至少符合0個
+:(quantifiers)至少符合1個

論文閱讀:Differential gene and transcript expression analysis of RNA-seq experiments with Tophat and Cufflinks

原始論文:

Cole Trapnell,Adam Roberts,Loyal Goff,Geo Pertea,Daehwan Kim,David R Kelley,Harold Pimentel,Steven L Salzberg,John L Rinn& Lior Pachter. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks.Nature Protocols.7,562–578(2012)Doi:10.1038/nprot.2012.016

 

RNAseq的資料分析主要有兩個大目的:

  • 找尋新的基因或是transcripts
  • 比較不同狀態下的基因表現亮

而在分析步驟上,可以把它拆成三個分析流程:

  • read alignment
  • transcript assembly 或是 genome annotation
  • transcript and gene quantification

本篇論文主要介紹用期開發的程式作為workflow,並且簡單介紹每個tool之間的異同。

所使用的工具和工作流程:

Read alignment with Tophat
Reads拼回去reference的步驟非常重要,他可以找出所有可能的insertion/deletion/indel,這些資訊可以讓我們了解辨識出相對於對照序列的polymorphism,另外,對不回去的reads,可能是新的protein-coding genes或是noncoding RNA。這邊拼回去的結果,也會影響之後如何去計算transcript abundance的數值。

螢幕快照 2016-01-14 下午3.50.21

 

Transcripts assembly with Cufflinks
為了要知道各種基因的表現量,我們必須知道各種reads所對應的isoforms transcripts,這邊就必須把Reads拼組成各種特定的transcipts,這步驟牽涉到如何辨別不同種transcripts variant ,所以最後產生transfrags來表達所有可能。

cufflinks會將一些low abundance的tansfrag當作是來自於immature transcripts而丟棄,這部分是要了解的。(或許這步驟可以思考一下!),有時候我們會將所有的bam pool在一起讓他找transcripts,但這種作法其實一方面增加電腦計算的loading,另一方面,則是會讓程式在計算isoform時增加很多變數,所以cuffmerge應映兒生,他可以把individual產出的transfrag merge在一起看!

螢幕快照 2016-01-14 下午4.14.47

Differeial Expression Analysis with Cuffdiff

cuffdiff可以輸出好幾個檔案描述不同基因和transcripts的表現量差別和p value以及其名字和在基因體上的位置,另外,cuffdiff也可以把同一組基因裡面不同TSS的transcripts分組去計算他們之間比例的變化!

 

視覺化資料CummeRbound

cuffdiff提供基因和transcripts的表現量分析,且這些資料是以tab-delimited的方式輸出很方便後續的分析,而CummeRbound則是提供一個更容易操作的視覺化分析,可以將cuffdiff的計算結果做更細緻的分析。

 

替代工具:

  • Read-alignment 
    • GSNAP
    • START
    • Map-Splice
  • Transcriptome reconstruction
    • De novo transcripts assembly
  • Quantification
  • Differential expression

有趣的相關論文:

Read-alignment

Transcript reconstruction

Quantification

Differential Expression

 

給生物學家的計算分析(二):用Regular Expression來搜尋及處理文本

閱讀自Practical Computing for biologists by Haddock Dunn

給生物學家的計算分析(一):文本乃資料分析之底

Regular Expression是什麼呢?對於生物背景的人來說,一開始對於這名詞總是霧茫茫的。想要深入研究的話,可以看這篇文章 。簡單來說,Regular Expression(正規表達式)是用來標記某種特定文字pattern的符號表達方式,讓我們能用更強大的方式搜索文本並且修改,要注意的地方是不同程式之間所使用的Regular Expression會有一些差異,這邊想要介紹是Regex

Regular Expression在表達一個搜尋pattern時,主要由三類型符號組成:

  • Anchors :代表要搜尋字串的相對位置
  • Character Sets:在此位置上的符號是什麼
  • Modifiers:多少次符合上面Character Sets的表達

威力強大的工具:搜尋+取代

搜尋字彙(character sets)
意思
\w
可代表文字、數字或是底標
\t
代表Tab符號
\s
代表空白鍵、Tab、和end of line
\r      \n
代表end of line
\d
代表數字0-9
.
代表任何符號除了end of line
Anchor
意思
第一個符號
最後一個符號
Modifiers
意思
[   ]
客製化wildcard
{   }
創建一個array
(    )
指定此內符號被置換

來看一下實際的例子

RE.001

論文閱讀Onco-proteogenomics:Cancer proteomics joins forces with genomics

Onco-proteogenomic research 是結合genome和proteome的資訊來回答癌症研究的問題。

在生物的中心法則central dogma: DNA -> RNA -> Protein,越往下越直接影響到生物體的功能,但也越難研究,自從Mass Spectrometry的技術不斷精進,目前研究蛋白體學的武器比之前跑電泳的時代不同了,但研究蛋白體學的技術難度終究還是遠大於基因體學的技術,隨者定序資料越來越便宜,使用定序資料來幫助蛋白體學的研究是另一個蠻新穎的方式,目前RNA-seq可以很容易找出大量的unknown transcripts或是chimeric transcripts,但是不是真的存在,是很難單純用定序資料可以回答的,要是能使用Proteomic data來validation,那麼會非常有說服力。

螢幕快照 2016-01-13 下午12.26.23

 

Onco-proteogenomics研究有許多待克服的挑戰,以下簡單介紹和搭配各自的reference:

挑戰一在癌症研究中,我們最在意的便是能否找到一些在癌症組織中基因、蛋白體的特異性改變(Tumor-specifci change in the proteome),像是:癌症是怎麼開始的(Tumor initiation)、癌症為何會惡化(Tumor progression)、癌細胞是如何對治療產生resistance(Adaptation to treatment),這過程中,研究者很容易會找到變異,但哪些是“passenger change”,哪些是“driven change”,背後要用什麼統計模型來解釋、如何處理多樣化的資料來源、資料如何呈現等等

挑戰二想要用MS/MS來看癌症組織中某個蛋白質或是變異蛋白質的量,來確認從Genomic data analysis的結果,也會遇到“不一定看得到想看的,沒看到也不代表沒有”的難題,畢竟蛋白體在細胞中的量是動態的,而一些基因的變異所產生的mutated protein產量不多(基因變異不一定都造成蛋白質的量下降,像是P53 mutated會造成相反的效果)

 

挑戰三從genomic data建置成reference data base時所使用的six-frame translation會產生比較多的可能,但實際上可能產生突變的蛋白質約莫數百,且大部分都是passenger mutation,從統計來看可能會造成false discovery rate上升。

挑戰四在proteogenomic 研究上有一個根本的問題,那便是如何將不同的資料結合,來闡述所有跟cancer-related phenotype的認知,這是非常困難的,但才是我們核心想知道的,到底病人怎樣才會回復健康。

挑戰五目前評估某一個新的Onco-proteogenomic技術都是看其發現多少基因變異為主這種思維很容易造成研究人員過度的追求找出新的mutation,但另一方面同時也是增加false discovery rate,造成一堆垃圾資料產生。

 

幾篇有趣的Onco-proteogenomic research paper:

Evans, V.C. et al. De novo derivation of proteomes from transcriptomes for transcript and protein identification. Nat. Methods 9, 1207–1211 (2012).

Li, H., Ruan, J. & Durbin, R. Mapping short DNA sequencing reads and calling variants using mapping quality scores. Genome Res. 18, 1851–1858 (2008).

Aquino, P.F. et al. Exploring the proteomic landscape of a gastric cancer biopsy with the Shotgun Imaging Analyzer. J. Proteome Res. 13, 314–320 (2014).

O’Rawe, J. et al. Low concordance of multiple variant-calling pipelines: practical implications for exome and genome sequencing. Genome Med. 5, 28 (2013).