0基础学习基因组三代测序组装（3）——全基因组Survey

之前说到如何对三代测序数据做污染评估，取随机序列做blastn比对nt库，确定物种分布情况。实际blast比对还要考虑比对的序列长度和ONT本身数据错误率，以及结合GC-depth确定是否有污染。基因组三代测序数据组装之前，我们还要做一个全基因组survey。主要是为了减少盲目性，先做低深度的基因组分析，也是初步了解物种基因组特征的有效方法，比如评估基因组大小和杂合情况，为后续全基因组de novo组装策略指定提供指导。

基因组复杂程度的经验性标准：

简单基因组: 单倍体；或纯合二倍体；或杂合度低于0.5%, 且重复序列低于50%, 且GC含量在35%-65%的二倍体。
复杂基因组: 杂合度在0.5%~1.2%之间，或重复序列高于50%，或GC含量异常(<35%或>65%)的二倍体，或者多倍体。复杂基因组可以采用“2+3”即二代和三代测序技术相结合，加之Hi-C辅助组装的组装策略。
高复杂基因组: 杂合度>1.2%；或重复序列占比大于65%。

有条件的话，也可以用流式细胞仪对基因组大小做个预估。我这里只有二代基因组测序数据，因此用基因组二代测序数据做全基因组survey。当然，这里要注意一点，做全基因组survey的样本和后续de novo组装的样本要来自同一个个体，避免个体间基因组特征的差异。

1 原始数据质控

因为是对二代测序数据进行分析，质控的过程本质上和之前处理转录组二代数据一样，这里只提下过程和结果。

1.1 fastqc生成质控报告

1	fastqc *.fq.gz -o ./

二代测序是双端测序结果，我这里只截图了部分qc报告，可以看出GC含量比较稳定，测序质量也比较高。

1.2 trim-galore数据过滤

报告中的结果虽然好，但是还是需要过滤一遍，把末端接头adapter序列过滤掉。

1	trim_galore -q 25 -phred33 -length 100 -stringency 1 -paired -o clean_data 1_raw_1.fq.gz 1_raw_2.fq.gz

参数在这篇博客转录组数据分析笔记（1）——如何用fastqc和trim-galore做测序数据质控有提到，这里不再赘述。

看下report文件，过滤了Q值25以下的reads和adapter序列

2 k-mer分析

先说一下k-mer的概念：k-mer在这里指将reads迭代拆分成包含k个碱基的序列（类似blast中的word length，蛋白质是3，核酸是11），我们后面要分析的基因组特征都是基于k-mer分布基础上进行的。

基因组大小可以通过总 （K-mer 数量）/（K-mer 期望测序深度）来估计
k-mer分布曲线的主峰所在横坐标可以作为期望的测序深度
测序覆盖均匀、不存在测序误差和基因组重复序列的理论条件下，K-mer分布曲线会符合泊松分布
单倍体或纯合基因组的 K-mer 分布曲线只有一个主峰
杂合二倍体基因组的 K-mer 分布曲线有两个峰，分别为杂合峰（主峰1/2处）和纯合峰（主峰），前者深度只有后者的一半
重复序列含量较高时会在主峰后面形成一个重复峰（主峰的2倍处）或者形成拖尾
一般选择17-mer评估基因组大小，因为ATCG组成长度为17的核酸序列，理论上有4的17次方种可能，足以覆盖一般的正常基因组。为了避免回文序列，K-mer分析选择K长度均为奇数。

2.1 安装jellyfish

根据上面说的k-mer概念，可以理解k-mer分析是非常耗计算资源的。我们要自己用脚本实现的话，需要将十几个G的reads分割成不同长度片段，再统计出现的次数，耗时而且麻烦。jellyfish是一款统计DNA序列中Kmer的分布的软件，它运行速度快，内存消耗低，支持并行，也是用的最多的统计k-mer的软件。

~~重点是可以通过conda直接安装……~~最好不要用conda安装，我之前运行了1天没出结果也没报错（一度怀疑我的参数设置是不是有问题），百思不得其解。后来从github上重新下载，编译和安装之后，不到10分钟就跑出结果了…我不知道两种安装方式有什么区别，这里就记录下自己踩的坑。

因为jellyfish不支持.gz的压缩文件，所以之前用tram galore过滤后得到的clean reads需要用gunzip命令解压。

1	conda install -c bioconda jellyfish # 可以用conda安装，我运行的时候出了问题，暂未解决，不推荐

点击这里进入jellyfish的github下载地址

我们用本地安装的方式，先下载tar.gz的源码包，tar -zxvf解压后进入jellyfish-2.3.0文件夹。

我是集群登录的，下面讲的步骤都是在集群上操作（非root账户）

# 第一步检测。本质上是一个shell脚本，根据系统环境产生合适的makefile文件或者C的头文件（.h结尾的文件），非root账户下--prefix后面接上自己账户的绝对路径。
./configure --prefix=/public/home/wlxie
# 第二步编译。对源代码包进行编译，如果有错误自己看是否有依赖库的缺失，主要是这个问题。
make
# 第三步安装。如果前面没有指定自己账户的路径，这一步是会报错没有权限的（用户不能向系统目录写入文件）。
make install
# 第四步自检。
make check

