EMBOSS

序列提取程序

• extractseq

   Example1:Extract the regions 30 to 45
   extractseq sequence -regions "30-45"
   Example2: extract the regions 1787-1912, 782-856;
   extractseq sequence -reg "1787..1912,782..856"
   EXAMPLE3: extract the regions 782-856, 1787-1912 all to separate output sequences:
   extractseq tembl:x65921 -reg "782..856,951..1095,1557..1612,1787..1912" stdout
 -separate

   其他高級(jí)檢索：
 -snucleotideq :如果序列是核苷酸
  -sprotein1 :如果序列是蛋白質(zhì)
  -slwer1 :make lower case
  -sid1 <entryname>
  -squery1 <query fields or ID list>
  -osformat1 <輸出序列格式>
  -ossingle2 :seperate file for each entry
  -ufo1 <UFO featrues>

雙序列比對(duì)

• needle
• water
• stretcher: 改進(jìn)的Needleman-Wunsch動(dòng)態(tài)規(guī)劃算法局相似性雙序列比對(duì)程序，占用內(nèi)存較少，運(yùn)行時(shí)間較長(zhǎng)。
  • 例：比較af164138時(shí)，needle顯示“Died: Sequences too big. Try 'stretcher'”

-datafile <矩陣>：選擇blosum做蛋白比對(duì)，或者EDNAFULL做核酸序列比對(duì)（事實(shí)上可自動(dòng)識(shí)別，不需再手動(dòng)設(shè)置）
-sbegin1 <數(shù)字>：作為序列起始
-send1 <integer>： 作為序列最后
-sprotein: 序列是否是蛋白，不是則可繼續(xù)讀取
-snucleotide: 序列是否是核苷酸，同上
-sreverse：默認(rèn)“-”前為mRNA，將其反轉(zhuǎn)成DNA后再比對(duì)（“-sreverse1”，即使在最后也將翻轉(zhuǎn)1）
-option：菜單模式

• matcher：局部相似性雙序列比對(duì)，通過設(shè)定參數(shù)可同時(shí)輸出多個(gè)相似性片段。
  • a rigorous algorithm based on Bill Pearson's lalign application
  • 用于相似度較低的矩陣（與water相似）

EXAMPLE:找到十個(gè)最佳a(bǔ)lignment；
matcher <seq1> <seq2> -alt 10;
額外的參數(shù)：
-alternatives <integer>（-alt）: 給出局部高分匹配序列，默認(rèn)1只給出得分最高的序列；但是在cDNA多域蛋白與基因組DNA比較中，需要修改可能會(huì)得到其他有趣且重要的片段（若序列過短則會(huì)給出全局比對(duì)）；（得分越高說明序列的相似度越高）
-gapextend（-gape）: 【所有序列默認(rèn)都是4】一般認(rèn)為取幾個(gè)長(zhǎng)的空位比去很多短空位要合理，所以gape的值一般較小（除去某些特殊情況，比如單端測(cè)序使得序列有誤時(shí)傾向于選擇多個(gè)短空位，可以通過調(diào)低gapo實(shí)現(xiàn)（罰分針對(duì)的是第一條的空位插入情況））；

• supermatcher：基于局部相似性的雙序列快速比對(duì)，適用于超長(zhǎng)序列之間或序列和數(shù)據(jù)庫(kù)條目之間相似性比對(duì)，所得結(jié)果為近似解而非最優(yōu)解

EXAMPLE:
supermatcher @eclac.list tembl:j016136 -word 50
可選參數(shù)：
-minscore <float>: 輸出的匹配最小得分
-width <integer>: alignment的寬度
-wordlen <integer>: 讀取步長(zhǎng)

• seqmaterall：基于局部相似性的雙序列快速比對(duì)，用于尋找一組序列中所有匹配字串（無文件）
• esim4：將mRNA序列定位于基因組序列：

esim4 <seq1> <seq2> :
seq1 = , 577 bp
seq2 =  ((no header)), 846 bp

1-132  (4-135)   100% ->
133-337  (253-457)   100% ->
338-577  (607-846)   100%

esim4 <mRNA> <genome>
可選參數(shù)：
-word <integer>: blast的word size參數(shù)設(shè)置
-extend <integer>: 設(shè)置 在3~10之間
-format <integer>: 0：only exon endpoints;
                   5:CDS(只顯示基因組上與mRNA相同的開始      和結(jié)束位置);
                   1：顯示序列對(duì)應(yīng)情況
-cutoff <integer>: integer 3~10            

• est2genome：將EST序列定位于基因組序列

序列變換：

