病原宏基因组测序分析（mNGS)¶

2024.10.3

测序原理¶

BGISEQ-50

BGISEQ-50是一款具有领先水平的高通量基因测序平台，以联合探针锚定聚合技术（cPAS）及改进的DNA纳米球（DNB）为核心技术支撑。采用加载系统集成化设计，制备好的DNB无需测序仪以外的配套设备即可自动化完成样本加载到芯片上。

BGISEQ-50

参考视频：https://www.bilibili.com/video/BV1Y7411Y7qb/

产品特点 PRODUCT FEATURE

8 Gb per run
仅需18小时以内即可完成SE50
支持SE35和SE50读长
支持科研，医学临床，农业等领域的测序应用和数据分析

行业应用

病原识别
科研
农业
肿瘤筛查

性能参数 PERFORMANCE PARAMETER

Q30≥75%（高于Q30的碱基百分比是特定标准文库通过整体运行平均所得）

支持读长	SE50、SE35	测序时间	10.5-15小时
有效Reads 数	160M	平均数据产量	5.6-8Gb/Run

分析流程¶

flowchart TD
    A[Rawdata] --bowtie2去除宿主reads--> B{Cleandata}
    B -- Kraken2基于reads --> C[细菌丰度表]
    B -- Spades组装 --> D[Contigs]
    C --差异分析--> F[Final Results]
    D --Prodigal注释--> E[基因]
    E --Salmon定量--> G[基因丰度表]
    G -->F[Final Results]

workflow

分析步骤¶

00 数据质控QC¶

fastp 自动去除低质量序列
multiQC 合并QC报告

01 去除宿主序列¶

bowtie2 比对人类基因组序列，保留不能比对上的序列（认为是微生物序列）

bowtie2 -q --very-sensitive-local -x /pub/genomes/hg38/index/bowtie2/hg38.fa -p 24 -U V350042627_L03_5.fq.gz --un-gz  ./cleandata/V350042627_L03_5.unmap.fq.gz >/dev/null  2>./cleandata/V350042627_L03_5.mappingrate.txt
# -q代表格式是fastq
# --very-sensitive-local 代表严格比对
# -x代表基因组index
# -p代表线程数
# -U代表是单端测序
# --un-gz代表输出单端位比对的序列并且压缩
# 2>把日志信息输出到./cleandata/V350042627_L03_5.mappingrate.txt
# 写一个循环
# 生成样品列表
ls *.gz|cut -d "." -f1 > id.txt # 列出所有.gz结尾的文件，并且按照“.”分割，取第一列，把结果赋予id.txt这个文件
# 通过id.txt写循环
 while read id; do echo "bowtie2 -q --very-sensitive-local -x /pub/genomes/hg38/index/bowtie2/hg38.fa -p 24 -U ${id}.fq.gz --un-gz 
 ./cleandata/${id}.unmap.fq.gz  >/dev/null  2>./cleandata/${id}.mappingrate.txt"; done<id.txt
# 读取每行，用id代替每行的内容，输出文字“bowtie2 -q --very-sensitive-local -x /pub/genomes/hg38/index/bowtie2/hg38.fa -p 24 -U ${id}.fq.gz --un-gz ./cleandata/${id}.unmap.fq.gz  >/dev/null  2>./cleandata/${id}.mappingrate.txt”

02 组装¶

spades 组装微生物基因组，把reads组装成contigs用于后续注释

spades.py --s1 test.fq.gz  -m 110 -o multi_metaspades --only-assembler -t 24
# --s1 单端测序
# -m 内存最大占用（G）
# -o 输出目录
# --only-assembler 仅组装
# -t 线程数

03 注释¶

为加快速度，这里只做基因结构注释，基因功能注释用其他方法

prodigal 用于原核生物的基因结构注释

prodigal -i all_406.fasta -o all_406.genes -a all_406.proteins.faa -d all_406.cds.fna -f gff -p meta
# -i 输入核酸序列
# -o 输出文件名
# -a 输出蛋白序列
# -d 输出mRNA序列
# -f 输出格式
# -p 运行模式
echo "prodigal -i all_406.fasta -o all_406.genes -a all_406.proteins.faa -d all_406.cds.fna -f gff -p meta" > 03prodigal.sh
nohup bash 03prodigal.sh &

04 去重¶

为了加快后续计算速度，这里用cd-hit对核酸序列进行去重复

Cd-hit

cd-hit-est -i all_406.cds.fna -o all_406.cds.cdhit.fna -c 0.95 -aS 0.9 -M 10000 -T 4

05 基因功能注释¶

eggnog （http://eggnog-mapper.embl.de）
KEGG （https://www.genome.jp/kaas-bin/kaas_main）

06 基因丰度定量¶

Salmon 这个软件最初用于转录组定量，后来发现也可以用于宏基因组的定量分析

salmon index -t all_406.cds.fna -i all_406_salmon_index
# 构建基因索引
# -t 输入基因序列
# -i 输出基因索引文件夹
salmon quant  --validateMappings --meta -p 12 -i all_406_salmon_index -l A  -r ～/mNGS/cleandata/N100012851_L01_35.unmap.fq.gz -o quants/N100012851_L01_35.quant
# 定量
# --validateMappings 启用选择比对，提高灵敏度和特异性
# 启用宏基因组模式
# -p 线程数
# -i 基因索引
# -l 测序模式，A代表自动识别
# -r 单端测序的序列
# -o 输出文件夹
## 合并表格
# 生成每个样品定量结果列表
quants=` ls quants/ | awk '{print "quants/"$1 }' | tr '\n' ' ' `
# 生成样品顺序
names=` ls quants/ | sed 's/.quant$//'| tr '\n' ' ' `
## 合并生成TPM定量文件
salmon quantmerge \
--quants $quants \
--names $names \
--column tpm \
-o all_406.tpm.txt