小耿的研究筆記

2014年2月20日星期四

Blast 單機版語法

最近在用blast比對TM protein database,
但下載了win單機版的blast[2.2.29+]版本後卻苦於不清楚指令囧,

參考這篇文章 終於找到方案了,
而更詳細的blast指令說明文件可以參考官網的資料
(網址:http://www.ncbi.nlm.nih.gov/books/NBK1763/)

Blast比對的概念是這樣的,你輸入一個要比對的資料庫然後系統會跟預先建好的參考資料庫作比對,所以進行比對前要先建立好自己的參考資料庫。

而實際執行的方法是這樣的:

0. 很重要,不要再點那個執行檔了,沒有反應的,請先開一個命令列視窗然後轉移到安裝blast的資料夾,後面以 ...\blast---\bin\ 代表,而所有的檔案都請放在這個資料夾,之後的輸入命令就在命令列視窗進行 (我居然困在這一步兩天阿囧rz)。

1. 要進行比對前要有一個比對的資料庫,其指令為

makeblastdb -in [輸入檔案] -dbtype [prot 為蛋白質資料庫; nucl為核酸資料庫] -out [輸出資料庫名]

例如:要從 seq_test.txt 這樣一個蛋白質序列資料庫檔生成一個叫 ooxx 的參考資料庫, 只要開命令列視窗在blast下的...\blast---\bin\ 中輸入

makeblastdb -in seq_test.txt -dbtype prot -out ooxx

執行就會根據seq_test.txt這個序列資料跑出 ooxx.phr ; ooxx.pin ooxx.psq 三個檔案

要注意的是seq_test.txt 的格式要是fasta檔, 我用的格式是由matlab讀檔生成,一行是檔名,下一行沒有分隔符號,全序列結尾也沒有空行,然後也放在...\blast---\bin\資料夾

fasta例:
======================================
>1ppj_I
MLSVAARSGPFAPVLSATSRGVAGALRPLVQAAVPATSESPVLDLKLSVLCRESLRG
>1ppj_H
GDPKEEEEEEEELVDPLTTVREQCEQLEKCVKARERLELCDERVSSR
>1ppj_G
GRQFGHLTRVRHVITYSLSPFEQRAFPHYFSKGIPNVLRRTRACILRVAPPFVAFYL

======================================

然後要做blastp的時候呢.. 可以用下面的指令就可以跑出來了

blastp -query [輸入檔] -db [參考資料庫名] -task [指定] -evalue [Evalue值] -outfmt [Num] -out [輸出檔名]

要注意的是參考資料庫名不需輸入副檔名, 然後outfmt是輸出格式0~9格各有不同格式如下表格

而 task可指定一些特別的功能如:

例如: 針對一個叫 seq_all_test.txt 的待測序列檔使用TM_5_test這個資料庫測試, 並指定對短序列最佳化,並設 Evalue>10^10 輸出格式為含表格介紹的表格(7) 存成ooxx.txt檔, 指令可以這樣下

blastp -query seq_all_test.txt -db TM_5_test -task blastp-short -evalue 10000000000 -outfmt 7 -out ooxx.txt

輸出的ooxx.txt檔的形式為

================================================

# BLASTP 2.2.29+

# Query: 1a0t_P

# Database: TM_5_test

# Fields: query id, subject id, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score

# 23 hits found

1a0t_P 1a0t_P 100.00 413 0 0 1 413 1 413 0.0 867

1a0t_P 1af6_A 27.05 451 257 22 4 413 2 421 2e-019 80.8

1a0t_P 1bcc_A 62.50 8 3 0 273 280 371 378 12 15.5

1a0t_P 1bcc_A 47.37 19 4 1 59 77 19 31 41 13.8

1a0t_P 1bcc_A 21.11 90 48 4 91 157 205 294 48 13.4

...

=================================================

2014 0220 小耿

註:目前跑出來的結果好像會把一些序列切成很多段分開比對, 我還找不到怎麼把兩個序列完整比對的參數 QQ

========
% 2014 0413加註:

如果對於單機版的blast使用不熟練也可以去myblast網站,那邊能夠線上執行blast,
(myblast 網址: http://mybioweb.nhri.org.tw/myblast/ )

只要註冊帳號後就可以用視窗界面上傳參考與目標資料庫即可進行大量序列的比對(當然使用前還是要先了解blast比對的參考與目標資料庫的概念跟格式就是),

myblast的優點是沒有blast 2.2.29版的序列切割問題,數據呈現較於簡潔,但是要注意的是至少在2014 0220時測試的結果該網站使用的還是舊版的blast, 版本為2.2.19。

而當下載大量的比對資料後,要把資料讀進matlab生成E value矩陣運算時,未免一次讀進太大的檔案造成matlab當機,我針對myblast跑出來未設定格式的輸出結果寫了一個逐行讀檔的讀取e value程式,可以在此下載( http://1drv.ms/1kTVQ5D ), 使用時用matlab打開該m檔案然後根據m檔案前面的說明1.設定要運算檔案的正確路徑後(預設資料夾為C:\blast-2.2.29+\bin\),執行即可自動抓檔讀取E value跟序列順序名單,輸出結果為:

序列名單: IDname_meta
E value矩陣: E_value_matrix

要注意的是此處為未處理的E value 矩陣,此處的E value不合於距離定義下該符合的 Eij =Eji, 為了達到此步驟我設定的修改後的E值為 E'ij=(Eij*Eji)^0.5 。

2013年12月19日星期四

[簡報]Ptt 推文資料庫可能性(2013 12 17)

此檔案為 2013 1217 在台大物理館 R516 的非正式討論會用的投影片檔，粗略的討論ptt推噓文系統所具有的價值跟可能性。

在Blogger上播放時，動畫與字體部分會出問題，有興趣的可以根據以下連結下載檔案
: http://sdrv.ms/1kXi2KV

2013年11月19日星期二

[小程式] 用 matlab 抓取 ptt上特定文章推文的時間序列資料 (Catch the "push" and "spurn" index time series data from a article in PTT by matlab )

動機:

ptt 是目前中文文化圈中最大的電子佈告欄，平均來說隨時隨地都有10萬名左右的使用者在線上瀏覽、發文、或參與討論。其討論區涵蓋的領域包羅萬象，幾乎你想的到的主題都有其討論區，沒有的話也可以自己申請創一個。而裡面大多數的發文以及後續討論資料皆會以純文字檔的形式被保存下來，久而久之其形成的巨大資料庫確實是一個分析以台灣為主的中文文化圈之社會行為的好資料。

推噓文系統則為ptt電子布告欄裡面獨有的文章底下評論模式，藉由推文可以肯定這篇文章的價值，而噓文則多半代表不認同 (唯此兩種情況又因各不同討論區文化有異)。但整體來說，一篇獲得眾多推文噓文或註解的文章可約略代表是一篇受到大眾所關注的資訊，較多推文的文章受到較多人賦予正面評價，較多噓文或註解的文章受到較多人賦予負面或持平評價。因此一篇文章的推文噓文或註解的數目嚴然可以成為一篇文章受到多大程度關注的指標。

而進一步來說，推噓文與註解隨著時間而變化的動力學過程，或許也可以展現出一則具有某種特色的資訊如何自某個地方發生後透過社會上人群的隨機關注語互相影響而漸漸發展成一個廣為人知的訊息。這些資料對於探討以台灣為主的中文文化圈中資訊如何透過社會發展自我關注度應該會是很有幫助的數據。

在2011/10/12的時候我曾經手動去統計ptt上面熱門排行版第一名的八卦版( Gossiping )上面當天全版的文章數目與推文數目對應的關係，發現文章數目與推文數之間有頗為明顯的乘冪(power law)關係，可很粗略的揭露出在一個最熱門看板中平均來看，一篇文章要取得高推文數的難易度大概如何(從關係可粗估推文數目增加兩倍的文章數目會是原先的1/5~1/4倍)。

但這樣粗糙的結果尚無法看見更精細的東西，若能取得單篇文章的推文噓文與註解的時間序列資料則可以進一步看到一篇文章受關注程度隨時間演變的資訊，或許就有更多可供分析的原始數據。

而今天中午因為新桌機正在進行持久的硬碟格式化過程，筆電也快跑不動目前的分析結果，該查的資料又毫無進展，一個翻桌之下我就來寫這個程式轉換心情了。

從今天中午寫的這個小程式可以藉由輸入ptt某篇文章的網址而得到該篇文章的推噓文時間序列資料，若進一步的運算與分析應可看到更多性質，而或許也可以根據這些數據的基本結果建立資訊受關注程度的時間演化模型來討論一則訊息是如何受到關注的唷!

然後這個大概會慢慢的玩，如果有人有興趣也可以來找我討論討論唷 ^^"

檔案與連結:

檔案名稱: push_time_sequence_data.m
檔案路徑: http://sdrv.ms/17IdlPe

使用說明:

0. 打開matlab,開啟本m檔案.
1. 請抓取目標文章網址(可在ptt介面下文章前用大寫 Q取得).
2. 將文章網址貼於下方 urlpath=後面, 記得網址前後要用 ' 符號包夾 ( 例: 'www.google.com.tw')
url_path='http://www.ptt.cc/bbs/StupidClown/M.1233060106.A.55E.html'

3. 儲存檔案後按F5開始跑
4. 跑完後資料存於tPEN_data檔:
tPEN_data(:,1) 為時間軸, 目前設定時間單位為小時(hour)
tPEN_data(:,2) 為推文資訊, 該時間有推文則值為1,否則為0
tPEN_data(:,3) 為註解資訊, 該時間有註解則值為1,否則為0
tPEN_data(:,4) 為噓文資訊, 該時間有噓文則值為1,否則為0
5. 資料取得完成

2013年7月13日星期六

小熊餅乾搖搖樂動態過程觀察實驗

原文網址: http://pipips.blogspot.tw/2013/07/blog-post.html
==========================================

這次實驗使用的是水源校區即品網買的牛奶口味小熊餅乾,因為沒拍照故取用網路圖片，
圖片來源: http://garywlee.blogspot.tw/2013/05/blog-post.html 。

實驗動機:

前陣子很流行小熊餅乾搖搖樂新吃法(他人的範例網址)，但是網路的影片幾乎都是放在小熊餅乾盒裡面搖，看不到裡面的變化情況，只能了解小熊餅乾實驗的初態跟末態，因此讓人很好奇整個過程中有怎樣的機制跟動力學的變化。然後剛好昨天的蘇力颱風預備存糧中有在水源校區即品網買的牛奶口味小熊餅乾，下午去聚餐的空檔突然想到這個有趣問題就馬上來做實驗了! ~~決不只是為了無聊打發時間而已唷~~。~~

實驗材料:

1. 小熊餅乾牛奶口味一盒 (水源校區即品網購得)
2. 蜜桃紅茶玻璃罐一個 (台北火車站爭鮮壽司購得並當水杯一陣子)

3. 單純拿來拍照用的難以啟齒之韓國S牌智慧型手機 (好像該換了?)

實驗方法:

為了解決搖動時看不到小熊們怎麼變化的缺點，並提升震動腔碰撞時能量傳遞的效率以節省實驗消耗的能量(or降低手部肌肉痠痛指數)，而選擇將小熊餅乾放進透明的玻璃罐來進行實驗。

接著為了紀錄過程，每搖動一百下(用手上下搖動，頻率約 100 次/分)就停下來拍一次照片當作實驗記錄，由於本實驗出於臨時性的興趣因此並沒有採取定性定量跟嚴謹的紀錄過程，僅呈現較粗略的結果以供大家參考。

實驗記錄:

1. 搖動次數零下，小熊們無辜的躺在玻璃罐裡不知道等等會發生什麼事情。