理论部分
该部分参考及引用陈连福的生信讲义，zhaofei的https://vip.biotrainee.com/d/103--。
由于仪器测序读长的限制等，在建库时会将DNA随机打断为小片段的序列，因此，基因组组装就是将小片段的序列连接起来，但是序列之间的联系十分复杂，常用构建Graph来进行表示，然后在对Graph进行简化、拼接，即reads→contigs→scaffolds。
下图即为简单的示意图，其中顶点A,B,C,D,E,F称为nodes，是6条小片段序列；连接2个nodes的有方向的箭头称为edges，这些边表示reads间的overlap;而数字代表着权重。通过对Graph进行分析，选择权重大的来简化Graph，得到ABCDEF这个序列，依此类推，将基因组产生ABCDEFGH...等若干的reads，再根据各reads间的重复区域，选出最优路径，形成contigs，再组成scaffolds，最终得到基因组序列。
组装示意图.png
寻找最优路径的常用算法：
早期使用的是贪婪法 ：先选定初始read，找到和其重叠区域最高的read进行延伸，直至拼接后的read两端不能再进行延伸为止。每次延伸都是从最优匹配开始的，贪婪法得到的是局部最优解，而不是全局最优解，因此，在遇到重复序列时会出现比较大的问题。
Overlap-Layout-Consensus(OCL) 常用语处理reads读长较大的测序数据，比如PacBio数据的组装。OLC算法分为三步：1）对所有的reads进行两两比对，得到overlap。2）根据overlap简化graph，建立overlap图，将reads组合成contig。3）得到consensus：将所有read序列排列起来，找到一条从起始节点到终止节点的最佳近似路径使得最终路径将会遍历一次重叠区域的每个节点，从而得到目标基因序列。
de Bruijin Graph（DBG）是目前常用的二代测序拼接算法。相关软件：ALLPATHS-LG、SOAPdenovo等。该算法和OLC类似，不同之处在于：该算法中的nodes是kmer序列，kmer和kmer必须仅有一个碱基差异才能相连，即相邻的kmer序列是通过k-1个碱基连接到一起的（每次只移动一个位置）。这种算法降低了重复区域的复杂度，降低了内存消耗。其步骤：1）构建DBG图，将reads分割为一系列连续的kmer；2）合并DBG图；3）构建contig：寻找最优路径（经过每一个节点且仅经过一次），最优路径对应的碱基序列构成一个contig；4)构建scaffold：通过PE reads位置确定contig之间的相对位置和方向，组装contig，填充contig之间的gap，得到scaffold序列。
由于DBG构建不需要reads具有endanger的长度，因此只适用于reads长度较短的Illumina测序数据。也因此DBG对于重复区域比较狠分析，进行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的区别：管道是实现“将前面的标准输出作为后面的标准输入”，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文件中第九列，统计每个插入片段长度的个数
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之间没有直接的联系，其会在运行过程中运用一系列的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软件中默认是将大片段文库的方向改变为inward。
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文件中几个目录的路径的设置。
接着就是漫长的结果等待啦~~~
由于组装需要的内存很大，一定要保证内存，否则会运行到一半会被删除。
网友评论