multiple whole genome alignment

寫在前面：在上億年的進化歷程中，基因組經(jīng)歷了大大小小的改變。從小的核苷酸突變、插入、缺失到大的基因缺失、重復(fù)、基因組重排和水平基因轉(zhuǎn)移?；蚪M比對可以通過比對序列中的同源區(qū)域，找到DNA中的改變。理解每種改變的速率和模式將有助于了解多種多樣的生物學(xué)過程。

經(jīng)過google，高分推薦的multiple whole genome alignment軟件有Mavue, Mugsy和Multiz三個。這三個軟件都可以通過conda安裝。
但相比于conda，源碼附帶的readme是我們了解這個軟件的重要參考，因此我更偏愛源碼。

Mugsy: fast multiple alignment of closely related whole genomes

補充：
過了兩個月，發(fā)現(xiàn)mauve和mugsy比對都有一些小問題。一些blast可以找出的homologous gene ，在比對結(jié)果中卻找不到···。我打算換方法，先Lastz比對，然后Multiz合并。

1. Mugsy

Mugsy是Angiuoli SV和Salzberg SL于2010年，開發(fā)的一款用于multiple whole genome alignment 的工具。

特點：

  • 可以接受組裝基因組草圖文件（contig/scaffold，一個文件中有多條fasta序列）
  • 不需要參考基因組
  • 非常適合closely related genomes
  • 可以識別到duplication、rearrangement和大規(guī)模的gain/loss。
  • 目前測試的是幾十個細菌基因組。

1.1 安裝

wget https://sourceforge.net/projects/mugsy/files/latest/download
tar xvzf download

#將Mugsy安裝路徑加到PATH里
vim ~/.bashrc
export PATH="/picb/evolgen/users/gushanshan/software/mugsy/mugsy_1.2.2:$PATH"
source ~/.bashrc

1.2. 運行

1.2.1 輸入

文件要求

  • 一個或多個FASTA文件
  • 每個文件應(yīng)該包括單個物種的所有序列
  • FASTA的header中不能包括：或-
  • 模糊字符將被轉(zhuǎn)換成N

選項

-p: 輸出文件前綴
--directory：用于儲存輸出和臨時文件的目錄，必須是全路徑，否則會報錯

1.2.2 輸出

MAF格式

1.2.3 例子

mugsy腳本可以完成multiple alignment 的所有過程。
其中，核心程序是：

mugsyWGA - whole genome aligner based on Seqan::TCoffee

synchain-mugsy - segmentation program to produce locally collinear
blocks LCBs from a set of anchors

nucmer - 3.20 release bundled for convenience with new utility
delta2maf and modified delta-filter to add support for reporting
duplications

例子：

mugsy --directory /local/scratch --prefix mygenomes genome1.fsa genome2.fsa genome3.fsa

1.3 預(yù)實驗

選擇10個Lp genome，進行預(yù)實驗，記錄時間。

10 genomes
Starting Nucmer: Wed Sep 23 12:43:43 CST 2020
.........
Finished Nucmer Wed Sep 23 12:51:10 CST 2020
Starting MUGSYWGA: Wed Sep 23 12:51:10 CST 2020

Finished MUGSYWGA: Wed Sep 23 15:55:35 CST 2020
Final output (MAF format): path/output/pre_experiment.maf
Finished Wed Sep 23 15:55:37 CST 2020

話說好久啊，10個基因組花了3個小時。那我有500多個···

1.4 maf格式轉(zhuǎn)為fasta

Mugsy安裝目錄下有轉(zhuǎn)換文件：

./maf2fasta.pl < a.maf > a.fasta

MAF格式：

##maf version=1 scoring=tba.v8 
# tba.v8 (((human chimp) baboon) (mouse rat)) 
                   
a score=23262.0     
s hg18.chr7    27578828 38 + 158545518 AAA-GGGAATGTTAACCAAATGA---ATTGTCTCTTACGGTG
s panTro1.chr6 28741140 38 + 161576975 AAA-GGGAATGTTAACCAAATGA---ATTGTCTCTTACGGTG
s baboon         116834 38 +   4622798 AAA-GGGAATGTTAACCAAATGA---GTTGTCTCTTATGGTG
s mm4.chr6     53215344 38 + 151104725 -AATGGGAATGTTAAGCAAACGA---ATTGTCTCTCAGTGTG
s rn3.chr4     81344243 40 + 187371129 -AA-GGGGATGCTAAGCCAATGAGTTGTTGTCTCTCAATGTG
                   
