理論部分

該部分參考及引用陳連福的生信講義沥曹，zhaofei的https://vip.biotrainee.com/d/103--。

由于儀器測序讀長的限制等窄锅，在建庫時會將DNA隨機打斷為小片段的序列藻肄，因此，基因組組裝就是將小片段的序列連接起來间学，但是序列之間的聯(lián)系十分復雜殷费，常用構建Graph來進行表示，然后在對Graph進行簡化低葫、拼接详羡，即reads→contigs→scaffolds。

下圖即為簡單的示意圖氮采，其中頂點A,B,C,D,E,F稱為nodes殷绍，是6條小片段序列；連接2個nodes的有方向的箭頭稱為edges鹊漠，這些邊表示reads間的overlap;而數字代表著權重主到。通過對Graph進行分析，選擇權重大的來簡化Graph躯概，得到ABCDEF這個序列登钥，依此類推，將基因組產生ABCDEFGH...等若干的reads娶靡，再根據各reads間的重復區(qū)域牧牢，選出最優(yōu)路徑，形成contigs姿锭，再組成scaffolds塔鳍，最終得到基因組序列。

組裝示意圖.png

尋找最優(yōu)路徑的常用算法：

早期使用的是貪婪法 ：先選定初始read呻此，找到和其重疊區(qū)域最高的read進行延伸轮纫，直至拼接后的read兩端不能再進行延伸為止。每次延伸都是從最優(yōu)匹配開始的焚鲜，貪婪法得到的是局部最優(yōu)解掌唾，而不是全局最優(yōu)解放前，因此，在遇到重復序列時會出現比較大的問題糯彬。

Overlap-Layout-Consensus(OCL) 常用語處理reads讀長較大的測序數據凭语，比如PacBio數據的組裝。OLC算法分為三步：1）對所有的reads進行兩兩比對撩扒，得到overlap似扔。2）根據overlap簡化graph，建立overlap圖却舀，將reads組合成contig虫几。3）得到consensus：將所有read序列排列起來，找到一條從起始節(jié)點到終止節(jié)點的最佳近似路徑使得最終路徑將會遍歷一次重疊區(qū)域的每個節(jié)點挽拔，從而得到目標基因序列辆脸。

de Bruijin Graph（DBG）是目前常用的二代測序拼接算法。相關軟件：ALLPATHS-LG螃诅、SOAPdenovo等啡氢。該算法和OLC類似，不同之處在于：該算法中的nodes是kmer序列术裸，kmer和kmer必須僅有一個堿基差異才能相連倘是，即相鄰的kmer序列是通過k-1個堿基連接到一起的（每次只移動一個位置）。這種算法降低了重復區(qū)域的復雜度袭艺，降低了內存消耗搀崭。其步驟：1）構建DBG圖，將reads分割為一系列連續(xù)的kmer猾编；2）合并DBG圖瘤睹；3）構建contig：尋找最優(yōu)路徑（經過每一個節(jié)點且僅經過一次），最優(yōu)路徑對應的堿基序列構成一個contig答倡；4)構建scaffold：通過PE reads位置確定contig之間的相對位置和方向轰传，組裝contig，填充contig之間的gap瘪撇，得到scaffold序列获茬。

由于DBG構建不需要reads具有endanger的長度，因此只適用于reads長度較短的Illumina測序數據倔既。也因此DBG對于重復區(qū)域比較狠分析恕曲，進行de novo組裝時需要構建大片段文庫，測MP數據渤涌，只要MP文庫長度大于重復序列長度码俩，則有利于重復序列的組裝。

kmer和內存

k值越大可辨別更多的小重復序列歼捏，越容易把DBG轉換為唯一的序列稿存，但得到的拼接過程含有更多的gaps；小的k值對應的DBG能夠得到較好的連通性瞳秽，但是算法的復雜度會提高瓣履，repeats序列處理會更復雜，增加了錯拼的可能性练俐。

kmer越大需要的內存就越多袖迎，所以計算機的內存大小也會限制kmer的取值。這里需要說明的是腺晾，輸入數據的多少不會影響memory用量燕锥，但是輸入數據的錯誤越少，占用的內存也就越少悯蝉，假設所有測序數據都沒有任何錯誤归形，那么DBG的大小并不會因為測序深度的增加，因為不需要將因為幾個堿基不一致的kmer存入到DBG中（下一部分會具體提到）鼻由。至于需要多大的RAM則取決于DBG的大小和組裝基因組的大小暇榴。

另外，在拼接的過程中盡量避免使用偶數kmer蕉世，否則容易是kmer產生回文序列蔼紧，特別是在鏈特意性的數據中。

在平時分析中狠轻，一般會設置一個kmer的梯度（21奸例，23，25向楼，27,2931）查吊，來解決DBG算法loss of read coherence的問題。然后從中選擇最好的結果蜜自。另外菩貌，還有一種說法是在進行拼接過程時，kmer應該選擇read長度的1/2到2/3大小重荠，否則可能拼接出過多的Contig箭阶。

