新冠参考基因组构建

生信喵

背景
      目前新冠病毒的基因组拼接主要采用与参考序列比对，生成一致性序列的方法。所以，参考序列就非常重要，那么参考序列从何而来，参考序列是否准备，遇到新物种如何构建参考序列？
      目前普遍使用的新冠病毒参考序列为 NC_045512.2，该序列为 2020 年 1 月 18 日第一株公布出来的新型冠状病毒序列。样品来自武汉采集样本，原始 GenBank accession number 为MN908947，refseq 库 accession number 为 NC_045512.2，长度 29903bp。该序列发表在文章《A new coronavirus associated with human respiratory disease in China》上。
      
      本节内容，我们将重新拼接该序列，并对该拼接结果进行验证。
1 原始数据处理
      目前该数据存储在项目 PRJNA603194 中，我们将下载该数据。
      
      测序数据样本来自于一个 41 岁男性患者，测序为宏基因组测序，里面除了包含新冠病毒序列，宿主人基因组之外，还包括其他一些微生物。
      数据 Accession ID 为 SRR10971381。测序平台为 MiniSeq，双末端 pairend 150bp 测序，数据量为 8G。
1.1 下载数据
      数据地址：https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA603194

1. #下载序列，由于 prefetch 速度较慢，我们直接使用 wget 进行下载
   
2. wget -O SRR10971381.sra https://storage.googleapis.com/nih-sequence-readarchive/run/SRR10971381/SRR10971381
   
3. #转换数据
   
4. fasterq-dump SRR10971381.sra

复制代码

1.2 质控过滤
      利用 fastqc 软件对原始数据进行质控过滤。

1. fastqc -f fastq SRR10971381_1.fastq SRR10971381_2.fastq

复制代码

     质控之后，利用 fastp 软件对数据过滤。

1. fastp -i SRR10971381_1.fastq -I SRR10971381_2.fastq -o clean.1.fq.gz -O
   
2. clean.2.fq.gz -z 4 -q 20 -u 30 -n 10 -h clean.html
   

复制代码

     过滤之后在利用 fastqc 进行质控。

1. fastqc -f fastq clean.1.fq.gz clean.2.fq.gz

复制代码

2 去除宿主
      由于测序数据中宿主人的细胞占据绝大部分比例，因此需要去除宿主。将测序数据与人基因组序列进行比对，比对之后，如果能够比对上则为人基因组序列。我们挑选比对不上参考序列的部分。也可以先拼接全部数据，在对拼接结果进行过滤，两种方法均可。
1、下载参考序列，建立索引

1. #下载人基因组序列建立索引
   
2. ~/.aspera/connect/bin/ascp -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh -
   
3. -overwrite=diff -QTr -l6000m
   
4. anonftp@ftp.ncbi.nlm.nih.gov:1000genomes/ftp/technical/reference/human_g1k_v37
   
5. .fasta.gz ./
   
6. gunzip -zxvf human_g1k_v37.fasta.gz
   
7. nohup time bwa-mem2 index human_g1k_v37.fasta &

复制代码

2、将测序数据与宿主基因组进行比对，过滤宿主，可以使用-f 4，但是使用该方法，对于pairend 序列，后面得到 ID 位置会不匹配。这里使用-G 2，过滤掉双末端比对上的，虽然这样可能会继续保留单端比对上的部分宿主序列，但可以保证 pairend reads 成对出现，方便下一步处理。

1. #bwa-mem2 加快比对速度
   
2. bwa-mem2 mem -t 24
   
3. /share/home/xiehs/04.ncov/28.data/reference/human_g1k_v37_bwamem2/human_g1k_v37.fasta
   
4. clean.1.fq.gz clean.2.fq.gz >ncov.sam
   
5. #过滤宿主序列
   
6. samtools fastq -G 2 ncov.sam -1 ncov.1.fq -2 ncov.2.fq
   
7. gzip ncov.1.fq ncov.2.fq
   
8. #统计过滤前后数变化
   
9. seqkit stat clean.1.fq.gz clean.2.fq.gz
   
10. seqkit stat ncov.1.fq.gz ncov.2.fq.gz
    

复制代码