a score=5062.0                    
s hg18.chr7    27699739 6 + 158545518 TAAAGA
s panTro1.chr6 28862317 6 + 161576975 TAAAGA
s baboon         241163 6 +   4622798 TAAAGA 
s mm4.chr6     53303881 6 + 151104725 TAAAGA
s rn3.chr4     81444246 6 + 187371129 taagga

a score=6636.0
s hg18.chr7    27707221 13 + 158545518 gcagctgaaaaca
s panTro1.chr6 28869787 13 + 161576975 gcagctgaaaaca
s baboon         249182 13 +   4622798 gcagctgaaaaca
s mm4.chr6     53310102 13 + 151104725 ACAGCTGAAAATA

FASTA格式：

>hg18.chr7 27578828 38 + 158545518
AAA-GGGAATGTTAACCAAATGA---ATTGTCTCTTACGGTG
>panTro1.chr6 28741140 38 + 161576975
AAA-GGGAATGTTAACCAAATGA---ATTGTCTCTTACGGTG
>baboon 116834 38 + 4622798
AAA-GGGAATGTTAACCAAATGA---GTTGTCTCTTATGGTG
>mm4.chr6 53215344 38 + 151104725
-AATGGGAATGTTAAGCAAACGA---ATTGTCTCTCAGTGTG
>rn3.chr4 81344243 40 + 187371129
-AA-GGGGATGCTAAGCCAATGAGTTGTTGTCTCTCAATGTG
=
>hg18.chr7 27699739 6 + 158545518
TAAAGA
>panTro1.chr6 28862317 6 + 161576975
TAAAGA
>baboon 241163 6 + 4622798
TAAAGA
>mm4.chr6 53303881 6 + 151104725
TAAAGA
>rn3.chr4 81444246 6 + 187371129
taagga
=
>hg18.chr7 27707221 13 + 158545518
gcagctgaaaaca
>panTro1.chr6 28869787 13 + 161576975
gcagctgaaaaca
>baboon 249182 13 + 4622798
gcagctgaaaaca
>mm4.chr6 53310102 13 + 151104725
ACAGCTGAAAATA
=

把比對信息抽出來，抽成一塊一塊的，不再是原來摞在一起的：

Args <- commandArgs()
a<-file(Args[6],"r")
output<-Args[7]
lines=readLines(a,n=1)
filename="1"
while(length(lines)!=0){
  if(lines !="="){
    cat(paste0(lines,"\n"),file=paste0(paste(output,filename,sep = "/"),".fasta"),append = T)
    lines=readLines(a,n=1)
  }else{
    filename=as.character(as.numeric(filename)+1)
    lines=readLines(a,n=1)
  }
}
close(a)

第一個參數(shù)是fasta文件，第二個是單一文件的輸出目錄

mkdir a #存儲輸出文件
Rscript ./singleFasta.R ./a.fasta a/

然后傳到MEGAX里就可以看了。

2. Mauve

Mauve分為原始Mauve算法和progressive Mauve算法。

original Mauve 的劣勢：

  • 比較適合closely-related 物種
  • 不比對部分基因組共享的大區(qū)域
  • 識別不到部分基因組共享的重排區(qū)域
  • 為了準(zhǔn)確估計基因組重排，必須人工調(diào)整最小LCB權(quán)重
  • 很難比對經(jīng)常出現(xiàn)片段重復(fù)的基因組

progressiveMauve算法的優(yōu)勢：

  • 能夠比對更多的基因組
  • 能夠比對分歧更大的基因組，核苷酸相似性可以低至50%
  • 不需要人工調(diào)整比對打分參數(shù)
  • 可以比對pan-genome
  • 準(zhǔn)確度更高
    progressiveMauve算法的劣勢：
  • 很慢
  • 很耗內(nèi)存

2.1 安裝

#下載地址：http://darlinglab.org/mauve/download.html
wget http://darlinglab.org/mauve/snapshots/2015/2015-02-13/linux-x64/mauve_linux_snapshot_2015-02-13.tar.gz
tar  -zxvf  mauve_linux_snapshot_2015-02-13.tar.gz 
rm mauve_linux_snapshot_2015-02-13.tar.gz
mv mauve_snapshot_2015-02-13/ mauve_2015_02_13

#把路徑加到環(huán)境變量里。因為要用命令行程序，所以路徑指定到安裝路徑下的linux-x64
export PATH="xxxxxx/mauve/linux-x64:$PATH"
source ~/.bashrc

