新冠参考基因组构建
本帖最后由 生信喵 于 2022-3-14 21:36 编辑背景
目前新冠病毒的基因组拼接主要采用与参考序列比对，生成一致性序列的方法。所以，参考序列就非常重要，那么参考序列从何而来，参考序列是否准备，遇到新物种如何构建参考序列？
目前普遍使用的新冠病毒参考序列为 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
#下载序列，由于 prefetch 速度较慢，我们直接使用 wget 进行下载
wget -O SRR10971381.sra https://storage.googleapis.com/nih-sequence-readarchive/run/SRR10971381/SRR10971381
#转换数据
fasterq-dump SRR10971381.sra
1.2 质控过滤
利用 fastqc 软件对原始数据进行质控过滤。
fastqc -f fastq SRR10971381_1.fastq SRR10971381_2.fastq 质控之后，利用 fastp 软件对数据过滤。
fastp -i SRR10971381_1.fastq -I SRR10971381_2.fastq -o clean.1.fq.gz -O
clean.2.fq.gz -z 4 -q 20 -u 30 -n 10 -h clean.html
过滤之后在利用 fastqc 进行质控。
fastqc -f fastq clean.1.fq.gz clean.2.fq.gz
2 去除宿主
由于测序数据中宿主人的细胞占据绝大部分比例，因此需要去除宿主。将测序数据与人基因组序列进行比对，比对之后，如果能够比对上则为人基因组序列。我们挑选比对不上参考序列的部分。也可以先拼接全部数据，在对拼接结果进行过滤，两种方法均可。
1、下载参考序列，建立索引
#下载人基因组序列建立索引
~/.aspera/connect/bin/ascp -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh -
-overwrite=diff -QTr -l6000m
anonftp@ftp.ncbi.nlm.nih.gov:1000genomes/ftp/technical/reference/human_g1k_v37
.fasta.gz ./
gunzip -zxvf human_g1k_v37.fasta.gz
nohup time bwa-mem2 index human_g1k_v37.fasta &
2、将测序数据与宿主基因组进行比对，过滤宿主，可以使用-f 4，但是使用该方法，对于pairend 序列，后面得到 ID 位置会不匹配。这里使用-G 2，过滤掉双末端比对上的，虽然这样可能会继续保留单端比对上的部分宿主序列，但可以保证 pairend reads 成对出现，方便下一步处理。
#bwa-mem2 加快比对速度
bwa-mem2 mem -t 24
/share/home/xiehs/04.ncov/28.data/reference/human_g1k_v37_bwamem2/human_g1k_v37.fasta
clean.1.fq.gz clean.2.fq.gz >ncov.sam
#过滤宿主序列
samtools fastq -G 2 ncov.sam -1 ncov.1.fq -2 ncov.2.fq
gzip ncov.1.fq ncov.2.fq
#统计过滤前后数变化
seqkit stat clean.1.fq.gz clean.2.fq.gz
seqkit stat ncov.1.fq.gz ncov.2.fq.gz
3 宏基因组拼接
数据质控，过滤，去除宿主之后就可以开始拼接了，由于去除了绝大部分宿主基因组，剩余数据量明显减少，极大降低了拼接的难度和资源消耗。但有些病毒序列也可以比对到人基因组，先去宿主可能影响基因组拼接，可以尝试先拼接，在去宿主的方式。比较二者之间的区别。
3.1 megahit 拼接
二代宏基因组测序可以使用 megahit，spades，SOAPdenovo 等软件进行拼接。不过由于测序读长短，物种组成复杂，拼接难度较大。本案例中，想要直接拼接得到新冠病毒的全基因组序列。根据以往经验，病毒基因组由于序列读长短，深度高，杂合度高，很难拼接出来，即使是病毒纯培养测序微生物，也很难拼接出完整基因组。不过该案例中只通过一次拼接，就拼出了一条完整的新冠病毒序列。
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，其中最长片段接近于冠状病毒长度，如果该条片段为冠状病毒则非常好，可能得到病毒完整基因组序列。不过目前只是潜在可能，还需做进一步验证。
统计每一条长度
seqkit fx2tab final.contigs.fa |awk -F"\t" '{print $1,length($2)}'
2、NCBI Blast 比对
将最长一条序列与 NCBI blast 数据库进行比对，验证序列来自于何种物种，这里注意，由于当前的最新的数据库已经包含新冠病毒数据库，因此，需要比对 2019 年 12 月之前的数据库，或者直接与 SARS 序列进行比对。比对结果显示，序列可以很好的比对到冠状病毒序列上。
blastn -query longest.fa -out blast.out -db /ifs1/MetaDatabase/nt_20200716/nt -outfmt 6 -evalue 1e-5 -num_threads 24
blast 比对结果（部分）
3、筛选最长序列
将最长一条序列挑选出来，进行验证，如果能够比对到新冠病毒序列，则该序列为冠状病毒。
#排序筛选最长序列
seqkit sort -l -r final.contigs.fa | seqkit seq -w 0 | head -2 >longest.fa#统计长度
seqkit stat longest.fa#统计结果file format type num_seqs sum_len min_len avg_len max_len
longest.fa FASTA DNA 1 29,863 29,863 29,863 29,863#修改名字
mv longest.fa ncov.fa
4 验证结果
接下来我们将验证结果的准确性和完整性，主要包括该序列是否完整，29863bp 是否已经是完整基因组，拼接结果中是否还包含其他冠状病毒序列，序列是否连接错误，测序位点是否正确？
4.1 与已发表 NC_045512 进行比对
比较拼接得到的 ncov.fa 与已发表 NC_045512 之间的差别，本次拼接结果比已发表的NC_045512 少 40bp，这 40bp 刚好在序列的头尾部分。
$ seqkit stat NC_045512.fa ncov.fa
file format type num_seqs sum_len min_len avg_len max_len
NC_045512.fa FASTA DNA 1 29,903 29,903 29,903 29,903
ncov.fa FASTA DNA 1 29,863 29,863 29,863 29,863 最后将重新拼接得到 ncov.fa 与已发表序列 NC_045512 进行比较，二者相似性达到 100%。其中 NC_045512 尾部 39bp 没有比对上 ncov。
$ lastz NC_045512.fa ncov.fa --out=lastz.axt
#比对结果局部
0 NC_045512.2 1 29862 k141_5526 2 29863 + 2819610
ATTAAAGGTTTATACCTTCCCAGGTAACAAACCAACCAACTTTCGATCTCTTGTAGATCTGTTCTCTAAACGAACTTT
AAAATTAAAGGTTTATACCTTCCCAGGTAACAAACCAACCAACTTTCGATCTCTTGTAGATCTGTTCTCTAAACGAAC
TTTAAA 通过直接观察两个比对结果首位序列，我们发现以下两点：
第一：NC_045512 头部比 ncov.fa 少了一个碱基 G；
第二：NC_045512 尾部多了 GAATGACAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA。
$ seqkit subseq -r 1:100 NC_045512.fa ncov.fa
>NC_045512.2 Severe acute respiratory syndrome coronavirus 2 isolate Wuhan-Hu-1, complete genom
ATTAAAGGTTTATACCTTCCCAGGTAACAAACCAACCAACTTTCGATCTCTTGTAGATCT
GTTCTCTAAACGAACTTTAAAATCTGTGTGGCTGTCACTC
>k141_5526 flag=1 multi=191.8580 len=29863
<font color="#ff0000">G</font>ATTAAAGGTTTATACCTTCCCAGGTAACAAACCAACCAACTTTCGATCTCTTGTAGATC
TGTTCTCTAAACGAACTTTAAAATCTGTGTGGCTGTCACT
(base) meta 15:13:43 /ifs1/User/meta/xinguan/sra
$ seqkit subseq -r -100:-1 NC_045512.fa ncov.fa
>NC_045512.2 Severe acute respiratory syndrome coronavirus 2 isolate Wuhan-Hu1, complete genome
TAATGTGTAAAATTAATTTTAGTAGTGCTATCCCCATGTGATTTTAATAGCTTCTTAGG<font color="#ff0000">A
GAATGACAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA</font>
>k141_5526 flag=1 multi=191.8580 len=29863
AGTGAACAATGCTAGGGAGAGCTGCCTATATGGAAGAGCCCTAATGTGTAAAATTAATTT
TAGTAGTGCTATCCCCATGTGATTTTAATAGCTTCTTAGG
4.2 PCR 产物验证
得到全基因组序列之后，可以进行 PCR 扩增，如果能够扩增出全基因组，则证明基因组没有问题。PCR 扩增产物可以覆盖全部基因组。对于尾部 A 的问题，如果 PCR 产物可以比对到 NC_045512 尾部 A，则可以证明基因组中包含该段区域。
#PCR 产物验证
#下载 PCR 扩增测序数据
prefetch SRR14867660 -O ./
locate SRR14867660.sra
fasterq-dump SRR14867660.sra
#与拼接结果比对
bwa-mem2 index ncov.fa
bwa-mem2 mem -t 12 ncov.fa SRR14867660_1.fastq SRR14867660_2.fastq >pcr.sam
#samtools 处理比对结果
samtools sort -@ 12 -o pcr.sorted.bam pcr.sam
samtools index pcr.sorted.bam
#对 ncov.fa 建立索引
samtools faidx ncov.fa
#tablet 可视化
#将文件拷贝至 windows 下使用 tablet 可视化
mkdir tablet
mv ncov.fa ncov.fa.fai pcr.sorted.bam pcr.sorted.bam.bai tablet
PCR 扩增测序序列 tablet 可视化
与 NC_045512 参考序列比对，尾部有 4 条 reads 覆盖到，基因组中包含该区域。
PCR 扩增序列可以比对到参考序列尾部 A 尾巴
页:
[1]