R 語言 判斷式語法

「短路」（short-circuit）運算子 &&（AND） 與 ||（OR） 常常用於 if 判斷式的條件控制部分，這裏要注意 & 與 | 是用於向量的所有元素，而 && 和 || 只能用於長度為 1 的向量，且第二判斷式只有在必要時才會被執行：&& 運算子若是發現第一個條件傳回的是 FALSE，那麼就不會執行第二個條件的的判斷；|| 運算子也有類似的情形，若是發現第一個條件傳回的是 TRUE，那麼就不會執行第二個條件的的判斷


KRHH <- repulsive[ ( ( V4 == 'KR' &  V5 == 'HH' ) | ( V4 == 'RK' &  V5 == 'HH' ) ) &  V6 < 15, ]

經典詞句

當斷不斷   反受其害

不怕狼一樣的隊手   就怕豬一樣的隊友

不可能三頭六臂做到面面俱到  讓每個人都喜歡你   又不是新台幣

無間道   我只想做一個好人  為什麼你們都不給我一次機會？


給錢的時機是學問
完事後給的是嫖娼
一周後給是性夥伴
按季定量給是包養
全年度都給是二奶
終生不給的是老婆
從不給是紅顏知己
總拖著不給是無賴

Dataset

1. Different packing of external residues can explain differences in the thermostability of proteins
from thermophilic and mesophilic organisms
A unique database of 373 structurally well-aligned protein pairs from thermophilic and mesophilic organisms is constructed.  

Link

2.

Pymol WCN color

用  WCN 算出 CA atom 值後
取代原有 bf values

利用 pymol 呈現出 WCN 差異性

參考這裡  Link
command line :

Assign color by B-factor

B-factor coloring can be done with the spectrum command. Example:

spectrum b, blue_white_red, minimum=20, maximum=50 

as cartoon 

cartoon putty

ps1. backgroud  Link

set ray_opaque_background, off 

ps2. pymol tips and tricks Link

ps3. Git Link 

ps4. http://kpwu.wordpress.com/?s=pymolhttp://kpwu.wordpress.com/?s=pymol 

hmmer

1. web download and install >> link


2. port install  >> port install hmmer @3.0

Align

Link  Structure

FAST is a program for aligning three-dimensional protein structures. mac download


Dali :  local run.   Link   


Link  Sequence

EMBOSS link

亦可由 port 來安裝


Needle:以Global Alignment進行序列的兩兩比對

		needle操作步驟：
			Go to ： needle
			輸入序列，有3種方法可供選擇：
			1.	由 file /database：鍵入database：accession number。
			2.	copy /paste：直接貼上序列。
			3.	擷取PC中的序列檔案：點選 "List of files"， 再將您電腦中適當的序列檔案輸入即可。
			以下列二條序列為例，選定由 file/database 輸入 輸入序列(1) embl:HSPDGA4 輸入序列(2) embl:HSA238420
			在 required section 設 gap opening penalty =10, gap extension penalty=0.5
			在 Advanced section 選擇scoring matrix (Blossum 62)
			在 out put section選擇輸出格式-- (default 為srspair)
			注意：needle程式一定要輸入兩條序列；並請記得在執行區右下角選擇batch (背景執行)
			按 Go 。
			檢視結果: 點選file, 選擇saved results, 選取本次分析結果，按display。 檢視結果後存檔：在結果畫面左上方點選file, 選擇save to local file, 輸入正確檔案名稱, 按「儲存」，例：存成hspdga4.needle。
			比對結果： 由於 "embl:HSPDGA4" 是 PDGF A chain exon4 的序列，而 "embl:HSA238420" 是 PDGF A chain gene exon1~4 的序列，在結果中很清楚的看到exon 4 對應在A chain 的第7758 至8005位置上。
		Water程式: 以local Alignment進行序列的兩兩比對
		water 操作步驟：
			Go to ：water
			輸入序列，有3種方法可供選擇：
			1.	由 file /database：鍵入database：accession number。
			2.	copy /paste：直接貼上序列。
			3.	擷取PC中的序列檔案：點選 "List of files"， 再將您電腦中適當的序列檔案輸入即可。
			以下列二條序列為例，選定由 file/database 輸入 輸入序列(1) embl:HSPDGA4 輸入序列(2) embl:HSA238420
			在 required section 設 gap opening penalty =10, gap extension penalty=0.5
			在 Advanced section 選擇scoring matrix (default=Blossum 62)
			在 out put section選擇輸出格式-- (default 為srspair)
			注意：請輸入兩條序列，並記得在執行區右下角選擇batch (背景執行)
			按 Go
			檢視結果: 點選file, 選擇saved results, 選取本次分析結果，按display。 檢視結果後存檔：在結果畫面左上方點選file, 選擇save to local file, 輸入正確檔案名稱, 按「儲存」。
			比對結果： 由於 "embl:HSPDGA4" 是 PDGF A chain exon4 的序列，而 "embl:HSA238420" 是 PDGF A chain gene exon1~4 的序列，在結果中很清楚的看到exon 4 對應在A chain 的第7758 至8005位置上，由於是local alignment，因此不match的序列就不列出來了，這也正是local alignment和global alignment最大的不同之處。

Pymol

Link