2.2 運行

2.2.1 輸入

文件格式
輸入文件可以是FastA，Multi-FastA或 GenBank。如果一個文件中有多個序列，即Multi-FastA，Mauve會先將它們合并起來，再將合并后的序列和其他序列比對。
選項

--output：輸出文件名稱
--collinear：假設(shè)輸入文件是共線性，即沒有重排

重要參數(shù)的默認(rèn)設(shè)置

--island-gap-size=<number> Alignment gaps above this size in nucleotides are considered to be islands [20]

似乎不允許換行？因為程序太長，我換行了。然后就報錯了···弄成一行后，又正常了

2.2.2 輸出

比對完成后，Mauve會產(chǎn)生多個比對相關(guān)文件。下面介紹一下每個文件包含的信息及對應(yīng)的格式。

2.2.2.1 .alignment文件和XMFA格式(.xmfa)

.alignment文件以 eXtended Multi-FastA格式存儲Mauve產(chǎn)生的比對數(shù)據(jù)。

XMFA格式中存儲了多個共線性塊的比對情況，不同共線性塊以=分割。

每個共線性塊中，一個基因組有一條對應(yīng)的fasta格式的序列。其中，定義行給出這條序列位于基因組的位置和方向（正負(fù)鏈）。

這些共線性塊共同組成了基因組比對結(jié)果。

>seq_num:start1-end1 ± comments (sequence name, etc.)
AC-TG-NAC--TG
AC-TG-NACTGTG
...

> seq_num:startN-endN ± comments (sequence name, etc.)
AC-TG-NAC--TG
AC-TG-NACTGTG
...
= comments, and optional field-value pairs, i.e. score=12345

> seq_num:start1-end1 ± comments (sequence name, etc.)
AC-TG-NAC--TG
AC-TG-NACTGTG
...

> seq_num:startN-endN ± comments (sequence name, etc.)
AC-TG-NAC--TG
AC-TG-NACTGTG
...
= comments, and optional field-value pairs, i.e. score=12345

Mauve略微改造了下xmfa格式，要求基因組序列中的核苷酸能且僅能出現(xiàn)一次。

  • 對于未對齊且在LCB之外的基因組片段，以XMFA中無gap的單條目形式出現(xiàn)。
  • XMFA中的第一個LCB條目會列出所有輸入的基因組序列。如果某些基因組未參與該lCB的比對，其對應(yīng)的坐標(biāo)范圍標(biāo)記為0-0。

2.2.2.2 Islands文件(.bbcols)

.island文件中存儲了比對中發(fā)現(xiàn)的genomic islands，以tab鍵分割。

Island指的是比對中一部分基因組有，而另一部分基因組沒有的區(qū)域。

一個island由一個序列的基因組坐標(biāo)定義，其中另一個基因組在比對的那部分中包含長度為n或更長的缺口。缺口的長度n可以人為定義。是不是很拗口，看下例子就明白了。

Genome 0: ACACGTTCGCTTCGAAA
Genome 1: ACAC------TTCGAA-
Genome 2: ATACGATCGCTTCGTAA

設(shè)定n=5，則稱基因組0和2在5-10位置上有個島
對應(yīng)的.islands文件中，一行記錄一個島。格式：基因組A??島最左側(cè)在A中的位置??島最右側(cè)在A中的位置??基因組B??島最左側(cè)在B中的位置??島最右側(cè)在B中的位置。

0 4 11 1 4 5
1 4 5 2 4 11

第一行記錄在基因組0中，核苷酸4至11與基因組1中的核苷酸4至5對齊。
島的長度=｜(rightA - leftA) - (rightB - leftB)｜。


https://www.mail-archive.com/search?l=mauve-users@lists.sourceforge.net&q=subject:%22%5C%5BMauve%5C-users%5C%5D+islands+file%22&o=newest&f=1

2.2.2.4 backbone文件(.backbone)

什么是骨架，支撐性的東西。

2.2.2.4.1 The original Mauve backbone file

.backbone文件存儲了在所有基因組中都是保守的比對區(qū)域。

保守區(qū)域定義為長度大于等于x，期間gap不超過y個核苷酸的alignment片段。

22256 22371 20147 20299 22255 22370

來自第一個基因組的片段[22,256-22,371]分別在第二個[20,147-20,299]和第三個基因組[22,255-22,370]中保守。

2.2.2.4.2 The Progressive Mauve backbone file

