新冠病毒快速鉴定

生信喵

背景
      目前新冠病毒的鉴定可以采用抗体抗原反应的快速鉴定，荧光定量 PCR 以及宏基因组测序等方法。这里我们主要介绍宏基因组测序的方法如何来鉴定新冠病毒。该方法无需扩增，通过测序的方法直接测序新冠病毒序列，可以得到全基因组序列，准确性更高。但该方法受限于成本，目前主要用于科学研究中。
      对于宏基因组测序的新冠病毒样本，可以直接进行物种分类鉴定。将测序数据与物种分类数据库直接比对即可。这里选择使用 centrifuge 工具进行物种分类鉴定。
      
      宏基因组测序分析新冠流程图
      Centrifuge 是一款快速有效的宏基因组物种组成分类的软件，采用了结合 BWT 变换（Burrows-Wheeler transform，BWT）和 FM 索引（Ferragina-Manzini ，FM）的策略对序列分类进行优化，通过基因组压缩策略有效降低了内存的需求，因此可以处理 NT 库级别的库索引。

1 软件安装
      Centrifuge 允许一条序列可以有多个 taxonomy 标签，并允许通过设置阈值将多个 hits 回归到 LCA 模式，针对 multi-hit 模式，通过 EM 算法可以进行丰度定量。centrifuge-kreport 提供了将 Centrifuge 的结果转换成 Kraken 风格的结果。
      官网：http://www.ccb.jhu.edu/software/centrifuge
      github 主页：https://github.com/infphilo/centrifuge

1. #下载解压缩软件
   
2. wget ftp://ftp.ccb.jhu.edu/pub/infphilo/centrifuge/downloads/centrifuge-1.0.3-beta-Linux_x86_64.zip
   
3. unzip centrifuge-1.0.3-beta-Linux_x86_64.zip

复制代码

2 下载数据库索引
      centrifuge 的数据库建库比较麻烦，所以可以选择一些公共数据库。对于新冠病毒鉴定，官方提供了一个 h+v+c 的数据库，该数据库里面包含了人基因组+病毒基因组+106 株SARS-CoV-2 基因组
      下载新冠病毒比对数据库：

1. mkdir centrifuge_h+v_v2_20200327/
   
2. wget -O h+v+c.tar.gz https://zenodo.org/record/3732127/files/h+v+c.tar.gz?download=1
   
3. tar -zxvf h+v+c.tar.gz

复制代码

3 物种分类鉴定
      centrifuge 的使用非常简单，输入数据包含测序的数据以及索引文件。可支持二代和三代测序数据，输入为 fastq 格式文件即可，也支持 fasta 格式以及原始 qseq 格式文件，同时支持 pairend 数据，也支持压缩格式，其中索引只写前缀名即可。如果有多个样本，可以进行批量操作。

1. #centrifuge 进行物种分类鉴定
   
2. <blockquote>centrifuge -x /share/home/xiehs/04.ncov/28.data/database/centrifuge_h+v_v2_20200327/hvc -1 /share/home/xiehs/04.ncov/28.data/sra/SRR10971381.sra_1.fastq.gz -2 /share/home/xiehs/04.ncov/28.data/sra/SRR10971381.sra_2.fastq.gz  --report-file report.tsv -S result.tsv -p 12 >centrifuge.log

复制代码

     选项参数：
      -x：索引，写前缀名即可，将索引添加到变量 CENTRIFUGE_INDEXES，就可以省去每次填写索引部分。
      -1：pairend 测序 reads1
      -2：pairend 测序 reads2
      -U:单端测序，若果有多个数据，用逗号分隔
      -S：输出结果文件，按照测序 reads 给出报告，
      --report-file ：输出结果文件，按照比对上物种报告
      -p：多线程，默认为 1
      --out-fmt：输出文件格式，列表格式或者 sam 格式
      --min-hitlen：最小比对长度，默认 22，最小 15
      -k：处理多重比对，如果一条 reads 比对到多个位置，如果小于这个值则保留，大于这个值，将合并所有比对到一个更高的物种分类层级。默认为 5。
      -u：一个整数，达到这个行数停止
      -s：跳过前面固定行
      -5：跳过 5‘端多少 bp
      -3：跳过 3‘端多少 bp
      --exclude-taxids：排除掉某个物种，例如人类宿主 9606

4 结果解读
      centrifuge 默认会输出两个文件，分别是按照 reads 进行统计的结果与按照物种进行统计的结果。
      1、按照 reads 进行统计的结果 result.tsv
      
      该文件一共分类 8 列 。
      1：原始 read ID ;
      2、比对到数据库中的序列 ID，如果使用的是 Refseq 数据库或者 nt 库，则是序列的Accession ID;
      3、物种分类 ID，第二列比对上序列对应的物种分类 ID;
      4、classification 的分值，比对上的序列之和;
      5、第二好比对结果分值;
      6、比对到序列的长度；
      7、比对的 reads 长度；
      8、这条 reads 比对上多少个物种序列；

      2、按照比对上的物种进行统计 report.tsv
      
      1、比对上物种名字，如果鉴定不到种，则上升一级；
      2、物种分类 ID；
      3、物种分类层级 rank；
      4、对应基因组大小；
      5、比对到的 reads 数目，包括多重比对的结果；
      6、唯一比对上的 reads 数目；
      7、比对的丰度，比对上区域/基因组长度。

5 过滤结果
      由于序列相似性的缘故，一条序列可能会比对到数据库中多个物种，Centrifuge 原始的结果会鉴定到很多物种，这就需要对原始数据进行过滤，通常选择每条序列最优的比对。然后根据每个物种比对上的 reads 数进行过滤，同时也可以根据鉴定到的物种水平进行筛选。可以根据比对上的 reads 数目或者唯一比对的 reads 数据进行过滤，过滤后的结果可以导入 megan 进行数据可视化。

1. awk -F "\t" '{if ($3=="species" && $6 >50) print $1"\t"$6}' report.tsv >filter.tsv

复制代码

6 pavian 结果可视化
      pavian 是一款基于 shinny 的 R 包，可以生成交互式的网页结果。也可以使用在线版本的pavian。支持 kraken，metaphlan 格式结果。如果要利用 pavian 可视化 centrifuge 结果，
      需要首先将其转换为 kraken 格式结果。
      github 主页：https://github.com/fbreitwieser/pavian
      在线工具：https://fbreitwieser.shinyapps.io/pavian/
centrifuge-kreport 结果转换：
      对 centrifuge 生成的 reads 统计的结果 report.tsv 进行重新计算。

1. centrifuge-kreport -x centrifuge_h+p+v_20200318/hvc report.tsv >report_kraken.tsv

复制代码

     统计结果：
      

桑基图：
      桑基图（Sankey diagram），即桑基能量分流图，也叫桑基能量平衡图。它是一种特定类型的流程图，图中延伸的分支的宽度对应数据流量的大小，通常应用于能源、材料成分、金融等数据的可视化分析。因 1898 年 Matthew Henry Phineas Riall Sankey 绘制的“蒸汽机的能源效率图”而闻名，此后便以其名字命名为“桑基图”。（百度百科）
      桑基图主要有边、流量和支点组成，其中边代表了流动的数据，流量代表了流动数据的具体数值，节点代表了不同的分类，边的宽度与流量成比例的显示，边越宽，数值越大。
      
       pavian 桑基图