>> set ray_trace_mode, 2

press "Ray" to show it like under that ! 

# Black and White Script
load /tmp/3fib.pdb;
show cartoon;
set ray_trace_mode, 2;  # black and white cartoon
bg_color white;
set antialias, 2;
ray 600,600
png /tmp/1l9l.png

Link 

amber

AMBER 是跑分子模擬的套裝軟體，比較普及被使用的 GROMACS 套裝軟體。

以下簡單介紹基本的四個步驟：

選擇力場=>讀取原子檔案=>能量最小化=>分子模擬開始。

步驟一 選擇力場
在我們進行計算之前，要先選擇我們計算的力場(如AMBER中的ff99)每一個分子動力學中，都會有一套參數來描繪我們實際現實的狀況，就像下面這個式子，它包含了兩個原子間鍵結所產生的位勢，三個原子形的角位勢，四個原子可以形成立體角位勢，加上弱作用力的凡德瓦力位勢，以及靜電力位勢。現在還是有很多的學者在鑽研如何把參數減少(參數越多，計算會越花時間)，又可以越接近現實環境的力場。

步驟二 讀取原子檔案內容
一般來說，都會直接讀取PDB檔，什麼是PDB呢？PDB就是Protein Data Bank的縮寫，中文來說，就是蛋白質資料庫，我們把實際自然狀態的蛋白質，用核磁共振(NMR= nuclear magnetic resonance)或X射線的方法把蛋白質的原子結構與特性紀錄下來(rcsb pdb)。

步驟三.1 設好需要參數
裡面最重要，大概有幾個參數，比如是否是算最小能量參數(imin=1，說明請看步驟三.2)，長程非鍵相互作用範圍大小參數(cut)。什麼是長程非鍵相互作用範圍呢？因為我們在做模擬的時候，主要是受到鍵結力影響，在來就是長距離的作用力，如凡德瓦力、靜電力，我們再計算一個原子受力時，不可能把所有的原子對它的影響都計算過(這樣只要算一個時間間隔可能就要一個月)，所以，一般模擬會用12埃(angstrom，光譜線波長單位)來估最重要的長距離非鍵結力的影響。

步驟三.2 跑最小能量調整
在跑分子模擬前，需要先做最小能量調整，因為我們運用任何原子結構來計算時，不見得是最適當的距離，所以通常我們都會先讓原子在原先給的位置做適當的微調，就有點像上火車當充滿人擠人的過程中，我們會先找到適當與其他乘客的距離一樣。

步驟四.1 設好需要參數
接下來我們要開始正式跑分子動力模擬囉！主要是最小能量參數imin要設為0，，起始溫度(tempi)，最後溫度(temp0)，時間間隔(dt)等等。設好後可以得到我們要的輸入檔(Input file)。

步驟四.2 開始跑分子動力模擬

另外，我再補充跑分子動力模擬中，需要用到的兩個重要的資訊「拓樸檔」與「GB模組」。

「拓樸檔」
通常一個蛋白質不是存有只有原子名稱與位置，還要有每個原子重要的性質與特性，比如電荷數、重量、編號、位勢常數等等非常多資訊，所以除了「空間結構檔外，還要有「拓樸檔」。

「GB模組」
GB模組，在英文也就是GENERALIZED BORN model。我們在模擬的時候，常常會計算蛋白質在溶液裡面的分子運動，但是往往這樣計算量會很大，所以，我們會以簡化的假想位勢，來取代我們真實的水分子。GB模組，可算是很普及被使用的假想位勢。


 

Link

DSSP

DSSP - 簡介

DSSP是用於對蛋白質結構中的氨基酸殘基進行二級結構構像分類的標準化算法，由Wolfgang Kabsch和Chris Sander設計。 DSSP數據庫是由此算法生成的一個存放蛋白質二級結構分類數據的數據庫，其中包括了PDB數據庫(Protein Data Bank)中的所有條目。  

算法名稱Define Secondary Structure of Proteins 由作者在其原始論文中，作為實現該算法的Pascal語言程序名稱所提及。 DSSP數據庫英文全稱為 Definition of Secondary Structure of Proteins。

DSSP - 算法原理

DSSP算法使用PDB格式的原子級分辨率的蛋白質三維結構坐標集數據，依靠以靜電學定義進行的氫鍵識別，以及對主鍊和側鏈二面角的計算，從而得到每個氨基酸殘基的二級結構構象參數。 算法使用的鍵能公式為：

E=q1q2{1/rON+1/rCH-1/rOH-1/rCN}·332 kcal/mol
 
其中q1=0.42，q2=0.20。
 
以上參數來自分別為-0.42e和+0.20e的羰基氧原子和氨基氫原子間的局部電荷、以及羰基碳原子和氨基氮原子間的互斥電荷。 當鍵能小於-0.5 kcal/mol時，DSSP算法將其定義為一個氫鍵。

DSSP根據氫鍵模式，可以識別八種類型的二階結構構像，各構像擁有各自的標識符。  