跟original mauve不同的是，progressive mauve的backbone文件不再要求在所有基因組中都是保守的，align regions conserved among two or more genomes。backbone文件按照seq_0_leftend列排序。

The Progressive Mauve backbone file

可以從backbone文件中推出island位置

2.2.2.5 相似度矩陣文件

行和列都是基因組文件，按照輸入順序排序。

0代表沒有一個相同的核苷酸，1代表每個位置的核苷酸都是相同的。

2.2.2.6 SNP文件

記錄了SNP模式（按照輸入順序排序）以及SNP在每個基因組中的位置。

SNP pattern     sequence_1      sequence_2      sequence_3
AAT     5276590 5246627 394
TTC     5276784 5246821 588
AAC     5277418 5247455 1222
MAA     5278225 5248262 2030
AAC     5282804 5252841 6609

2.2.2.7 orthologs文件

0:Z03:2818-3750  1:c04:3512-4444  2::2801-3733    3:ECSE_03:2800-3732
0:Z04:3751-5037  1:c05:4445-5731  2::3734-5020    3:ECSE_04:3733-5019
0:Z05:5251-5547  1:c07:5945-6241  1:c08:6021-6269  2::5234-5530  3:ECSE_05:5233-5529
0:Z06:5700-6476  1:c10:6301-7077  3:ECSE_06:5682-6458

第一行列出了4個同源基因，每個基因分別來自一個基因組。0:Z03:2818-3750中：

  • 0代表基因組index，index跟輸入順序一致。
  • Z03代表注釋到的基因的locus_tag。
  • 2818-3750代表在基因組上的位置范圍。

第2行中，2::2801-3733代表該區(qū)域沒有注釋到的基因，因此沒有l(wèi)ocus_tag，但是是同源的。

第3行中，基因組1上注釋到了兩個基因。

第4行中，基因組2上沒有注釋到基因。

2.2.3 例子

# 比對1, 2 , 3這3個基因組，將輸出保存到threeway.xmfa文件中
progressiveMauve --output=threeway.xmfa genome_1.gbk genome_2.gbk genome_3.gbk

相同的測試文件，Mauve用了大約45分鐘。這明明比Mugsy快好吧···


3. 先pairwise alignment，后合并

有LASTZ+ Multiz和LAST+ Multiz兩種方案。鑒于UCSC用lastz，我就先試試第一種方案。


  1. 受網(wǎng)友Mr_我愛讀文獻啟發(fā)，我才知道還有這一種辦法，感謝感謝
  2. 與LASTZ相比，LAST的輸入基因組不用提前屏蔽掉重復(fù)序列。
    What distinguishes LAST from BLAST and similar tools (e.g. BLAT, LASTZ, YASS)? The main difference is that it copes more efficiently with repeat-rich sequences (e.g. genomes). For example: it can align reads to genomes without repeat-masking, without becoming overwhelmed by repetitive hits.
  3. Multiz是內(nèi)嵌在TBA中的一個reference-biased 多基因組比對工具。UCSC用它來生成多基因組比對。通常，先讓其他序列和reference genome做兩兩比對，再以兩兩比對的結(jié)果為輸入，供Multiz處理。
  4. 參考鏈接：
    很大一部分參考了ucsc genomewiki，寫得很清楚。
    https://darencard.net/blog/2019-11-01-whole-genome-alignment-tutorial/
    http://genomewiki.ucsc.edu/index.php/Whole_genome_alignment_howto
    https://www.biostars.org/p/43225/
    http://seqanswers.com/forums/showthread.php?t=29383
    http://www.bx.psu.edu/miller_lab/dist/tba_howto.pdf

https://github.com/mcfrith/last-genome-alignments/issues/7

3.1 安裝

因為用的是UCSC的自動化流程-doBlastzChainNet.pl，所以需要配置一些東西

  • lastz
  • parasol
  • kent command line binaries
  • 其他不在庫里的一些小程序
    除了lastz，UCSC上有后三種的詳細安裝教程。安裝的時候很容易搜到。

3.2 pipeline

3.2.1 調(diào)整lastz的score文件

lastz_D target[multiple] query --inferonly=lastzd.control --infscores=q 


--inferonly: 表明該程序只用于推斷參數(shù)文件
--infscores: 存儲生成的score文件

control文件選的默認(rèn)的文件：見https://lastz.github.io/lastz/
cat lastzd.control
# base the inference on alignments in the middle half by identity
min_identity       = 25.0%    # 25th percentile
max_identity       = 75.0%    # 75th percentile

