拼接病毒基因组

生信喵

背景
      得到参考序列之后，由于病毒基因组较小，不容易从头拼接，因此可以采用与参考序列比对生成一致性序列的方法。无论是 PCR 扩增产物还是宏基因组测序都可以使用该方法得到新冠基因组。
1 病毒基因组拼接原理
      由于不能直接使用原始测序数据进行拼接，目前病毒基因组主要采用基于参考序列指导（reference guide）的方法，所谓生成一致性序列（consensus）的方法。
      首先将测序数据与最近源的参考序列进行比对，得到每个位点的覆盖情况，然后利用软件，挑选每个位点最佳的碱基。比如参考序列为 A，比对到这个位点上有 1000 个碱基，其中900 个为 T，100 个为 A，那么我们就选择 T 作为待拼接基因组的碱基类型。最终我们会得到一条与参考基因组长度一致，只是部分位点有所差异的基因组，作为最终得到的病毒基因组序列。例如下图中，参考序列的基因组为 GTCTG，一致性序列的结果为 GACTC。
      
      生成 consensus 序列原理
      这种生成一致性序列的方法最开始是用于变异检测分析中的方法，后来用于病毒基因组拼接。这种方法显然有很大的问题。因为该方法有个前提，就是与参考序列一定要非常近源，比如参考序列是 29903 个碱基，新测序的样品是否有可能是 30000 个碱基呢，但是利用这种方法得到的基因组序列必然还是 29903 个碱基，剩下的 97 个碱基信息就没有了。所以你看到目前拼接出来的基因组都是 29903，虽然后面通过碱基 Polish 可以插入或者删除部分序列，但整体长度变化不大。如果第一个参考序列拼接质量有问题，后面基于这个参考序列得到的新数据，都存在问题。

2 拼接新冠病毒
      接下来我们将采用一致性序列的方法新冠病毒，该方法主要分为三步：
      1、下载近源参考序列和模板文件；
      2、将测序数据与参考序列进行比对；
      3、得到每个位点最大概率位点信息，生成一致性序列。
2.1 下载测序数据
      下载参考序列时需要注意一点，目前新冠病毒通用的参考序列为 NC_045512，这是 refseq数据库的 ID，原始 genbank 库 ID 为 MN908947.3，二者序列是完全一致的。不过由于后面使用的引物文件 bed 文件中 ID 为 MN908947.3，所以，参考序列最好使用 ID 为MN908947.3 的序列。

1. #下载参考序列
   
2. efetch -db nuccore -format fasta -id MN908947.3 > MN908947.3.fa
   
3. #下载illumina测序数据
   
4. wget https://storage.googleapis.com/nih-sequence-read-archive/run/SRR14867660/SRR14867660 -O SRR14867660.sra
   
5. fasterq-dump SRR14867660.sra

复制代码

     下载引物 bed 文件，引物文件可以使用 Artic Network 提供的文件。

1. git clone https://github.com/artic-network/artic-ncov2019.git

复制代码

2.2 比对
      比对就是将测序数据比对到参考序列上，形成一种堆叠效果，可以检测到参考序列每个位点的比对细节，用于下一步生成一致性序列。由于测序技术的差别，比对也有一些差别。对于短读长测序，可以使用 bwa，bowtie2 等比对工具，当然也可以使用 minimap2 的-x sr比对方法。对于 nanopore 以及 pacbio 等长读长测序，推荐使用 minimap2 软件。

1. #illumina 测序数据比对，排序一步完成
   
2. bwa-mem2 index MN908947.3.fa
   
3. bwa-mem2 mem -t 12 MN908947.3.fa SRR14867660_1.fastq
   
4. SRR14867660_2.fastq | samtools sort | samtools view -F 4 -o
   
5. ncov.sorted.bam

复制代码

     比对之后 samtools 的处理方法一致，主要是转换 bam，排序，建立索引，如果使用 PCR扩增，还可以增加一步剪切掉引物，此步骤需要一个 bed 区间文件，给出 PCR 具体位置信息。PCR 引物分为 V1，V2，V3 三个版本，如果是使用 V1 版本的 bed 文件，由于里面分隔符不是制表符，还需要一步转换过程，不过 V2，V3 版本的 bed 文件不需要此步骤。目前的数据都是使用 V3 版引物。对于宏基因测序，无需该步骤，直接生成一致性序列即可。

1. #切出引物部分
   
2. ivar trim -e -i ncov.sorted.bam -b /share/home/xiehs/04.ncov/31.assembly/artic-ncov2019/primer_schemes/nCoV-2019/V3/nCoV-2019.bed -p ncov.primertrim
   
3. #排序建立索引
   
4. samtools sort -@ 12 -o ncov.primertrim.sorted.bam ncov.primertrim.bam
   
5. samtools index ncov.primertrim.sorted.bam
   
6. #统计覆盖度
   
7. samtools coverage ncov.primertrim.sorted.bam -o ncov.samcov.txt -m

复制代码

2.3 生成一致性序列
      生成一致性序列是病毒基因组拼接中最重要的一个步骤，由于目前还没有同一的标准，使用软件众多，就会造成结果的不一致。尤其是在处理比对过程中的插入缺失问题，不同的软件可能采用不同的策略。另外，有些人采用类似基因组“纠错”的方法生成一致性序列，只是对与原参考序列不一致的位点进行替换，其他位点不处理，这样生成的一致性序列与参考序列差别较小，也不包含 N 碱基。另外一种方法类似，比对之后将每个位点出现频率最高的碱基提取出来，没有比对的位置选择 N 碱基替换，最终结果长度会与参考序列不同，里面也有可能包含很多 N 碱基。目前采用的软件有 bcftools，ivar，medaka 等方法。下面分别来试一下，然后分别比较一些结果。

1. samtools mpileup -A -d 1000 -B -Q 0 --reference MN908947.3.fa
   
2. ncov.primertrim.sorted.bam | ivar consensus -p ivar_consensus -n N

复制代码

2.4 统计长度
      分别比较参考序列，ivar软件生成的一致性序列。

1. $ seqkit stat MN908947.3.fa ivar_consensus.fa

复制代码

1. blat MN908947.3.fa ivar_consensus.fa -out=axt blat.axt

复制代码

     生成blast格式直接看

1. blat MN908947.3.fa ivar_consensus.fa -out=blast blat.out

复制代码