make和make check这两步都会因为动态链接库命名不同，导致报错无法找到动态库；以及我在检测通过之后，用集群运行程序仍然出现了动态库的某个模块无法调用的情况。这里统一说下解决方法。

前面configure会在我们的家目录下生成bin、lib和share目录，这里比较重要的是bin和lib目录。我们运行的命令在bin目录里，对应要改环境变量PATH；而需要调用的动态库是在lib目录下，对应要改环境变量LD_LIBRARY_PATH。家目录下的.bashrc文件加入以下内容

1 2	export PATH="/public/home/wlxie/bin:$PATH" export LD_LIBRARY_PATH="/public/home/wlxie/lib:$LD_LIBRARY_PATH"

添加之后保存退出，并且source ~/.bashrc刷新一下系统环境变量。

我碰到的报错是libcrypto.so.1.0.0和libstdc++.so.6这两个动态库找不到，但是locate命令查看这两个动态库，在系统目录/lib64/下都能找到文件，因此将这两个动态库文件直接复制到家目录的lib文件夹，问题就全部解决了。

如果libstdc++.so.6报错某版本的文件不存在，可以先到动态库目录下，运行strings命令查看动态库中是否有对应的文件：

strings /lib64/libstdc++.so.6 | grep CXXABI   # 比如找不到GLIBCXX_3.4.26，看看动态库中是否存在这个版本的文件，如果不存在，更新动态库；如果存在但是找不到，建议直接拷贝到自己的lib目录下

make check之后会生成一个日志文件test-suite.log，没有fail的项目说明软件安装成功，没有问题。

2.2 k-mer频数分布

# k-mer计数
jellyfish count -m 17 -s 300M -t 50 -C -o 17-mer.jf ./1_raw_1_val_1.fq ./1_raw_2_val_2.fq
# k-mer频数统计
jellyfish histo -t 4 17-mer.jf > 17-mer.histo
# 统计总k-mer数和特征k-mer数等
jellyfish stats 17-mer.jf -o counts_stats.txt

记录一下各个参数的意义：

-m # k-mer长度设置为17bp，进行计数

-s # 存储用的hash表大小，说实话我没看懂什么意思，基因组估计有多大就用多大就是了，单位是M或者G

-t # 使用的线程数，也就是cpu核数

-C # 大写的C，对正负链reads都进行统计，双端测序一定要加这个参数

-o # 结果文件的前缀名，结果文件是一个二进制文件

正常来说，10分钟就能跑完程序并给出k-mer计数结果文件。我用conda安装的jellyfish同样条件运行了20个小时没有结束……而且还不报错！第一次运行这个软件，没有人参考和交流，百度到的教程都是抄来抄去的也没有人说明时间的问题……以后还是去官网安装生信软件了，虽然麻烦一点但是靠谱……

这个软件的帮助文档在/jellyfish-2.3.0/doc目录下，所有功能和参数都有英文详解。

k-mer频数统计是在计数结果文件上进一步统计各个k-mer出现的次数，频数统计结果文件17-mer.histo将k-mer从1统计到10000，最后一行是10001以后对应的总频次。counts_stats.txt是总的统计结果，包括k-mer总数（Total），特异的k-mer数目（Distinct）只出现过一次的k-mer数量（Unique），频数最高的k-mer数量（Max_count）四项。

有了频数统计结果文件17-mer.histo就可以用R作图了，以下R作图代码来自于CSDN博主生信技工

kmer <- read.table('17-mer.histo')
kmer <- subset(kmer, V1 >=5 & V1 <=500)     # 只取5-500bp长度的k-mer统计频次
Frequency <- kmer$V1
Number <- kmer$V2
png('kmer_plot.png')
plot(Frequency, Number, type = 'l', col = 'blue')
dev.off()       # 保存png

k-mer分布图如下，当然这只是一个简略图，上面R作图代码还有很多细节可以补充

2.3 基因组大小、重复率、杂合率估算

横坐标表示k-mer深度，纵坐标为k-mer数量，可以看得出来测序的样本是个杂合二倍体。主峰坐标（116，2584902），杂合峰坐标（57，1188461），也就是说k-mer期望深度为116；k-mer总数为34655456060；主峰2倍深度也就是232之后的k-mer为重复序列k-mer，总数可以通过导出17-mer.histo文件进行统计（改成csv格式直接两步算出），共16361378388

k-mer分布曲线中无异常峰，说明二代测序提取的DNA纯度较高，没有被污染
根据（K-mer 数量）/（K-mer 期望测序深度）估算基因组大小为298M。去除深度小于5的错误k-mer，估算基因组大小为292M.
根据（重复序列的k-mer总数）/（K-mer 期望测序深度）估计重复序列大小为141M，即重复率48.29%
单拷贝序列大小U=292-141=151M，要计算杂合率，需要统计非重复k-mer的总数，也就是计算杂合峰面积，建议还是用软件或者在线工具比如genomescope2.0

jellyfish + GenomeScope是一套应用非常广泛的基因组survey方法，GenomeScope2.0适合用于分析二倍体生物。