• revseq:將輸入序列轉(zhuǎn)換成反向互補(bǔ)序列

EXAMPLE1:給出seq1的反義互補(bǔ)序列的sev文件
revseq seq1 seq1.sev
EXAMPLE2: 只輸出seq1的互補(bǔ)序列sev文件
revseq seq1 seq1.sev -norev
EXAMPLE3: 輸出seq1自身的反義序列
reseq seq1 seq1.sev -nocomp

-其他命令：
-notag: 輸出文去除標(biāo)題

• msbar：對(duì)序列進(jìn)行模擬突變
  • the number, size and type of mutation may be specified
  • 堿基替換是隨機(jī)單位點(diǎn)的堿基替換

交互式，主要命令都會(huì)顯示并解釋
高級(jí)參數(shù)：
-othersequence <seq>: 輸出的突變序列不會(huì)和othersequece的序列相同

• shuffleseq：對(duì)輸入序列進(jìn)行變換，產(chǎn)生隨機(jī)新序列（只改變堿基或氨基酸殘基的順序）

e.g.給出輸入序列的兩個(gè)隨機(jī)拷貝：
shuffleseq -shuffle 2
交互式命令

核酸序列CpG島分析

• cpgplot: 預(yù)測(cè)核酸序列中的CpG島，用圖形方式輸出結(jié)果

交互式命令：
-window <integer>: CG百分?jǐn)?shù)和觀察到的CG頻率以此參數(shù)設(shè)定的窗口大小來計(jì)算，并且窗口在序列上移動(dòng)，并將數(shù)值累加。
-minlen <integer>: CpG島的最小長(zhǎng)度
-minoe <float>: 設(shè)置了一組10個(gè)窗口中，(C+G)觀察值與CpG的期望值最小比率的平均值
-minpc <float>: 設(shè)置G和C在一組十個(gè)窗口中的最小平均百分?jǐn)?shù)（期望值）
-graph <格式>：ps,hpgl,meta,cps（彩色）,x11,tek,none,data,xterm,png,gif,pdf,svg

手動(dòng)添加：
-（no）plot: 作為開關(guān)，可以繪出得分

• cpgreport: 識(shí)別核酸序列中富含CpG雙核苷酸區(qū)域
  • 可人為改變CpG島輸出的閾值
  • 由于該條命令算法是只要遇到CG就會(huì)有一個(gè)得分，并且將序列上的所有得分相加，因此會(huì)發(fā)生過度預(yù)測(cè)的情況，但是可以發(fā)現(xiàn)主外顯子附近的較小的CpG島
  • 會(huì)出現(xiàn)的錯(cuò)誤
    [圖片上傳失敗...(image-1c222c-1588692442306)]

e.g.報(bào)告seq1中CpG富集區(qū)域,輸出閾值為28
cpgreport seq1 -score 28 -outflie <> -outfeat <>
(-outfeat 是特征值輸出格式；或者可直接使用交互命令)

• newcpgreport：識(shí)別核酸序列中富含CpG雙核苷酸區(qū)域新方法
  • 輸出顯示：位置、長(zhǎng)度、C+G的總量，CG百分比和期待值百分率

相較于cpgplot的特殊之處：
-shift <integer>: 改變每次位移長(zhǎng)度，當(dāng)integer=1時(shí)，與cpgplot選擇的cpg島的片段是相同的（可見二者得分算法應(yīng)該是相同的）

• newcpgseek: 識(shí)別核酸序列中富含CpG雙核苷酸區(qū)域新方法，靈敏度高。
  • 與cpgreport的算法基本相同，但是會(huì)自動(dòng)忽略把CG直接當(dāng)成一個(gè)得分17的CpG島的情況

-score <integer>: CpG被確認(rèn)的得分閾值

讀碼框分析程序

• plotorf:根據(jù)起始密碼子和終止密碼子位置用圖形方式顯示DNA序列開放讀碼框
• showorf：按一定格式顯示DNA序列及其翻譯所得蛋白質(zhì)序列
• getorf：從DNA序列中提取開放讀碼框序列，或其編碼的氨基酸序列

-find <number>:0，終止密碼子之間的翻譯；1，起始密碼子和終止密碼子之間的翻譯；2，終止子之間的核苷酸序列；3，啟動(dòng)子和終止子之間的核苷酸序列；4，啟動(dòng)子旁側(cè)序列；5，最初終止子的旁側(cè)序列；6，最后終止子的旁側(cè)序列；

• sixpack：顯示DNA序列6個(gè)開放讀碼框和翻譯所得氨基酸序列