這八種構像為：310螺旋， α-螺旋 ，π螺旋三種螺旋，分別以G、H和I標識，它們的識別特徵是殘基與主鏈中分別3、4、5個後續殘基間形成氫鍵並產生重複序列；兩種β-折疊片 中的氫鍵對類型，平行與反向平行橋中，單獨的橋結構標識為B，即β橋，含有β凸起的折疊片標識為E；在其餘構像種，包含即具有螺旋典型特徵，連接螺旋與折 疊結構的轉角區域表示為T，具有高曲率，即從第i個主鏈α碳原子指向第i+2個的向量與從第i-1個α碳原子指向第i個的向量間夾角小於70°的區域標識 為S，這也是程序中唯一不依靠氫鍵分類的構像，剩餘的環狀區域標識為L。

這八種構像在實用中又常被歸為三個大類：螺旋(G、H和I)、折疊(E和B)以及環狀結構(其餘所有標識符)

Source

ATCC 

Lists of specialized materials
Extremophiles (36.2KB/.pdf)




DSMZ

Masker 安裝

INSTALLATION:
(0) Unpack masker  >> tar -zxvf masker.tgz
(1) go to masker directory  >>  cd masker
(2) clean >> make clean
(3) set resolution of masks >> make setxxx,  xxx =  128, 256, 512, 1024, 2048 or 4096
(4) compile (make sure g95 is installed first) >> make all
(5) Install masks >> make install
(6) Run masker  >>  xpdbmask


缺 g95 ( fortran )

>> port install g95 @0.92
--->  Computing dependencies for g95
--->  Fetching archive for g95
--->  Attempting to fetch g95-0.92_2+gcc42.darwin_11.x86_64.tbz2 from http://packages.macports.org/g95
--->  Fetching g95
--->  Attempting to fetch gcc-core-4.2.4.tar.bz2 from http://jog.id.distfiles.macports.org/macports/mpdistfiles/gcc42
--->  Attempting to fetch g95_source.tgz from http://jog.id.distfiles.macports.org/macports/mpdistfiles/gcc42
--->  Verifying checksum(s) for g95
--->  Extracting g95
--->  Applying patches to g95
--->  Configuring g95
--->  Building g95
--->  Staging g95 into destroot
--->  Installing g95 @0.92_2+gcc42
--->  Activating g95 @0.92_2+gcc42
--->  Cleaning g95

python 格式

list  [ ]   為有序的資料結構，且可以動態新增或是刪除，是可改變的序列


tuple (  ) 和 List 很像，但是是個不可改變的序列, 用小括號對 (parentheses) 


dict  {  } 鍵(Key)與值(Value)的對應關係




Link  Link

Consurf server install

0. copy file and mkdir bin

1. tar -xvf Python-2.6.6.tar
    cd Python-2.6.6
    ./configure --prefix=/home/kevet/bin/pkg/ --exec-prefix=/home/kevet/bin/pkg
    make
    make install
    export PATH="/home/kevet/bin/pkg/bin:$PATH"
   
2.    Install by python2.6
    Scipy Numpy Biopython

3.    Install Muscle / blast / Rate4 / DSSP / CDhit

4.  Rename the source directory in setup file : "single_seq_rate4site.ini"

5.  Run python2.6 1_parsepdb_557_consider_whole_chains_case_heavy_bu_normal.py 1CRN_A Y

    Y: mean to calculate the conservation values

6.    Final data find in folder PDB " 1CRNA.profile_chs "

7.    Field mean: code - 1_parsepdb_557_consider_whole_chains_case_heavy_bu_normal.py - line 917




PS.  mac 安裝 rate4site

gunzip Rate4Site.tar.gz   
tar -xvf Rate4Site.tar   
g++ -o r4s.exe -O3 -ftemplate-depth-250 *.cpp

ps. Server run

kevet@c-105:~/ConSurf$ python2.6 1_parsepdb_557_consider_whole_chains_case_heavy_bu_normal.py 1BGA_A Y

1BGA_A

 !!! Residue ARG  448 A has  8 instead of expected   7 sidechain atoms.
     last sidechain atom name is  OXT
     calculated solvent accessibility includes extra atoms !!!

 !!! HEADER-card missing !!!
 !!! COMPOUND-card missing !!!
 !!! SOURCE-card missing !!!
 !!! AUTHOR-card missing !!!
447
447
447
1788
main::logevent() called too early to check prototype at single_seq_rate4site.pl line 144.
[START] --------------------------------------------------------
[START] single_seq_rate4site version 1.0 starting...
[INFO] blasting database #1. minimum clusters: 5
[INFO] blast #1: cmd line: /home/kevet/ConSurf/ncbi-blast-2.2.25+/bin/psiblast -query pdb/1BGA/A/1BGAA.fasta -out pdb/1BGA/A/seq.blast -db /temp/db/nr -num_iterations 3 -out_ascii_pssm pdb/1BGA/A/seq.pssm -evalue 0.001 -num_threads 4 -num_descriptions 4000 -num_alignments 4000 1>>/home/kevet/ConSurf/pdb/1BGA/A/single_seq_rate4site.log 2>>/home/kevet/ConSurf/pdb/1BGA/A/single_seq_rate4site.log'
[INFO] blast #1: done
[INFO] blast #1: blast round 3 saved to 'pdb/1BGA/A/seq_last_round.blast'
[INFO] blast #1: choosing homologs from blast
[INFO] blast #1: 4117 homologs chosen from blast result
[INFO] blast #1: screen homolougs by redundancy rate, according to clusters/home/kevet/ConSurf/cd-hit-v4.5.4-2011-03-07/
[INFO] blast #1: 2142 homolog clusters found by cd-hit.
[INFO] building final sequences file 'seq_final.fasta'
[INFO] running MUSCLE, command line: '/home/kevet/ConSurf/muscle3.8.31_i86linux64 -in pdb/1BGA/A/seq_final.fasta -out pdb/1BGA/A/msa.aln -clwstrict -quiet 1>>/home/kevet/ConSurf/pdb/1BGA/A/single_seq_rate4site.log 2>>/home/kevet/ConSurf/pdb/1BGA/A/single_seq_rate4site.log'
[INFO] running rate4site, command line: /home/kevet/ConSurf/rate4site.3.2.source/rate4site -s msa.aln -o r4s.res -a 1BGAA -bn 1>rate4site_stdout.txt
Optimizing branch lengths and alpha...