String Graph能較好的組裝散在重復序列。

目前構建Graph的方法主要有3種：Overlap-Layout-Consensus(OCL),de Bruijin Graph（DBG）和String Graph戈鲁。

實操部分

使用ALLPaths-LG組裝進行基因組組裝仇参，適合于短reads數據，也是現在公認的進行基因組De novo 組裝效果最好的軟件婆殿。但是該軟件十分消耗內存和計算诈乒。

1. fastqc軟件使用查看數據原始質量

find ./ -name "*.gz" |xargs fastqc -t 3 -o ./result 

xargs命令的用法：

其中管道和xargs的區(qū)別：管道是實現“將前面的標準輸出作為后面的標準輸入”，xargs是實現“將標準輸入作為命令的參數”婆芦∨履ィ可以試試下面兩代碼的結果

echo "--help"|cat #此處結果是顯示 “--help”
echo "--help"|xargs cat #此處結果是顯示cat的幫助說明文檔 相當于 cat --help喂饥，即xargs是將內容作為普通的參數傳遞給程序
cat gzip.sh|xargs -i echo "quick_qsub {}"

• 結果說明

image

一般是查看網頁版結果，結果解釋說明fastqc中文結果說明肠鲫。

fastqc 軟件主要是針對全基因組測序的员帮，并且各建庫方法不同，其評判標準也會有所差距导饲；不能只是一味的尋求全部結果通過捞高。

ls *zip|while read id;do unzip $id;done #批量解壓壓縮結果

從文件夾中批量抓取里面的%GC,Total sequences等信息

Q20：質量值大于20的堿基數目占總共堿基數目的比例。

• 使用multiqc批量查看fastqc結果

multiqc *fastqc.zip --pdf #

image

2. nxtrim文庫分選

nxtrim分開PE和MP文庫渣锦。
該軟件會用于去除Nextera Mate Pair 文庫中接頭并且根據接頭位置的方向分類reads硝岗。
其中的每個read都會搜索接頭信息（CTGTCTCTTATACACATCT）和他的反向序列（AGATGTGTATAAGAGACAG），因此這個完整的連接接頭是CTGTCTCTTATACACATCT+AGATGTGTATAAGAGACAG袋毙。

nxtrim -1 reads1 -2 read2 -O folderfile --joinreads --preserve-mp --separate 1 > read_nxtrim.log 2 >&1

注：該軟件的默認參數--rf 會將MP文庫的reads方向變換回去型檀。

nxtrim.png

其中的MP和PE 的文件即為所得。

3. bwa比對得到插入片段文庫

bwa mem -M ref.fa reads1.pe.fq.gz reads2.pe.fq.gz > read.pe.sam
bwa mem -M ref.fa reads1.mp.fq.gz reads2.mp.fq.gz > read.mp.sam
perl count_num_sam.pl read.pe.sam
perl count_num_sam.pl read.mp.sam 

SAM文件中第9列是建庫時候的打斷的片段長度娄猫，如若使用的是PE150的數據贱除，那么打斷成350bp，則這里的數據應該是350個字符左右媳溺。

#!/usr/bin/perl   #count_num_sam.pl;提取sam文件中第九列月幌，統(tǒng)計每個插入片段長度的個數
open FH,'>',"$ARGV[0].count.out" or die;
my @num;
my $i = 0;
while(<>){
        next if /^\@/;
        if(/(?:.*?\s){8}(-?\d+)\s.*/){
                $num[$1]++ if($1>=0);
        }
}
$num[0] = $num[0]/2;
for($i=0; $i<20001; $i++){
        if($num[$i] == 0){
                print FH "$i\t0\n";
        }else{
                print FH "$i\t$num[$i]\n";      
        }        
}

用excel得到插入片段的圖，平均值計算方法是：(A+B)/2;偏差值=平均值-A悬蔽。

插入片段長度圖.png

根據實驗中選擇文庫的原理（如下圖）扯躺，對于同一個Well的數據，其index相同蝎困，則其MP文庫插入片段分布一致录语，此類數據只需分析一個插入片段的圖即可。對于同一個 Fragment禾乘，其PE文庫插入片段分布一致澎埠，同理此類數據只需分析一個插入片段的圖即可。

文庫構建.jpg

4. ALLPaths-LG組裝（參考陳連福生信講義）

使用ALLPATHS-LG的條件比較嚴格始藕，有以下的注意事項：

1.不能只使用一個library數據進行組裝蒲稳；

2.必須有一個“overlapping”的片段文庫的PE數據；

3.必須有jumping library 數據（也就是Illumina Mate-pair測序數據）

4.基因組組裝需要有100X或者以上基因組覆蓋度的堿基伍派，這里的覆蓋度是指raw reads數據的覆蓋度江耀；

5.可以使用PacBio數據；

6.不能使用454數據和Torrent數據诉植；

7.輸入10G的堿基數據量祥国，大約需要17G內存；