3 宏基因组拼接
      数据质控，过滤，去除宿主之后就可以开始拼接了，由于去除了绝大部分宿主基因组，剩余数据量明显减少，极大降低了拼接的难度和资源消耗。但有些病毒序列也可以比对到人基因组，先去宿主可能影响基因组拼接，可以尝试先拼接，在去宿主的方式。比较二者之间的区别。
3.1 megahit 拼接
      二代宏基因组测序可以使用 megahit，spades，SOAPdenovo 等软件进行拼接。不过由于测序读长短，物种组成复杂，拼接难度较大。本案例中，想要直接拼接得到新冠病毒的全基因组序列。根据以往经验，病毒基因组由于序列读长短，深度高，杂合度高，很难拼接出来，即使是病毒纯培养测序微生物，也很难拼接出完整基因组。不过该案例中只通过一次拼接，就拼出了一条完整的新冠病毒序列。

1. megahit -t 24 -o megahit/ -1 ncov.1.fq.gz -2 ncov.2.fq.gz 1>megahit.log 2>megahit.err

复制代码

     选项参数
      -1：reads1
      -2：reads2
      -o：数据文件夹，自动创建
      -h 显示参数详细
      --k-min 27 --k-max 191 --k-step 20 # 手动设置 kmer
      -r 单端
      -t 设置线程数，默认全用
      --use-gpu 支持 GPU 运算
      --continue 支持中断继续运行
      结果解读：其中 final.contigs.fa 为拼接得到的基因组序列。


3.2 拼接结果评估
1 、统计结果
      对拼接统计结果进行统计，一共拼接了 14946 条序列，其中最长序列长度为 29863bp，最短序列为 200bp，其中最长片段接近于冠状病毒长度，如果该条片段为冠状病毒则非常好，可能得到病毒完整基因组序列。不过目前只是潜在可能，还需做进一步验证。
      
      统计每一条长度

1. seqkit fx2tab final.contigs.fa |awk -F"\t" '{print $1,length($2)}'

复制代码

2、NCBI Blast 比对
      将最长一条序列与 NCBI blast 数据库进行比对，验证序列来自于何种物种，这里注意，由于当前的最新的数据库已经包含新冠病毒数据库，因此，需要比对 2019 年 12 月之前的数据库，或者直接与 SARS 序列进行比对。比对结果显示，序列可以很好的比对到冠状病毒序列上。

1. blastn -query longest.fa -out blast.out -db /ifs1/MetaDatabase/nt_20200716/nt -outfmt 6 -evalue 1e-5 -num_threads 24

复制代码

     
      blast 比对结果（部分）

生信喵

3、筛选最长序列
      将最长一条序列挑选出来，进行验证，如果能够比对到新冠病毒序列，则该序列为冠状病毒。
#排序筛选最长序列

1. seqkit sort -l -r final.contigs.fa | seqkit seq -w 0 | head -2 >longest.fa

复制代码

#统计长度

1. seqkit stat longest.fa

复制代码

#统计结果

1. file format type num_seqs sum_len min_len avg_len max_len
   
2. longest.fa FASTA DNA 1 29,863 29,863 29,863 29,863

复制代码

#修改名字

1. mv longest.fa ncov.fa

复制代码

4 验证结果
      接下来我们将验证结果的准确性和完整性，主要包括该序列是否完整，29863bp 是否已经是完整基因组，拼接结果中是否还包含其他冠状病毒序列，序列是否连接错误，测序位点是否正确？
4.1 与已发表 NC_045512 进行比对
      比较拼接得到的 ncov.fa 与已发表 NC_045512 之间的差别，本次拼接结果比已发表的NC_045512 少 40bp，这 40bp 刚好在序列的头尾部分。

1. $ seqkit stat NC_045512.fa ncov.fa
   
2. file format type num_seqs sum_len min_len avg_len max_len
   
3. NC_045512.fa FASTA DNA 1 29,903 29,903 29,903 29,903
   
4. ncov.fa FASTA DNA 1 29,863 29,863 29,863 29,863

复制代码

     最后将重新拼接得到 ncov.fa 与已发表序列 NC_045512 进行比较，二者相似性达到 100%。其中 NC_045512 尾部 39bp 没有比对上 ncov。

