如何检查叶绿体基因组组装结果的准确性v1

  有小伙伴之前问过我，怎么判断叶绿体基因组组装结果是否正确。其实判断方法主要还是依靠测序数据，基因注释只能作为辅助判断，无法起决定作用。将测序数据比对到组装的序列上，看看有没有没有覆盖到的情况，但是即使全部覆盖到了，也可能有问题，比如从一个完整的组装结果上任意截取两个不连续的子序列，然后拼接到一起，你会发现新的序列中每个碱基都会有覆盖，但是实际上连接的地方是有问题的，没有reads能够将两者连接到一起。所以直接看覆盖深度的图是无法判断每个位置是否正确，这种可以放到IGV等可视化软件中放大些去看，虽然基因组不大（150k左右）但是放大看还是比较辛苦的。为此，写了一个简单的脚本来帮助判断覆盖情况和连续性。思路如下：

1. 以步长为1，将组装的序列打断成长度为k的子序列（kmer）
2. 搜索测序数据中包含上述子序列的reads，统计该kmer出现的频数
3. 在原始序列上展示每个位点的覆盖情况

  如果有懂组装原理的小伙伴就可以发现其实这是个逆向组装的过程，通过k-1长度kmer的overlap把序列连接连接到一起，如果连接不起来，说明序列拼接有可能是错误的。当然，如果连接起来了，也不能说明你的序列一定对（比如可以发生重组，或者较长重复片段的地方，或者低覆盖区域的地方，但是对于大部分简单的序列如果能连接起来，说明序列是没问题的）。为了方便使用，脚本采用纯perl语言编写，不依赖其他软件。有条件的小伙伴可以先把序列mapping下组装结果，用mapping上的序列做检查速度会更快一些。
  用法：

perl check_ass_v1.pl   -h

Usage:
                perl check_ass_v1.pl -g your_genome.fasta -d NGS.1.fq NGS.2.fq -o check -k 75
  Options:
        -g [str]        genome fasta file
        -d [str]        NGS data fastq format(.fq/.gz)

        --------
        -k [int]        kmer length, default: 125
        -c!     genome sequence is not circler
        -f!     fast mode
                        The speed is fast, but the number of kmer obtained will be less
        -o [str]        outdir, default : ./check
        
        -h      Help
---------------------------------
必选参数：
	-g : 输入的组装结果，fasta格式
	-d : 输入的测序数据，fastq格式，支持gz格式，采用PerlIO::gzip模块读取
可选参数：
	-k : kmer大小，默认125，一般默认即可
	-c : 如果加了 -c 选项（不接具体参数），则把输入的序列当成线性序列处理，否则当成环状
	-f : -f 选项（不接具体参数），较快的模式，如果读入的一条测序reads前三个kmer都不在组装的序列中，则跳过该reads。建议数据量较大的时候使用（双端大于4G）。
	-o : 输出的文件夹，默认为：check。里面包含了覆盖情况，以及用到的reads
	

  结果解读：
  

K125_check.xls:
	kmer length:125
	#NLT_S30_L003_R1_001.fastq.gz:
	#seq_id base    base_in_data    kmer_nu
	Tn      0       C       C       15
	Tn      1       T       T       18
	Tn      2       T       T       18
	Tn      3       A       A       20
	Tn      4       A       A       21
	Tn      5       T       T       24
	Tn      6       G       G       25
	Tn      7       T       T       34
	Tn      8       T       T       39
	Tn      9       A       A       43
	Tn      10      A       A       49
	Tn      11      T       T       54
	Tn      12      G       G       56
	Tn      13      T       T       56
	Tn      14      G       G       58
	Tn      15      T       T       58
	Tn      16      T       T       62
#共五列，第一列是序列id；第二列是碱基位置，从0开始；第三列是序列中的碱基；
#第四列是测序reads中对应的碱基（对应kmer的第一个碱基），如果不存在，则用N；
#第五列是该kmer出现的频率。
#所以最终搜索第四列有没有N即可。如果中间有N（kmer_nu=0），并且N上下非N的碱基对应的kmer_nu很高，由一个比较高的值突然变成0，那么该处的连接是有问题的。如果kmer_nu是缓慢变成0，然后缓慢递增，说明该处可能只是覆盖深度比较低。如果是末尾有N，说明序列不成环（不加-c选项）
K125_mapped.fa：
	#比对上的reads，暂时用处不大，后面打算写个v2版本，把每个kmer对应的四种碱基频率都列出来会用到这个~

  因为测序数据较大，所以总体运行时间可能较长，3m~1h都有可能，快速模式可以跳过很多无用reads，节省了不少时间，但是也可能会漏掉一些有用的reads（问题不大），建议使用-f参数，提高速度。此外，kmer较大，循环的次数会减少，也会提高一些效率，用125就行。
  该脚本适用于初学者，为了便于交流，加Q群：936427018，免费获取脚本和帮助。