8.對于試探性的參數晾腔，比如K舌稀，原則上可以調整啊犬；但是一般不自行調整，也不推薦扩借。本軟件中Kmer大小的參數K和read之間沒有直接的聯(lián)系椒惨，其會在運行過程中運用一系列的K值。

4.1 準備in_groups.csv和in_libs.csv文件

in_groups.csv用于指出測序數據的存放路徑潮罪，

其中file_name：數據文件所存放位置，文件名可以包含‘*‘和’?’领斥，從而代表paired數據嫉到。支持的文件類型有“.bam,.fasta,.fa,.fastq,.fq.,fastq.gz,.fq.gz”。

file_name, library_name, group_name
HoS150_peQ20-1_R?.cut.fq.trimmed.paired.fastq,     PE1,       PE01
HoS_mp1_R?.fastq, MP6-1,      HoS250_6-1
HoS_mp2_R?.fastq, MP6-2,      HoS250_6-2
HoS_mp3_R?.fastq, MP6-3,      HoS250_6-3

in_libs.csv用于給出文庫的特性月洛。

其中：library_name指文庫的名字何恶，和In_groups.csv相匹配。type文庫類型：fragment→PE測序嚼黔；jumping→MP測序细层；long→Pacbio測序。read_orientation reads的方向唬涧，小片段文庫為inward疫赎，大片段文庫為outward。但是需要注意的是nxtrim軟件中默認是將大片段文庫的方向改變?yōu)閕nward碎节。

library_name, project_name,     organism_name,        type,     paired,   frag_size, frag_stddev, insert_size, insert_stddev,    read_orientation, genomic_start, genomic_end
PE1, HoS, RFgenome, fragment, 1, 287, 50,, , inward,   0, 0
MP6-1,  HoS, RFgenome,  jumping, 1,  , ,  2290,220,inward,  0, 0
MP6-2,  HoS, RFgenome,  jumping, 1,  , ,  2808,  318,  inward,  0, 0
MP6-3,  HoS, RFgenome,  jumping, 1, , , 3954, 750,  inward,  0, 0

4.2 將Fastq轉換為ALLPATH-LG支持的輸入格式

#!/bin/sh
#這個為prepare.sh文件
# ALLPATHS-LG needs 100 MB of stack space.  In 'csh' run 'limit stacksize 100000'.
ulimit -s 100000 #程序中需要較大的堆椗醺悖空間，使用ulimit將計算資源放寬松些狮荔。

mkdir -p /02allpaths/data #用于將轉換后的數據文件放入到此目錄下胎撇。

# NOTE: The option GENOME_SIZE is OPTIONAL. 
#       It is useful when combined with FRAG_COVERAGE and JUMP_COVERAGE 
#       to downsample data sets.
#       By itself it enables the computation of coverage in the data sets 
#       reported in the last table at the end of the preparation step. 

# NOTE: If your data is in BAM format you must specify the path to your 
#       picard tools bin directory with the option: 
#
#       PICARD_TOOLS_DIR=/your/picard/tools/bin

PrepareAllPathsInputs.pl\

 DATA_DIR=$PWD/genome/data \
 PLOIDY=2\ #生成ploidy文件，1表示基因組為單倍體型殖氏；2為雙倍體型晚树。
 IN_GROUPS_CSV=in_groups.csv\
 IN_LIBS_CSV=in_libs.csv\
 OVERWRITE=True\
 | tee prepare.out 

運行以上命令，將fastq文件轉成運行ALLPATH-LG所需要的文件雅采，并存放到/02allpaths/data文件夾下爵憎。

4.3 ALLPATHS-LG主程序的使用

使用RunAllPathsLG這個命令來進行基因組組裝，該命令有很多參數总滩，但是在一般情況下不要隨意使用纲堵，使用默認設置即可。

#!/bin/sh
# assemble.sh
# ALLPATHS-LG needs 100 MB of stack space.  In 'csh' run 'limit stacksize 100000'.
ulimit -s 100000

RunAllPathsLG\
 PRE=$PWD\ #程序運行的根目錄闰渔，所有的其他目錄全在該目錄下
 REFERENCE_NAME=.genome \ #參考基因組目錄名稱席函，位于PRE目錄下，若有參考基因組冈涧，可放于此目錄下茂附。
 DATA_SUBDIR=data\ #DATA子目錄正蛙，位于REFERENCE_NAME目錄下，程序從該目錄下讀取數據营曼。
 RUN=run\  #位于DATA_SUBDIR目錄下乒验，將生成的中間文件和結果文件存儲與該目錄中。
 SUBDIR=test\
 TARGETS=standard\
 OVERWRITE=True\
 | tee -a /02allpaths/assemble.out

注意：assemble.sh文件中幾個目錄的路徑的設置蒂阱。

接著就是漫長的結果等待啦~~~

由于組裝需要的內存很大锻全，一定要保證內存，否則會運行到一半會被刪除录煤。