前面我们介绍过，病毒的基因组虽然很小，但确是最难拼接的，这看起来似乎不合逻辑，按理说只有几十K应该很容易拼接才是，但如果你真正尝试去拼接就知道了。因为病毒基因组突变率高，并且测序覆盖度高，一般都1000X以上，在加上测序错误，这几个原因综合起来，就让两个片段连接出现太多可能性，导致无法连接，也就连接不起来了。举个简单的例子，我不怕妖魔鬼怪，但是却被一直小小的蚊子搞得一夜难眠。
病毒基因组拼接方法
由于不能直接使用原始测序数据进行拼接，目前病毒基因组主要采用基于参考序列指导（reference guide）的方法，所谓生成一致性序列（consuses）的方法。首先将测序数据与最近源的参考序列进行比对，得到每个位点的覆盖情况，然后利用软件，挑选每个位点最佳的碱基。比如参考序列为A，比对到这个位点上有1000个碱基，其中900个为T，100个为A，那么我们就选择T作为待拼接基因组的碱基类型。最终我们会得到一条与参考基因组长度一致，只是部分位点有所差异的基因组，来作为最终得到的病毒基因组序列。例如下图中，参考序列的基因组为GTCTG，一致性序列的结果为GACTC。
这种生成一致性序列的方法最开始是用于变异检测分析中的方法，后来用于病毒基因组拼接。这种方法显然有很大的问题。这种方法有个前提，就是参考序列一定要非常近源，比如参考序列是29903个碱基，新测序的样品是否有可能是30000个碱基呢，但是利用这种方法得到的基因组序列必然还是29903个碱基，剩下的97个碱基信息就没有了。所以你看到目前拼接出来的基因组都是29903，虽然后面通过碱基Polish可以插入或者删除部分序列，但整体长度变化不大。如果第一个参考序列拼接质量有问题，后面基于这个参考序列得到的新数据，都存在问题。
但是，但是，因为目前只能使用这种方法，没有更好的方法，只能这样。
你也可以拼接基因组
https://github.com/CDCgov/SARS-CoV-2_Sequencing/tree/master/protocols/BFX-UT_ARTIC_Illumina
安装软件
conda install -c bioconda entrez-direct sra-tools bwa samtools ivar quast
1、创建工作目录
mkdir ncov
cd ncov
2、下载数据
#非常关心美国情况，这里下载一个美国的测序数据efetch -db=nuccore -format=fasta -id=MN908947.3 >MN908947.3.fa prefetch SRR11247077 -O ./ fasterq-dump SRR11247077.sra
3、bwa比对
bwa index MN908947.3.fa
bwa mem -t 12 MN908947.3.fa SRR11247077.sra.fastq | samtools sort | samtools view -F 4 -o usa.sorted.bam
4、切除引物
#下载引物列表
wget https://raw.githubusercontent.com/artic-network/artic-ncov2019/master/primer_schemes/nCoV-2019/V1/nCoV-2019.bed
#转换适合ivar格式文件
perl -ne 'my @x=split m/t/; print join("t",@x[0..3], 60, $x[3]=~m/LEFT/?"+":"-"),"n";' nCoV-2019.bed >ARTIC-V1.bed
ivar trim -i usa.sorted.bam -b ARTIC-V1.bed -p usa.primertrim
samtools sort usa.primertrim.bam -o usa.primertrim.sorted.bam
5、生成一致性序列
samtools mpileup -A -d 6000000 -B -Q 0 --reference MN908947.3.fa usa.primertrim.sorted.bam | ivar consensus -p usasus -n N
6、quast评估
#需要下载MN908947.3的基因列表文件gff3格式
quast.py usa.consensus.fa -r MN908947.3.fa --featuresMN908947.3.gff3 --ref-bam usa.sorted.bam --output-dir quast
7、计算覆盖深度
samtools bedcov nCoV-2019.bed.txt usa.sorted.bam