新冠病毒数据分析
-
nCoV-2019 novel coronavirus bioinformatics protocol
https://artic.network/ncov-2019/ncov2019-bioinformatics-sop.html
1.硬件配置
ARTIC 针对的设备是MinION，推荐的配置为：
Intel i7 or Xeon processor、16GB RAM、1TB SSD hard drive、USB 3
这个是一个主流笔记本的配置2.软件环境
2.1 操作系统
64-bit UNIX, Linux or similar environment
Mac OS X (Yosemite or later), Linux (e.g., Ubuntu 16 or later), or Windows 10 Subsystem for Linux。2.2 应用环境
https://conda.io/docs/user-guide/install/
64-bit Python 3.6 version of Miniconda2.3 安装artic
git clone --recursive https://github.com/artic-network/artic-ncov2019.git
conda env create -f artic-ncov2019/environment.yml
conda activate artic-ncov2019
去激活：conda deactivate 删除：conda remove --name artic-ncov2019 --all更新artic
cd artic-ncov2019
conda env remove -n artic-ncov2019
conda env create -f environment.yml2.4 开始分析
创建分析目录
mkdir analysis
cd analysis
mkdir run_name
cd run_name
2.4.1 激活环境
source activate artic-ncov2019
2.4.2 basecalling
快速模式
guppy_basecaller -c dna_r9.4.1_450bps_fast.cfg -i /path/to/reads -s run_name -x auto -r
高精度模式
guppy_basecaller -c dna_r9.4.1_450bps_hac.cfg -i /path/to/reads -s run_name -x auto -r
2.4.3 Demultiplexing
guppy_barcoder --require_barcodes_both_ends -i run_name -s output_directory --arrangements_files "barcode_arrs_nb12.cfg barcode_arrs_nb24.cfg"
--require_barcodes_both_ends:Reads will only be classified if there is a barcode above the min_score at both ends of the read.
-i Path to input fastq files
-s Path to save fastq files.
--arrangements_files Files containing arrangements.2.4.4 过滤
由于有嵌合体，需要按照长度进行过滤。
artic guppyplex --min-length 400 --max-length 700 --directory output_directory/barcode03 --prefix run_name2.4.5 分析
artic minion --normalise 200 --threads 4 --scheme-directory ~/artic-ncov2019/primer_schemes --read-file run_name_barcode03.fastq --fast5-directory path_to_fast5 --sequencing-summary path_to_sequencing_summary.txt nCoV-2019/V3 samplename
产生如下文件：
samplename.rg.primertrimmed.bam - BAM file for visualisation after primer-binding site trimming
samplename.trimmed.bam - BAM file with the primers left on (used in variant calling)
samplename.merged.vcf - all detected variants in VCF format
samplename.pass.vcf - detected variants in VCF format passing quality filter
samplename.fail.vcf - detected variants in VCF format failing quality filter
samplename.primers.vcf - detected variants falling in primer-binding regions
samplename.consensus.fasta - consensus sequence -
minimap2 -R '@RG\tID:sras19\tPU:nanopore\tSM:sars19\tPL:minion\tLB:sqklsk109' -a -x map-ont -t 24 ~/ref/h1n1/ref.fasta barcode$i.fastq | samtools view -bS -F 4 - | samtools sort -o barcode$i.bam -
samtools index barcode$i.bam
medaka consensus barcode$i.bam barcode$i.hdf
medaka variant ~/ref/h1n1/ref.fasta barcode$i.hdf barcode$i.vcfbgzip barcode$i.vcf
tabix -p vcf barcode$i.vcf.gzlongshot -P 0 -F -A --no_haps --bam barcode$i.bam --ref ~/ref/h1n1/ref.fasta --out barcode$i.longshot.vcf --potential_variants barcode$i.vcf.gz
artic_vcf_filter --longshot barcode$i.longshot.vcf barcode$i.pass.vcf barcode$i.fail.vcf
artic_make_depth_mask ~/ref/h1n1/ref.fasta barcode$i.bam barcode$i.coverage_mask.txtbgzip -f barcode$i.pass.vcf
tabix -p vcf barcode$i.pass.vcf.gzartic_mask ~/ref/h1n1/ref.fasta barcode$i.coverage_mask.txt barcode$i.fail.vcf barcode$i.preconsensus.fasta
bcftools consensus -f barcode$i.preconsensus.fasta barcode$i.pass.vcf.gz -m barcode$i.coverage_mask.txt -o barcode$i.consensus.fastahttps://github.com/nanoporetech/medaka/issues/149
longshot只是做了一些统计https://artic.readthedocs.io/en/latest/release-notes/
Longshot added to Medaka to permit filter VCF on depth
Longshot似乎只是做了一些过滤https://bedtools.readthedocs.io/en/latest/content/tools/maskfasta.html
咱们的过程没有用到 medaka tools annotate 步骤 这个步骤和longshot是二选一的过程
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7480024/ -
-
(artic-ncov2019) bioinfo@genostack:~/artic_ncov_2019/analysis/run_name$ artic minion --scheme-directory ~/artic_ncov_2019/artic-ncov2019/primer_schemes/ nCoV-2019/V3 nCov19 --read-file guppy_fastq/all_sars19.fastq --fast5-directory /mnt/sdf1/sars_2019/ --sequencing-summary guppy_fastq/sequencing_summary.txt
Running: nanopolish index -s guppy_fastq/sequencing_summary.txt -d /mnt/sdf1/sars_2019/ guppy_fastq/all_sars19.fastq
[readdb] indexing /mnt/sdf1/sars_2019/
[readdb] num reads: 2072359, num reads with path to fast5: 2072359
Running: minimap2 -a -x map-ont -t 8 /home/bioinfo/artic_ncov_2019/artic-ncov2019/primer_schemes//nCoV-2019/V3/nCoV-2019.reference.fasta guppy_fastq/all_sars19.fastq | samtools view -bS -F 4 - | samtools sort -o nCov19.sorted.bam -
[M::mm_idx_gen::0.0071.36] collected minimizers
[M::mm_idx_gen::0.0112.78] sorted minimizers
[M::main::0.0122.75] loaded/built the index for 1 target sequence(s)
[M::mm_mapopt_update::0.0152.37] mid_occ = 3
[M::mm_idx_stat] kmer size: 15; skip: 10; is_hpc: 0; #seq: 1
[M::mm_idx_stat::0.0172.15] distinct minimizers: 5587 (99.93% are singletons); average occurrences: 1.004; average spacing: 5.332
[M::worker_pipeline::80.0646.71] mapped 550538 sequences
[M::worker_pipeline::152.0505.23] mapped 548324 sequences
[M::worker_pipeline::209.0364.73] mapped 548462 sequences
[M::worker_pipeline::259.8403.83] mapped 425035 sequences
[M::main] Version: 2.17-r941
[M::main] CMD: minimap2 -a -x map-ont -t 8 /home/bioinfo/artic_ncov_2019/artic-ncov2019/primer_schemes//nCoV-2019/V3/nCoV-2019.reference.fasta guppy_fastq/all_sars19.fastq
[M::main] Real time: 259.915 sec; CPU: 995.230 sec; Peak RSS: 2.570 GB
[bam_sort_core] merging from 3 files and 1 in-memory blocks...
Running: samtools index nCov19.sorted.bam
Running: align_trim --start --normalise 100 /home/bioinfo/artic_ncov_2019/artic-ncov2019/primer_schemes//nCoV-2019/V3/nCoV-2019.scheme.bed --report nCov19.alignreport.txt < nCov19.sorted.bam 2> nCov19.alignreport.er | samtools sort -T nCov19 - -o nCov19.trimmed.rg.sorted.bam
Running: align_trim --normalise 100 /home/bioinfo/artic_ncov_2019/artic-ncov2019/primer_schemes//nCoV-2019/V3/nCoV-2019.scheme.bed --remove-incorrect-pairs --report nCov19.alignreport.txt < nCov19.sorted.bam 2> nCov19.alignreport.er | samtools sort -T nCov19 - -o nCov19.primertrimmed.rg.sorted.bam
Running: samtools index nCov19.trimmed.rg.sorted.bam
Running: samtools index nCov19.primertrimmed.rg.sorted.bam
Running: nanopolish variants --min-flanking-sequence 10 -x 1000000 --progress -t 8 --reads guppy_fastq/all_sars19.fastq -o nCov19.nCoV-2019_1.vcf -b nCov19.trimmed.rg.sorted.bam -g /home/bioinfo/artic_ncov_2019/artic-ncov2019/primer_schemes//nCoV-2019/V3/nCoV-2019.reference.fasta -w "MN908947.3:1-29904" --ploidy 1 -m 0.15 --read-group nCoV-2019_1
[post-run summary] total reads: 6604, unparseable: 0, qc fail: 26, could not calibrate: 22, no alignment: 102, bad fast5: 0
Running: nanopolish variants --min-flanking-sequence 10 -x 1000000 --progress -t 8 --reads guppy_fastq/all_sars19.fastq -o nCov19.nCoV-2019_2.vcf -b nCov19.trimmed.rg.sorted.bam -g /home/bioinfo/artic_ncov_2019/artic-ncov2019/primer_schemes//nCoV-2019/V3/nCoV-2019.reference.fasta -w "MN908947.3:1-29904" --ploidy 1 -m 0.15 --read-group nCoV-2019_2
[post-run summary] total reads: 5890, unparseable: 0, qc fail: 24, could not calibrate: 15, no alignment: 86, bad fast5: 0
Running: artic_vcf_merge nCov19 /home/bioinfo/artic_ncov_2019/artic-ncov2019/primer_schemes//nCoV-2019/V3/nCoV-2019.scheme.bed nCoV-2019_1:nCov19.nCoV-2019_1.vcf nCoV-2019_2:nCov19.nCoV-2019_2.vcf
Running: artic_vcf_filter --nanopolish nCov19.merged.vcf nCov19.pass.vcf nCov19.fail.vcf
Running: artic_make_depth_mask --store-rg-depths /home/bioinfo/artic_ncov_2019/artic-ncov2019/primer_schemes//nCoV-2019/V3/nCoV-2019.reference.fasta nCov19.primertrimmed.rg.sorted.bam nCov19.coverage_mask.txt
Running: artic_plot_amplicon_depth --primerScheme /home/bioinfo/artic_ncov_2019/artic-ncov2019/primer_schemes//nCoV-2019/V3/nCoV-2019.scheme.bed --sampleID nCov19 --outFilePrefix nCov19 nCov19.depths
Running: bgzip -f nCov19.pass.vcf
Running: tabix -p vcf nCov19.pass.vcf.gz
Running: artic_mask /home/bioinfo/artic_ncov_2019/artic-ncov2019/primer_schemes//nCoV-2019/V3/nCoV-2019.reference.fasta nCov19.coverage_mask.txt nCov19.fail.vcf nCov19.preconsensus.fasta
Running: bcftools consensus -f nCov19.preconsensus.fasta nCov19.pass.vcf.gz -m nCov19.coverage_mask.txt -o nCov19.consensus.fasta
Note: the --sample option not given, applying all records regardless of the genotype
Applied 2 variants
Running: artic_fasta_header nCov19.consensus.fasta "nCov19/ARTIC/nanopolish"
Running: cat nCov19.consensus.fasta /home/bioinfo/artic_ncov_2019/artic-ncov2019/primer_schemes//nCoV-2019/V3/nCoV-2019.reference.fasta > nCov19.muscle.in.fasta
Running: muscle -in nCov19.muscle.in.fasta -out nCov19.muscle.out.fastaMUSCLE v3.8.1551 by Robert C. Edgar
http://www.drive5.com/muscle
This software is donated to the public domain.
Please cite: Edgar, R.C. Nucleic Acids Res 32(5), 1792-97.nCov19.muscle.in 2 seqs, lengths min 29903, max 29903, avg 29903
00:00:00 20 MB(-2%) Iter 1 100.00% K-mer dist pass 1
00:00:00 20 MB(-2%) Iter 1 100.00% K-mer dist pass 2
00:01:16 1004 MB(-116%) Iter 1 100.00% Align node
00:01:16 1004 MB(-116%) Iter 1 100.00% Root alignment详细过程分析
1.nanopolish index -s guppy_fastq/sequencing_summary.txt -d /mnt/sdf1/sars_2019/ guppy_fastq/all_sars19.fastq
nanopolish需要使用原始数据，index命令会将原始数据和guppy产生的reads关联起来2.minimap2 -a -x map-ont -t 8 /home/bioinfo/artic_ncov_2019/artic-ncov2019/primer_schemes//nCoV-2019/V3/nCoV-2019.reference.fasta guppy_fastq/all_sars19.fastq | samtools view -bS -F 4 - | samtools sort -o nCov19.sorted.bam -
samtools index nCov19.sorted.bam
将序列和参考序列进行比对 生成bam文件 此步骤也可以使用bwa mem进行比对3.align_trim --start --normalise 100 /home/bioinfo/artic_ncov_2019/artic-ncov2019/primer_schemes//nCoV-2019/V3/nCoV-2019.scheme.bed --report nCov19.alignreport.txt < nCov19.sorted.bam 2> nCov19.alignreport.er | samtools sort -T nCov19 - -o nCov19.trimmed.rg.sorted.bam
目的：这一步严重依赖primer的相关知识 相关论文和工具可以从下面链接查到
https://primalscheme.com/
The purpose of alignment post-processing is:
assign each read alignment to a derived amplicon
using the derived amplicon, assign each read a read group based on the primer pool
softmask read alignments within their derived amplicon4.突变分析 medaka-longshot 流程 在这个流程中Longshot is used in the ARTIC pipeline simply to annotate the VCF with various statistics (like read support for ALTs). The newest medaka v1.0.0 has a tool to perform this operation:可以被medaka tools annotate --help 取代
medaka consensus
medaka variant or snps (if pipeline has been told not to detect INDELS via --no-indels)
medaka tools annotate (if --no-longshot has been selected)
longshot (if --no-longshot not selected)另外一个突变流程是nanopolish 但是medaka宣称 50X faster than Nanopolish (and can run on GPUs).
5.生成一致性序列
-
https://sci-hub.st/10.3390/cimb43020061
Next-Generation Sequencing (NGS) in COVID-19: A Tool for
SARS-CoV-2 Diagnosis, Monitoring New Strains and
Phylodynamic Modeling in Molecular Epidemiology这个文章里面对武汉新冠当时的情况做了一个介绍 5个样本中国疾控处理、4个样本华大处理的
华大用bwa和hg19人的基因组对比 去除了宿主 然后和NCBI的冠状病毒（具体是哪个序列还不知道）做了对齐 然后使用SPAdes做了一个一致性序列
中国疾控用的是 CLCBio （就是我研究的竞品 CLC workbench）software version 11.0.1 was used for de novo assembly, variant calling, and alignment
有了这些组装的结果之后，就可以做进化分析（ phylogenetic analysis）该文章对covid-19 NGS实验和生信做了综述 里面引用了一些文章还是有价值的：
1.https://pubmed.ncbi.nlm.nih.gov/29154853/
Standards and Guidelines for Validating Next-Generation Sequencing Bioinformatics Pipelines: A Joint Recommendation of the Association for Molecular Pathology and the College of American Pathologists
这个文章里面总结了一个专家共识，提出了17条临床NGS生物信息学管道验证的最佳实践共识建议
属于解释
Terminology
Description
BAM
Compressed (binary) format of a SAM file intended for faster random access (search) of aligned and unaligned sequences and its related metadata. Compression enables smaller file size and storage efficiency which makes it popular over the SAM format. This is a default output of many alignment and post-alignment softwares used in bioinformatics pipeline and commonly used as an input by many variant callers*.
FASTQ
It is a de facto, human readable, file format that stores nucleotide sequences and corresponding quality (PHRED) scores for each nucleotide as an ASCII encoded character.4 This is commonly used for storing unaligned short sequence reads after the steps of base calling and is a typical starting point for NGS bioinformatics pipeline.
PHRED Score
It is a per base (nucleotide) quality score that is defined as an estimated probability for a called base to be incorrect (erroneous call). Mathematically, it is expressed as4where Q is Phred quality score, Pe is the probability for an erroneous base call. The Pe is typically generated by the base calling software which is sequence instrument specific. Therefore, Q values in isolation cannot be used to compare sequence quality across different sequencing platforms.
SAM
Stands for Sequence Alignment/Map format. It is a human readable (text file) file format specification for storing information on aligned sequence. This is a default output of many alignment softwares used in bioinformatics pipeline. Given the large file size and slower random access, BAM format is preferred for routine bioinformatics data processing. This format is helpful for technical troubleshooting when manual review of the stored information is necessary*.
Variant - horizontally complex
When two or more sequence alterations are present on the same read in close proximity such that they may represent a single complex variant. These variants are frequently represented as deletion-insertions and may result in ambiguous sequence description or HGVS nomenclature.
Variant - Left-aligned
If there are multiple potential VCF entries of the same allele length that represent the same variant, then left-alignment refers to the VCF entry with the smallest base position. The base position is typically represented in genomic coordinate for a given primary assembly (eg, GRCh38) and represents the most 5’ position.5
Variant - Normalized
A normalized variant must be parsimonious as well as left-aligned.5
Variant - Parsimony
If there are more than one way to represent the same variant in a VCF file, parsimony refers to the representation with the shortest possible allele length5 (positive and non-zero length).
Variant - vertically complex
A vertically complex variant occurs when three or more alleles are represented by different sequence reads, typically with or uncommonly without a reference (normal) allele, at the same genomic coordinate or set of coordinates.
VCF
Variant Call Format is a versioned, text-file (human readable) specification for storing sequence variant calls. The file contains meta-information containing various details of the variant calling process and definition of headers and format tags, a header line and data lines. Each data line represents a sequence variant defined using a combination of chromosome, position, reference allele, and alternate allele†. -
https://pubmed.ncbi.nlm.nih.gov/31978945/
A Novel Coronavirus from Patients with Pneumonia in China, 2019
这个文章提到了最初武汉的几个样本采用的二代和三代混合测序的方式 -
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7929396/
Bioinformatics resources for SARS-CoV-2 discovery and surveillance
几种实验方法
The workflow of different NGS sequencing approaches currently available for virus discovery and genomic surveillance. The library construction scheme employed in (A) metatranscriptomic sequencing, (B) a hybrid capture-based approach based on a metatranscriptomic library, (C) multiplex PCR amplification for NGS platforms and (D) the Oxford Nanopore sequencing platform.
新病毒发现的基本过程和工具
在去宿主步骤 提到了要去掉rRNA 而且也提到病毒载量比较低的情况下 可以不去宿主 -
covid19.sfb.uit.no
新冠数据库 -
-
https://www.frontiersin.org/articles/10.3389/fmicb.2021.665041/full
A Novel SARS-CoV-2 Viral Sequence Bioinformatic Pipeline Has Found Genetic Evidence That the Viral 3′ Untranslated Region (UTR) Is Evolving and Generating Increased Viral Diversity这个文章对新冠的UTR部分的突变做了分析
这个文章的附件有一张表 里面有SRR号 序列数很多 对我们来说需要支持这种批量下载数据的情况
-
https://www.cdc.gov/amd/pdf/slidesets/toolkitmodule_3.5-508c.pdf
新冠病毒的仓库 gisaid 和ncbi -
https://www.mdpi.com/2075-1729/12/1/69/htm
Direct RNA Nanopore Sequencing of SARS-CoV-2 Extracted from Critical Material from Swabs
直接用Nanopore RNA测序来检测新冠-
basecalling
Nanopore Guppy base caller (v3.4.4) tool
“flow cell = FLO-MIN106” and “kit = SQK-RNA002”. -
质控
PycoQC (v2.5.0.21) software -
过滤
NanoFilt (v2.7.0)
minimum read length ≥500 nt and read quality ≥8. -
去宿主和其他微生物(去污染)
去除人GRCh38 (hg38) fungal and bacterial genome
minimap2 (v2.17–r941) -
提取新冠病毒
samtools (v1.7) view unmapped reads and reads with mapping quality lower than 10 -
对齐 call 突变
minimap2
BCFtools mpileup\call
使用Integrative Genomic Viewer (IGV) (v2.8.2) 查看突变
-
-
新冠的命名
https://covariants.org/
https://cov-lineages.org/ Pango的官网 -
-
https://terra.bio/workflow-updates-to-the-covid-19-workspace-better-viral-assembly-and-phylogenetics-with-nextstrain/
terra 对sars 流程的更新说明
公开流程仓库
https://app.terra.bio/#workspaces/pathogen-genomic-surveillance/COVID-19https://support.terra.bio/hc/en-us/articles/360041068771
2020.08.23.20178236v1.full.sars.kraken2.terra.pdf这个文章里面提到 使用kraken2检测其他病毒 主要用于排除新冠和其他病毒的交叉感染
We used Kraken2 (46) to identify other viral taxa present in NP swab samples from COVID positive patients, excluding those removed by filters i and ii described above. To do so, we ran the classify_single workflow on all reads from all samples (with kraken2_db_tgz=”gs://pathogen-public-dbs/v1/kraken2-broad-20200505.tar.zst”, krona_taxonomy_db_kraken2_tgz=”gs://pathogen-public-dbs/v1/krona.taxonomy-20200505.tab. zst”, ncbi_taxdump_tgz=”gs://pathogen-public-dbs/v1/taxdump-20200505.tar.gz”, trim_clip_db=”gs://pathogen-public-dbs/v0/contaminants.clip_db.fasta”, spikein_db=”gs://pathogen-public-dbs/v0/ERCC_96_nopolyA.fasta”). Our kraken2 database wasTerra的这个流程是针对illumina二代的情况
https://app.terra.bio/#workspaces/pathogen-genomic-surveillance/COVID-19_Broad_Viral_NGS -
-
https://artic.readthedocs.io/en/latest/primer-schemes/
artic的引物设计
该图来自该nature的文献 对多重PCR做了详细说明
https://www.nature.com/articles/nprot.2017.066
Multiplex PCR method for MinION and Illumina sequencing of Zika and other virus genomes directly from clinical samples
文档中提到了这个方法主要是用于在临床中 宏基因测序的病毒载量很低
genome sequencing directly from clinical samples (i.e., without isolation and culture) remains challenging for viruses such as Zika, for which metagenomic sequencing methods may generate insufficient numbers of viral reads.
这个文章还提到了一个在线引物设计工具（引物设计是实验的一个关键环节 这类工具我们可以做到系统里面 甚至做成一个app）
-
-