# scale scores so max substitution score will be 100, and only use
# alignments scoring at least as well as 20 ideal matches
inference_scale    = 100      # max substitution score
hsp_threshold      = 20*inference_scale
gapped_threshold   = hsp_threshold

# allow substitution score inference to iterate at most 20 times;
# don't perform gap penalty inference -- instead hardwire gap penalties
# relative to max substitution
max_sub_iterations = 20
max_gap_iterations = 0
gap_open_penalty   = 4*inference_scale
gap_extend_penalty = 0.3*inference_scale

# use all seedword positions (don't sample)
step               = 1

# adjust for entropy when qualifying HSPs
entropy            = on

3.2.2 UCSC-doBlastzChainNet.pl

ucsc上的doBlastzChainNet.pl可以自動完成pairwise alignement的所有步驟，包括partition, blastz, cat, chainRun, chainMerge, net, load, download, cleanup, syntenicNet。最后直接得到可用于Multiz的結(jié)果文件。

3.2.2.1數(shù)據(jù)準(zhǔn)備

基因組文件需要先壓縮成gz格式，然后再用faToTwoBit將基因組文件轉(zhuǎn)為.2bit格式，用twoBitInfo獲取各染色體的信息。

# example
gzip ps128.fna
faToTwoBit ps128.fna.gz ps128.2bit
twoBitInfo ps128.2bit stdout | sort -k2,2nr > ps128.chrom.sizes

3.2.2.2 參數(shù)文件準(zhǔn)備

創(chuàng)建一個名為DEF的文件，里面為lastz所用參數(shù)。

TMPDIR目錄要提前創(chuàng)建

參數(shù)意義詳見：http://genomewiki.ucsc.edu/index.php/DoBlastzChainNet.pl

# ps128 vs 299v
PATH=/data/bin:/data/scripts # kent command line binaries存儲位置
BLASTZ=/data/bin/lastz-1.04.00 # lastz
BLASTZ_H=2000
BLASTZ_Y=3400
BLASTZ_L=4000
BLASTZ_K=2200
BLASTZ_Q=/data/lastz/HoxD55.q # score文件，這里要換成我們調(diào)整好的score文件

# TARGET: ps128
SEQ1_DIR=path/ps128/ps128.2bit # target genome
SEQ1_LEN=path/ps128/ps128.chrom.sizes # target chromosome size
SEQ1_CHUNK=32100000 #將target genome切成多大一份
SEQ1_LAP=10000 # target chunk之間的重疊程度
SEQ1_LIMIT=18

# QUERY: 299v
SEQ2_DIR=path/ps128/trackData/299v/299v.2bit
SEQ2_LEN=path/ps128/trackData/299v/299v.chrom.sizes
SEQ2_CHUNK=1000000
SEQ2_LIMIT=2000
SEQ2_LAP=0

BASE=/path/ps128/trackData/299v
TMPDIR=path/dev/shm

3.2.2.3 比對

為了縮短比對時間，將genome切成了chunk。為了充分利用多核系統(tǒng)，UCSC開發(fā)了parasol系統(tǒng)，用于管理任務(wù)。
a. 比對前，要先開啟parasol

initParasol start

b. 比對中

nohup bash -c 'time (doBlastzChainNet.pl `pwd`/DEF -verbose=2 -noDbNameCheck \
 -workhorse=localhost -bigClusterHub=localhost -skipDownload \
   -dbHost=localhost -smallClusterHub=localhost \
     -trackHub -fileServer=localhost -syntenicNet)' > do.log 2>&1 &

c. 比對后

initParasol stop

最后 BASE/axtChain/target.query.synNet.maf為結(jié)果文件，作為Multiz的輸入文件。

3.3.3 Multiz

合并pairwise的結(jié)果：

multiz ps128_299v.maf.gz ps128_ATG-K2.maf.gz 1 out1 out2 > multipleAlign.maf

寫在最后：

  • 第三部分安裝寫的很簡略，實際步驟還是挺多的。
  • 第三部分的parasol只能配置在組里的服務(wù)器上，ParasolLSF可以將它配置在slurm上，但我還沒學(xué)會···

  • 除了上面兩個缺點，doBlastzChainNet.pl還是很好用噠。


最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。
禁止轉(zhuǎn)載，如需轉(zhuǎn)載請通過簡信或評論聯(lián)系作者。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容