1. $ lastz NC_045512.fa ncov.fa --out=lastz.axt
   
2. #比对结果局部
   
3. 0 NC_045512.2 1 29862 k141_5526 2 29863 + 2819610
   
4. ATTAAAGGTTTATACCTTCCCAGGTAACAAACCAACCAACTTTCGATCTCTTGTAGATCTGTTCTCTAAACGAACTTT
   
5. AAAATTAAAGGTTTATACCTTCCCAGGTAACAAACCAACCAACTTTCGATCTCTTGTAGATCTGTTCTCTAAACGAAC
   
6. TTTAAA

复制代码

     通过直接观察两个比对结果首位序列，我们发现以下两点：
      第一：NC_045512 头部比 ncov.fa 少了一个碱基 G；
      第二：NC_045512 尾部多了 GAATGACAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA。

1. $ seqkit subseq -r 1:100 NC_045512.fa ncov.fa
   
2. >NC_045512.2 Severe acute respiratory syndrome coronavirus 2 isolate Wuhan-Hu-1, complete genom
   
3. ATTAAAGGTTTATACCTTCCCAGGTAACAAACCAACCAACTTTCGATCTCTTGTAGATCT
   
4. GTTCTCTAAACGAACTTTAAAATCTGTGTGGCTGTCACTC
   
5. >k141_5526 flag=1 multi=191.8580 len=29863
   
6. <font color="#ff0000">G</font>ATTAAAGGTTTATACCTTCCCAGGTAACAAACCAACCAACTTTCGATCTCTTGTAGATC
   
7. TGTTCTCTAAACGAACTTTAAAATCTGTGTGGCTGTCACT
   
8. (base) meta 15:13:43 /ifs1/User/meta/xinguan/sra
   
9. $ seqkit subseq -r -100:-1 NC_045512.fa ncov.fa
   
10. >NC_045512.2 Severe acute respiratory syndrome coronavirus 2 isolate Wuhan-Hu1, complete genome
    
11. TAATGTGTAAAATTAATTTTAGTAGTGCTATCCCCATGTGATTTTAATAGCTTCTTAGG<font color="#ff0000">A
    
12. GAATGACAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA</font>
    
13. >k141_5526 flag=1 multi=191.8580 len=29863
    
14. AGTGAACAATGCTAGGGAGAGCTGCCTATATGGAAGAGCCCTAATGTGTAAAATTAATTT
    
15. TAGTAGTGCTATCCCCATGTGATTTTAATAGCTTCTTAGG

复制代码

4.2 PCR 产物验证
      得到全基因组序列之后，可以进行 PCR 扩增，如果能够扩增出全基因组，则证明基因组没有问题。PCR 扩增产物可以覆盖全部基因组。对于尾部 A 的问题，如果 PCR 产物可以比对到 NC_045512 尾部 A，则可以证明基因组中包含该段区域。

1. #PCR 产物验证
   
2. #下载 PCR 扩增测序数据
   
3. prefetch SRR14867660 -O ./
   
4. locate SRR14867660.sra
   
5. fasterq-dump SRR14867660.sra
   
6. #与拼接结果比对
   
7. bwa-mem2 index ncov.fa
   
8. bwa-mem2 mem -t 12 ncov.fa SRR14867660_1.fastq SRR14867660_2.fastq >pcr.sam
   
9. #samtools 处理比对结果
   
10. samtools sort -@ 12 -o pcr.sorted.bam pcr.sam
    
11. samtools index pcr.sorted.bam
    
12. #对 ncov.fa 建立索引
    
13. samtools faidx ncov.fa
    
14. #tablet 可视化
    
15. #将文件拷贝至 windows 下使用 tablet 可视化
    
16. mkdir tablet
    
17. mv ncov.fa ncov.fa.fai pcr.sorted.bam pcr.sorted.bam.bai tablet
    

复制代码

     
       PCR 扩增测序序列 tablet 可视化
      与 NC_045512 参考序列比对，尾部有 4 条 reads 覆盖到，基因组中包含该区域。

      
      PCR 扩增序列可以比对到参考序列尾部 A 尾巴