likelihoodComputation::getLofPos: likelihood of pos was zero!
rate4site: errorMsg.cpp:40: static void errorMsg::reportError(const std::string&, int): Assertion `0' failed.
sh: line 1:  9243 Aborted                 /home/kevet/ConSurf/rate4site.3.2.source/rate4site -s msa.aln -o r4s.res -a 1BGAA -bn >rate4site_stdout.txt
[WARNING] error found in rate4site: likelihoodComputation::getLofPos: likelihood of pos was zero!

[WARNING] rate4site output file is empty
[INFO] running rate4site SLOW version, command line: /home/kevet/ConSurf/rate4site.3.2.source/rate4site_slow -s msa.aln -o r4s.res -a 1BGAA -bn 1>rate4site_stdout.txt
Optimizing branch lengths and alpha...
Computing the rates...
[INFO] Done :)
[END]
[END]   Run Time           : 0h 54m 58s

Clean_Uniprot

What is the 'Clean_Uniprot' database?

'Clean_Uniprot' is a modified version of the UniProt database aimed to screen the more reliable sequences based on two criteria:

1. If the "Decription" (DE) field contain "Disease", "RIKEN", "variant", "mutation", "mutant" or "whole genome shotgun sequence" the sequence is removed;
2. If the database is "TrEMBL" and the "Comments" (CC) lines contain the word "CAUTION" the sequence is removed.

NCBI dataset

ftp://ftp.ncbi.nlm.nih.gov/blast/db/

Find "指令"

Linux 環境變數

lab 的環境變數是   vi ~/.profile   or  ~/.bashrc

1. ~/.profile

export PATH="$PATH:~/bin:/hom:/home/kevet/proflex/prog/first/bin/:/home/kevet/first:/home/kevet/first/lib:/home/kevet/bin/pkg/bin:"

#export GMXLIB="/pkg/x86_64/gromacs-3.3.1/share/gromacs/top"

#export FRODANHOME="/home/kevet/frodaN/bin/"

export g03root="/pkg/`uname -i`/chem/g03"
source "$g03root/g03/bsd/g03.profile"

export GAUSS_SCRDIR="/temp/$USER/g03"
source /home/tthuang/Programs/gromacs-4.0.7/bin/GMXRC.bash

export PROFLEX_HOME="/home/kevet/proflex/prog"

alias nwchem='/pkg/x86_64/chem/nwchem-4.7o/bin/nwchem'
ulimit -s unlimited

PS1='\u@\h:\w\$ '
set prompt $PS1


2. ~/.bashrc
/home/kevet/.bashrc: No such file or directory





echo $PATH

1. 直接用export命令
   export PATH="/home/kevet/bin/pkg/bin:$PATH"

2. 修改profile文件

#vi /etc/profile 
在里面加入: export PATH="/home/kevet/bin/pkg/bin:$PATH"

3. 修改.bashrc文件

# vi /root/.bashrc
在里面加入：export PATH="/home/kevet/bin/pkg/bin:$PATH"


 “/bin”、“/sbin”、“/usr/bin”、“/usr/sbin”、“/usr/local/bin”等路徑已經在系統環境變量中了，如果可執行文件在這幾個標準位置，在終端命令行輸入該軟件可執行文件的文件名和參數(如果需要參數)，回車即可。
如果不在標準位置，文件名前面需要加上完整的路徑。 不過每次都這樣跑就太麻煩了，一個“一勞永逸”的辦法是把這個路徑加入環境變量。 命令“PATH=$PATH:路徑”可以把這個路徑加入環境變量，但是退出這個命令行就失效了。

要想永久生效，需要把這行添加到環境變量文件裡。 有兩個文件可選：
“/etc/profile”和用戶主目錄下的“.bash_profile”，

“/etc/profile”對系統裡所有用戶都有效，
用戶主目錄下的“.bash_profile”只對這個用戶有效。


“PATH=$PATH:路徑1:路徑2:...:路徑n”，意思是可執行文件的路徑包括原先設定的路徑，也包括從“路徑1”到“路徑n”的所有路徑。 當用戶輸入一個一串字符並按回車後，shell會依次在這些路徑裡找對應的可執行文件並交給系統核心執行。 那個“$PATH”表示原先設定的路徑仍然有效，注意不要漏掉。 某些軟件可能還有“PATH”以外類型的環境變量需要添加，但方法與此相同，並且也需要注意“$”。

注意，與DOS/Window不同，UNIX類系統環境變量中路徑名用冒號分隔，不是分號。 另外，軟件越裝越多，環境變量越添越多，為了避免造成混亂，建議所有語句都添加在文件結尾，按軟件的安裝順序添加。

格式如下：

# 軟件名-版本號
PATH=$PATH:路徑1:路徑2:...:路徑n

其他環境變量=$其他環境變量:...

在“profile”和“.bash_profile”中，“#”是註釋符號，寫在這裡除了視覺分隔外沒有任何效果。
設置完畢，註銷並重新登錄，設置就生效了。 如果不註銷，直接在shell裡執行這些語句，也能生效，但是作用範圍只限於執行了這些語句的shell。
 
相關的環境變量生效後，就不必老跑到軟件的可執行文件目錄裡去操作了。

轉貼 BLAST E value

一般的，当我们使用BLAST（是一种用于在数据库当寻找任何蛋白质或者基因序列与你的目标序列一致的程序）时，我们会注意到这里有一个E值。

那么这个E value是什么呢？怎么来理解这个值呢？

下面是一个平常的blast结果，

Sequences producing significant alignments:	Score (S)	E
gi|83574104|Moth_2374|sporulation – prote…	202	2e-53
gi|83573446|Moth_1696|Sporulation – prote…	112	1e-26
gi|83571874|Moth_0087|sporulation – prote…	95	3e-21
gi|83573435|Moth_1685|Substrate-binding -…	27	1.0

后面有两个值，一个是S值，一个E值。可以发现，结果是依据S值的高低来显示的。
S值表示两序列的同源性，分值越高表明它们之间相似的程度越大。

E值就是S值可靠性的评价。它表明在随机的情况下，其它序列与目标序列相似度要大于这条显示的序列的可能性。所以它的分值越低越好。

E值的计算：   E = Kmn (e^-lambda*S)

其中，K和lambda与数据库和算法有关，是个常量；m代表目标序列的长度，n代表数据库的大小，S就是前面提到的S值。

通常来讲，我们认为E值小于10^-5就是比较可性的S值结果。我们可以想象，相同的数据库，E=0.001时如果有1000条都有机会S值比现在这个要高的话，那么不E设置为10^-6时可能就会只得到一条结果，就是S值最可靠的那个。

但是E值也不是万能的。它在以下几个情况下有局限性：

1. 当目标序列过小时，E值会偏大，因为无法得到较高的S值。
2. 当两序列同源性虽然高，但有较大的gap（空隙）时，S值会下降。
   这个时候gap scores就非常有用。
3. 有些序列的非功能区有较低的随机性时，可能会造成两序列较高的同源性。

BLAST试图去避免这些问题，但是还是应该自己有个清晰的概念。
E值总结：

E值适合于有一定长度，而且复杂度不能太低的序列。
当E值小于10^-5时，表明两序列有较高的同源性，而不是因为计算错误。
当E值小于10^-6时，表时两序列的同源性非常高，几乎没有必要再做确认。










 1、基本概念
   相似性(Similarity)
是指序列比对过程中用来描述检测序列和目标序列之间相同或相似碱基或氨基酸残基占全部比对碱基或氨基酸残基的比例的高低，属于量的判断。
   同源性(Homology)
是指从某一共同祖先经趋异进化而形成的不同序列。只有当两个蛋白质在进化关系上具有共同的祖先时，才可称它们为同源的，属于质的判断。

相似性和同源性的关系

当相似程度高于50%时，比较容易推测检测序列和目标序列可能是同源序列；
而当相似性程度低于20%时，就难以确定或者根本无法确定其是否具有同源性。
总之不能把相似性和同源性混为一谈。所谓“具有50%同源性”，或“这些序列高度同源”等说法，都是不确切的，应避免使用。

序列相似性比较和同源性分析

序列相似性分析：就是用来计算待研究序列与某序列之间的相似性程度，常用的软件包有BLAST、FASTA等；
序列同源性分析：是将待研究与来自不同物种的序列中进行进化分析，以确定该序列与其它序列间的亲源关系。常用的程序包有Phylip及Mega等进化分析软件；
全局比对与局部比对
   全局比对
寻找序列在全长范围内最佳比对。 
常用算法如：Needleman-Wunsch algorithm(Needle)
在线程序如： Needle
      局部比对
寻找序列在局部区域的最高比对打分。
常用算法如：Smith-Waterman algorithm， blast,fasta等
在线程序如： Water
   Needle及Water的在线程序
http://bioweb2.pasteur.fr/alignment/intro-en.html
也可以本地安装Emboss执行以上程序

局部相似性比对的生物学基础
蛋白质功能位点往往是由较短的序列片段组成的，尽管在序列的其它部位可能有插入、删除等突变，但这些关键的功能部位的序列往往具有相当大的保守性。
而局部比对往往比整体比对对这些功能区段具有更高的灵敏度，因此其结果更具生物学意义。
BLAST程序常用的两个评价指标
Score：
使用打分矩阵对匹配的片段进行打分，这是对各对氨基酸残基（或碱基）打分求和的结果，一般来说，匹配片段越长、相似性越高则Score值越大，结果越可信。
E-value：
BLAST程序在搜索空间中可随机找到获得这样高分的序列的可能性(期望值)，因此E-value越高，则代表结果越有可能是随机获得的，也就越不可信。搜寻空间大小约略等于查询序列的长度乘以全部database序列长度的总和，再乘以一些系数。

我们在获得一个Blast结果时需要看这两个指标。
如果Blast获得的目标序列的Score值越高并且E-value越低表明结果越可信，反之越不可信。
其它的一些重要关键概念
HSP(High Scoring Pair)：
在局部比对时，得分高的匹配序列被称为高分值片段。
LCRs(low compositional complexity regions)：
低复杂度区域，即这些区域的组成有某些偏好，比如DNA中的简单重复序列。在蛋白质中一些残基过多表现。在进行BLAST比较时，将会把LCRs屏蔽掉，防止它们过高评价匹配的显著性。在核酸中用n，在蛋白质中用X代替。
gi(GenBank Index)
特定于GenBank数据库中所赋予每一条序列的特定索引数字。
nr(non-redundant database)
非冗余数据库，该库信息多，并且无冗余序列
2、常用BLAST程序
BLAST(Basic Local Alignment Search Tool) 基于匹配短序列片段，并用一种强有力的统计模型来确定未知序列与数据库序列的最佳局部联配的一种程序。
主要的BLAST程序
PSI-BLAST(位置相关的迭代 BLAST)
这个程序主要用来搜索蛋白质的“远亲”。
首先，用户提交的蛋白质序列的所有“近亲”的列表被建立起来，然后这些蛋白质被结合成一种平均的“特征序列” 。
再用这个特征序列在蛋白质数据库中进行搜索，就会找出更大的一组蛋白质的列表。再将这个蛋白质列表生成一个不同的特征序列，这个序列被用来迭代地运行上述过程。
通过在搜索中包含相关的蛋白质，PSI-BLAST对于寻找已知蛋白进化上的“远亲”的灵敏度要比一般的blastp高很多。
其它的一些BLAST子程序
Gapped BLAST
允许在它产生的比对（alignments）中存在缺口。
Megablast
该程序使用“模糊算法”加快了比较速度，可以用于快速比较两个较长的序列。
discontiguous megablast
与megablast不同的是主要用来比较来自不同物种之间的相似性较低的分歧序列。
PHI-BLAST
模式发现迭代BLAST。
Bl2seq
给定两个序列，相互进行BLAST比对，快速检查两个序列是否存在相似性片断
Specialized Blast Specialized BLAST pages
CD - Search
是使用RPS - BLAST程序以一个蛋白质序列与保守结构域数据库(Conserved Domain Database) 做比较。
Pairwise BLAST
Pairwise BLAST是用BLAST程序实现两个序列之间的比较。
IgBLAST
IgBLAST被开发出来以便於分析在GenBank中的免疫球蛋白的序列。它允许用blastp或blastn来搜索nr资料库或一个由免疫球蛋白生殖系变化区基因的特殊的资料库。
3、BLAST算法简介
BLAST 是一种基本局域联配搜索工具，主要用来搜索数据库中相似序列。
它的搜索速度快并且把数据库搜索建立在了严格的统计学基础之上，是目前最常用的同源检索工具，是由Altschul SF et al(1990)提出的一种算法。
BLAST的基本步骤
将待检索序列分割成长度为w的连续子串
快速找出数据库中所有与固定长度w完全配对的位置
以此位置为起点进行延伸比对，并计算出最高分数
将最高分标准化，并按此分数进行排序
换算成期望值(E-VALUE)
显示出符合Score及E-value的序列
4、BLAST常用参数设置
在NCBI进行BLAST的操作程序非常简单，只要将你的序列贴进去，点几下鼠标就会得到结果，但是如果能正确的修改一下BLAST的参数，可能你会得到更好的结果！以下我们一起讨论一下如何来修改BLAST的参数！
BLAST的具体过程：
登陆NCBI的BLAST主页http://www.ncbi.nlm.nih.gov/BLAST/
根据序列类型及目的选择合适的程序
填写表单信息
提交任务
查看和分析结果
BLAST程序的选择
在BLAST程序选择上，应尽可能地利用blastp从蛋白质水平进行检索，然后用blastx、tblastn、tblestx从DNA或蛋白质翻译水平进行检索，最后才用blastn进行DNA水平进行检索。
当然如果为非编码序列只有采用blastn进行检索。
E-value的设置
如果检索的序列较短，可适当的提高E值,否则可能会找不到目的序列，反之如果序列较长可适当提高E值。
通常无论是从DNA水平，还是蛋白质水平进行检索，E值设为1通常可满足要求。
Word size的选择
BLAST算法将查询序列分割成一系列具有字段长度的小的序列段进行数据库搜索，因此当此值越小得到的搜索结果越多，但假阳性也越多，服务器负担也越重。
对于蛋白质搜索,窗口大小可设置为3或2,默认为3；对于核酸搜索，默认的字段长度是11，可选择7，11或15。
因此如果你对搜索的结果不满意时可以试着降低Word size的值。
打分矩阵的选择
比对所选用的记分矩阵对最终结果影响也很大。
一般高值BLOSUM矩阵和低值PAM矩阵最适合于研究近相关的蛋白质序列。低值BLOSUM矩阵和高值PAM矩阵最适合于研究远相关的蛋白质序列。
一般情况BLOSUM62检测各种蛋白的效果比BLOSUM60和BLOSUM70稍好，比PAM矩阵好得多。
在BLAST五个程序中只有BLASTN不需要这些矩阵，搜索时不必选定。
空位罚分的选择
严紧的罚分很难将本来很相似的序列对准；而松弛的罚分甚至可以使两个无关的序列达到100%的相似性。
一般情况下程序默认的空位罚分(11/1)基本能满足检索要求，但对具体的查询序列，采用不同的空位罚分方法会取得不同的检索效果。
低复杂区域及重复区域的处理
无论是DNA序列类似性检索，还是蛋白质序列类似性检索，一般都应该去除查询序列中的低复杂区域。
就蛋白质序列检索而言，不必去除序列中的重复片段，但对DNA序列检索，就必须去除序列中的重复片段。
5、本地BLAST的安装
大家一般都做过基于网络的BLAST ，但网络BLAST一般只能搜索一个序列，要搜索多个序列，特别是做大量的数据比较时，网络BLAST几乎是不可能的，这个时候我们就可以考虑做本地BLAST了。
使用本地BLAST的原因
特殊的数据库要求
涉及序列的隐私与价值
批量处理
与其它本地程序配合使用
其他原因？？
本地BLAST构建步骤
下载BLAST的安装程序
将BLAST保存到适当的位置
点击安装程序来安装BLAST
设置BLAST的相应参数
下载BLAST的本地安装程序
可以到NCBI的官方网站下载最新的BLAST程序。
下载网址：ftp://ftp.ncbi.nih.gov/blast/executable
注意一定要选择和你的计算机操作系统相匹配的程序，如Windows系统要下载“blast-2.2.18-ia32-win32.exe”。由于大部分螺友都在用Windows系统我们下面的讲座也是在Windows系统下进行
本地BLAST的相应参数设置
告诉BLAST程序你的数据及数据库放在哪
1.建立一个新的文件并命名为：ncbi.ini
2.在该文件中输入四行数据如下所示：
[NCBI]
Data=“C:\ncbi-blast\data”　（你的数据存放的文件夹）
[BLAST] BLASTDB=“C:\ncbi-blast\db” （数据库存放在的文件夹）
3.将该文件拷贝到你的Windows或Winnt目录里
路径设置步骤
右键点击我的电脑
选择属性
再选择系统属性
选择高级标签
选择环境变量
双击path
在路径中填入你的BLAST的可执行文件所在目录
有的时候还需要重新启动电脑
6、本地BLAST的使用
数据库的获取
最简单的方法是直接到NCBI或别的网站去下载
也可以将自己的序列，或与自己工作相关的序列进行整理构建成一个小型的数据库
注意：以上文件格式一般可存为fasta格式
构建BLAST用的数据库
将已构建好的数据拷贝到你所设定的数据库所在文件夹
运行cmd命令
在cmd环境中输入如下所示命令
formatdb –i inseqs.fa –p F –o T –n db_name
命令结束后你会发现在你的数据库文件夹里多了一些以db_name开头的文件，这些就是BLAST所需要的一些文件
输入过程
Formatdb的一些参数说明
-i 输入文件，只能是一个文件
-o Parse options （默认是F) T - True: Parse SeqId and create indexes. F - False: Do not parse SeqId. Do not create indexes
-p 文件类型 （默认是T) T - protein F - nucleotide [T/F] Optional
-n 数据库名称 不指出的话默认为输入文件名
更多选项请参阅 解压后的doc文件夹的formatdb.html文件
进行BLAST搜索
在命今行下录入blastall 命令及相应的参数
打开输出文件分析结果，如果结果不好可以试着调整参数再次进行BLAST
如下所示命令：
blastall -p blastn -d db_name -i QUERY -o out.QUERY
Blastall的一些参数说明(1)
-p 程序名 包括
blastp: 通过蛋白质序列搜索蛋白质序列数据库
blastn: 通过核酸序列搜索核酸数据库
blastx: 通过翻译后的核酸序列搜索蛋白质数据库
tblastn: 通过蛋白质序列搜索翻译后的核酸数据库
tblastx: 通过翻译后的核酸序列搜索翻译后的核酸数据库
-d 数据库名称 与formatdb中-n选项一致
-i 输入文件 不指明的话默认为STDIN
-o 输出文件 不指明的话默认为STDOUT
Blastall的一些参数说明(2)
-e:设置e-value
-m:比对结果显示格式选项，缺省值为0 ，即pairwise格式。另外还可以根据不同的需要选择1~6等不同的格式。
-I:在描述行中显示gi号[T/F]，缺省值F
-b:显示的比对结果的最大数目，缺省值250
-F:对于要查询的序列做低复杂度区域(low complexity regions, LCR)的过滤[T/F]，缺省值T。对blastn用的是DUST程序，其他比对用的是SEG程序。
Blastall的一些参数说明(3)
-G: 打开一个gap的罚分（0表示使用缺省设置值），默认0
-E: 扩展一个gap的罚分（0表示使用缺省设置值），默认0
-q: 一个核酸碱基的错配(mismatch)的罚分（只对blastn有效），缺省值-3
-r : 一个核酸碱基的正确匹配(match)的奖分（只对blastn有效），缺省值1
-M: 所使用的打分矩阵，缺省值BLOSUM6244
Bl2seq(两条序列的BLAST)
bl2seq的绝大部分参数是与通用检索程序blastall一致的，只是没有了-d 的选项。另外增加了两个输入选项：
-i：第一个输入序列文件
-j：第二个输入序列文件
注：这两个输入序列都应该是FASTA格式，各自的序列类型--核酸或蛋白--应由所选择的-p 参数决定
命令如下所示：
bl2seq -i query.fa -j sbjct.fa -e 0.01 -o out
Psi-BLAST
Psi-BLAST是由blastpgp命令实现的，它的大部分参数是与blastall一致的，只有少数与迭代检索相关的选项是特别的：
-j: 最大迭代检索的次数，缺省值1，即等同与在blastall中所使用blastp程序
-h: 在每轮检索后构建新的打分矩阵时所选择的序列的期望值(E value)的阈值，缺省值0.001
-C: 将生成的位点特异性的打分矩阵输出到一个文件（二进制格式）
-R: 从文件读取一个原先输出的位点特异性的打分矩阵，然后使用这个矩阵来继续进行以后的检索比对
-Q: 输出一个可读的文本（ASCII）格式的PSI-BLAST的打分矩阵
-B: 设置让blastpgp读取一个已经存在的多重比对文件来构建位点特异性的打分矩阵而进行以后的检索
如下命令所示：
blastpgp -i query.faa -d db_name -o query_out
Fastacmd从数据库中提取序列
-a：是否提取重复accession号的序列[T/F]
-l ：设置输出的序列文件每行的字符数
-t ：设置在FASTA格式的序列的描述行中只包含gi号[T/F]
-o：输出文件名
命令如下：

fastacmd -d db_name -s p38398



From： http://www.biostatistician.cn/thread-467-1-1.html



Using BLASTClust to Make Non-redundant Sequence Sets

BLASTClust is a program within the standalone BLAST package used to cluster either protein or nucleotide sequences. The program begins with pairwise matches and places a sequence in a cluster if the sequence matches at least one sequence already in the cluster. In the case of proteins, the blastp algorithm is used to compute the pairwise matches; in the case of nucleotide sequences, the Megablast algorithm is used.

In the simplest case, BLASTClust takes as input a file containing catenated FASTA-format sequences, each with a unique identifier at the start of the definition line. BLASTClust formats the input sequence to produce a temporary BLAST database, performs the clustering, and removes the database at completion. Hence, there is no need to run formatdb in advance to use BLASTClust. The output of BLASTClust consists of a file, one cluster to a line, of sequence identifiers separated by spaces. The clusters are sorted from the largest cluster to the smallest.

BLASTClust accepts a number of parameters that can be used to control the stringency of clustering including thresholds for score density, percent identity, and alignment length. The BLASTClust program has a number of applications, the simplest of which is to create a non-redundant set of sequences from a source database. As an example, one might have a library of a few thousand short nucleotide sequence reads and wish to replace these with a non-redundant set. To produce the non-redundant set, one might use:

blastclust -i infile -o outfile -p F -L .9 -b T -S 95

The sequences in "infile" will be clustered and the results will be written to "outfile". The input sequences are identified as nucleotide (-p F); "-p T", or protein, is the default. To register a pairwise match two sequences will need to be 95% identical (-S 95) over an area covering 90% of the length (-L .9) of each sequence (-b T) . Using "-b F" instead of "-b T" would enforce the alignment length threshold on only one member of a sequence pair. The parameter "S", used here to specify the percent identity, can also be used to specify, instead, a "score density." The latter is equivalent to the BLAST score divided by the alignment length. If "S" is given as a number between 0 and 3, it is interpreted as a score density threshold; otherwise it is interpreted as a percent identity threshold.

To create a stringent non-redundant protein sequence set, use the following command line:

blastclust -i infile -o outfile -p T -L 1 -b T -S 100
In this case, only sequences which are identical will be clustered together. The “blastclust.txt” file in the standalone BLAST package details the full range of BLASTClust parameters.

Linux 常用指令

rm -fr * 刪除目前目錄下所有檔案及子目錄(-r)，不經確認(-f)
rm -fr back 一次刪除目錄back及其下所有檔案和子目錄

mv  移動檔案或更改檔名

cp -vr public_html/* test 將public_html下所有檔案及子目錄複製到test中(-r)，並顯示出複製的過程(-v)

ping 163.17.51.1 -c 5 ping 163.17.51.1 五次

tar -czvf test.tar.gz    * 將目前目錄所有檔案及子目錄打包成一個檔案test.tar.gz 並呼叫gzip壓縮指令壓縮
tar -xzvf junior_0727.tar.gz 將junior_0727.tar.gz解開到目前目錄下

( c建一個新的tar檔，x解開tar檔，z呼叫gzip，v顯示製作過程，f指定檔案名稱)