GATK的人类宿主的微生物检测流程PathSeq和在空转上的运用

作者，Evil Genius

今天我们在外显子数据分析的基础上看一看最近很火的微生物检测的分析，主要还是用到GATK，其中PathSeq 是一个 GATK pipeline，用于检测取自宿主生物体（当然我们主要是做人的）的短读长深度测序样本中的微生物。比如人类肿瘤测序数据，就可以使用它看看是否有微生物序列以及微生物的种类。

分析示意图如下

先对reads进行质量过滤，减去来自宿主的reads，将剩余的（非宿主）reads与微生物参考基因组比对，并生成检测到的微生物的表。结果可用于确定微生物的存在和丰度以及发现新的微生物序列。

那么我们首先第一步，准备参考基因组，其中包括宿主的参考基因组（主要是人）和微生物的参考基因组。一般而言，数据库往往从NCBI进行下载并加以整理。

建议下载最新的微生物参考文件

gsutil cp gs://gcp-public-data--broad-references/hg38/v0/CrossSpeciesContamination/CrossSpeciesContaminant/pathseq_microbe.fa ./

或者 下载资源包中的数据，二选一

gsutil  cp gs://gatk-best-practices/pathseq/resources/pathseq_microbe.tar.gz ./

pathseq_microbe.tar.gz这个压缩包包含了以下四个文件，总文件大小 90G

如果分析得过程有新的想加入的微生物参考基因组，可以直接添加并创建字典即可。

创建示例如下

gatk BwaMemIndexImageCreator -I host.fasta
gatk BwaMemIndexImageCreator -I microbe.fasta

Generate the host k-mer library file

gatk PathSeqBuildKmers \
    --reference host.fasta \
    -O host.hss

Build the taxonomy file（构建微生物的分类文件）

Download the latest RefSeq accession catalog RefSeq-releaseXX.catalog.gz, where XX is the latest RefSeq release number, at:
ftp://ftp.ncbi.nlm.nih.gov/refseq/release/release-catalog/
Download NCBI taxonomy data files dump (no need to extract the archive):
ftp://ftp.ncbi.nlm.nih.gov/pub/taxonomy/taxdump.tar.gz Assuming these files are now in your current working directory, build the taxonomy file using PathSeqBuildReferenceTaxonomy:

gatk PathSeqBuildReferenceTaxonomy \
    -R microbe.fasta \
    --refseq-catalog RefSeq-releaseXX.catalog.gz \
    --tax-dump taxdump.tar.gz \
    -O microbe.db

整体的过程如下

#!/bin/bash
set -eu
GATK_HOME=/path/to/gatk
REFSEQ_CATALOG=/path/to/RefSeq-releaseXX.catalog.gz
TAXDUMP=/path/to/taxdump.tar.gz

echo "Building pathogen reference..."
$GATK_HOME/gatk BwaMemIndexImageCreator -I microbe.fasta
$GATK_HOME/gatk PathSeqBuildReferenceTaxonomy -R microbe.fasta --refseq-catalog $REFSEQ_CATALOG --tax-dump $TAXDUMP -O microbe.db

echo "Building host reference..."
$GATK_HOME/gatk BwaMemIndexImageCreator -I host.fasta
$GATK_HOME/gatk PathSeqBuildKmers --reference host.fasta -O host.hss

运行PathSeq pipeline，运行的命令也很简单，会分析外显子panel的对于PathSeq应该是轻车熟路,下面是官方的运行示例

gatk PathSeqPipelineSpark \
    --input test_sample.bam \ #输入样本的bam,这里面包含了人和微生物的reads信息。
    --filter-bwa-image hg19mini.fasta.img \ #人类参考基因组的BWA索引镜像
    --kmer-file hg19mini.hss \ #根据人类参考基因组构建的k-mer库
    --min-clipped-read-length 70 \ #设置排除假阳性的阈值，越高则比对到的外源序列越少
    --microbe-fasta pathseq_microbe.fa \ #待检测微生物的参考基因组
    --microbe-bwa-image pathseq_microbe.fa.img \ #待检测微生物参考基因组的BWA索引镜像
    --taxonomy-file microbe.db \ #待检测微生物的分类学文件
    --output output.pathseq.bam \ # 包含与微生物参考对齐的所有高质量非宿主读取。
    --scores-output output.pathseq.txt # 输入样本的微生物组成表

conda环境下运行示例(参考生信技能树)

gatk PathSeqPipelineSpark \
    --input test_sample.bam \ #输入样本的bam
    --filter-bwa-image hg38.fasta.img \ #人类参考基因组的BWA索引镜像
    --kmer-file hg38.hss \ #根据人类参考基因组构建的k-mer库
    --min-clipped-read-length 70 \ #设置排除假阳性的阈值，越高则比对到的外源序列越少
    --microbe-bwa-image microbe.fasta.img \ #待检测微生物参考基因组的BWA索引镜像
    --microbe-dict  microbe.fasta.dict \ #待检测微生物参考基因组的字典文件
    --taxonomy-file microbe.db \ #待检测微生物的分类学文件
    --output output.pathseq.bam \ # 包含与微生物参考比对到的所有非宿主读取。
    --scores-output output.pathseq.txt \# 输入样本的微生物组成表
    --read-filter WellformedReadFilter \ #开启过滤器保证输入SAM格式正确
    --divide-by-genome-length true \ #根据参考基因组长度对score进行标准化
    --java-options "-Xmx48G" #指定该程序的运行内存,宜大不宜小,实际运行过程中spak会根据实际进行调整

结果解读

PathSeq 输出文件是：

output.pathseq.bam：包含与微生物参考对齐的所有高质量非宿主读取。
output.pathseq.txt：输入样本微生物组成表，可以将其导入 Excel 查看。

tax_id	taxonomy	type	name	kingdom	score	score_normalized	reads	unambiguous	reference_length
1	root	root	root	root	189580	100	189580	189580	0
131567	root cellular_organisms	no_rank	cellular_organisms	root	189580	100	189580	189580	0
2	... cellular_organisms Bacteria	superkingdom	Bacteria	Bacteria	189580	100	189580	189580	0
1224	... Proteobacteria	phylum	Proteobacteria	Bacteria	189580	100	189580	189580	0
1236	... Proteobacteria Gammaproteobacteria	class	Gammaproteobacteria	Bacteria	189580	100	189580	189580	0
91347	... Gammaproteobacteria Enterobacterales	order	Enterobacterales	Bacteria	189580	100	189580	189580	0

解释说明

每行提供分类树中单个节点的信息。始终列出与树顶部相对应的“根”节点。分类信息右侧的列是：

• Score ：根据与该分类单元中对齐的read数量，指示该分类单元存在的证据量。这通过将读数的权重除以每个可能的命中来考虑由于模糊映射读数而导致的不确定性。它还通过基因组长度进行标准化（需要通过参数-divide-by-genome-length true 将原有的score*（1,000,000/bacterial genome size)。
• Score_normalized ：与得分相同，但在每个kingdom内标准化为总和为 100。
• reads ：映射读取的数量（模糊+明确）
• unambigously : 明确映射的读取数
• Reference_length ：所发现微生物的参考基因组长度.

今天我们先初步了解一下，关于病原微生物的检测过程，还是有很多需要注意的地方，尤其是病原微生物如何指导临床用药，需要更多的关注。

其中这个方法在空间转录组上已经有了很好的运用，在文章Effect of the intratumoral microbiota on spatial and cellular heterogeneity in cancer中，就运用到了该方法检测微生物的空间分布特点。

文章解读在单细胞空间宏基因组揭秘微生物群对癌症空间和细胞异质性的作用2。

生活很好，有你更好