GATK的人类宿主的微生物检测流程PathSeq

2023-09-19 19:34:20 浏览数 (1)

PathSeq 是一个 GATK 管道，用于检测取自宿主生物体（例如人类）的短读长深度测序样本中的微生物。比如人类肿瘤测序数据，就可以使用它看看是否有微生物序列！ 下图总结了它的工作原理。该管道先对reads进行质量过滤，减去来自宿主的reads，将剩余的（非宿主）reads与微生物参考基因组比对，并生成检测到的微生物的表。结果可用于确定微生物的存在和丰度以及发现新的微生物序列。

https://drive.google.com/uc?id=1DBCy-ZySO8G3MVfB839vBIaXce307N4w

PathSeq 管道图：蓝色虚线框代表输入文件。绿色框描绘了管道的三个阶段：reads质量过滤/宿主除去、微生物比对。

软件安装

PathSeq是目前集成于GATK包之中，其运行基于Unix环境

①检查Java版本

从版本 3.6 开始，GATK 需要 Java 运行时环境版本 1.8 (Java 8)。2.6 以下的早期版本需要 JRE 1.7，而早期版本则需要 1.6。所有 Linux/Unix 和 MacOS X 系统都应预安装 JRE，但版本可能有所不同。要测试您的 Java 版本，请在 shell 中运行以下命令来查看Java版本：

代码语言：javascript复制
java -version

②安装GATK

虽然官方给出了最佳安装的实践：

(How to) Install all software packages required to follow the GATK Best Practices

但是笔者亲身体验发现，通过conda已经可以安装最新版的GATK，并且顺带安装了其他必须软件（强烈推荐conda！！！使用conda安装后运行命令可以避免自己直接书写Java命令）

代码语言：javascript复制
conda install -c bioconda gatk4

③安装samtools

在conda环境中要单独安装samtools，建议仍是conda安装

代码语言：javascript复制
conda install -c bioconda samtools

运行PathSeq pipeline

官方运行命令示例：

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

conda运行示例

注意：具体写什么参数需要根据当前安装版本的输入，具体参考gatk PathSeqPipelineSpark --help

代码语言：javascript复制
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会根据实际进行调整

输入文件准备

输入样本的bam

管道接受 BAM 格式的读取

How to Generate an unmapped BAM from FASTQ or aligned BAM - Legacy GATK Forum

使用FastqToSam将FASTQ转换为uBAM并添加读取组信息

FastqToSam函数的文档：

Tool documentation

官方示例：

代码语言：javascript复制
java -Xmx8G -jar picard.jar FastqToSam  #目前该函数已经集成在conda安装的GATK中
    FASTQ=6484_snippet_1.fastq  #first read file of pair# required!
    FASTQ2=6484_snippet_2.fastq  #second read file of pair
    OUTPUT=6484_snippet_fastqtosam.bam  #required!
    READ_GROUP_NAME=H0164.2  #required; changed from default of A
    SAMPLE_NAME=NA12878  #required!
    LIBRARY_NAME=Solexa-272222  #required
    PLATFORM_UNIT=H0164ALXX140820.2 
    PLATFORM=illumina  #recommended
    SEQUENCING_CENTER=BI 
    RUN_DATE=2014-08-20T00:00:00-0400

conda运行示例（对于双端测序）：

代码语言：javascript复制
gatk FastqToSam 
    --FASTQ test_R1.fastq.gz 
    --FASTQ2 test_R2.fastq.gz 
    --OUTPUT test.bam 
    --SAMPLE_NAME test_sample #样本名需要自己定义且是必须参数

有关选择参数的一些信息：

  • 对于配对的reads，分别为第一个读取文件和第二个读取文件指定每个 FASTQ 文件 FASTQ 和 FASTQ2 。
  • 对于单端reads，请使用 FASTQ 指定输入文件。
  • 如果未指定，则自动检测 QUALITY_FORMAT 。

人类参考基因组/微生物参考基因组及相关文件

GATK认为的“正确“的参考基因组应包括：

  • 主 FASTA文件
  • 附有以 .dict 结尾的字典文件
  • 以 .fai 结尾的索引文件

常见微生物参考基因组下载链接

参考基因组后缀是：genomic.fna.gz，PS：参考序列、NCBI 目录和分类档案之间经常存在不一致，导致构建失败。为了最大限度地减少此问题，请确保在同一日期检索输入文件。

①PathSeq的best-practice提供了包含细菌，病毒，古生菌，真菌和原生动物的参考基因组文件

需要通过Google Cloud SDK下载

代码语言：javascript复制
gs://gatk-best-practices/pathseq/resources/pathseq_microbe.tar.gz
gs://gatk-best-practices/pathseq/resources/pathseq_microbe_list.txt
gs://gatk-best-practices/pathseq/resources/pathseq_taxonomy.tar.gz

②refseq/release提供的是各个kingdom下所有的参考基因组，并把他们打包分开，如真菌就分为了27份：

Index of /refseq/release

代码语言：javascript复制
#批量下载
wget https://ftp.ncbi.nlm.nih.gov/refseq/release/fungi/fungi*genomic.fna.gz
#批量解压
gunzip fungi*genomic.fna.gz
#fna文件合并
for i in $(ls *.fna);
do
  cat ${i}
  # 在每一个fasta输出之后再输出一个空行
  echo
done > ./fungi.genomic.fasta

③refseq/genome提供的是每个物种各自的参考基因组，如真菌就包含为了约540个种的参考基因组：

Index of /genomes/refseq

代码语言：javascript复制
#使用ncbi-genome-download批量下载
#ncbi-genome-download可以通过conda安装从而实现refseq的批量下载
#（本次运行最终只下载了313个，失败原因未知)
ncbi-genome-download -s refseq -F fasta -o ~/fungi_genome fungi

④人类hg38版本的参考基因组从此链接下载：

⑤PathSeq的best-practice同样提供了人类参考基因组（从运行结果判断和上述数据来源应该一致）

同样需要通过Google Cloud SDK下载

代码语言：javascript复制
gs://gatk-best-practices/pathseq/resources/pathseq_microbe.tar.gz
gs://gatk-best-practices/pathseq/resources/pathseq_microbe_list.txt
gs://gatk-best-practices/pathseq/resources/pathseq_taxonomy.tar.gz
gs://gatk-best-practices/pathseq/resources/taxdump.tar.gz

生成参考基因组的字典和索引文件

创建 FASTA 序列字典文件

使用 CreateSequenceDictionary 工具从 FASTA 文件创建 .dict 文件。仅需要指定输入，该工具将自动适当地命名输出。

代码语言：javascript复制
gatk CreateSequenceDictionary -R ref.fasta #现在此工具已经整合至conda的GATK工具内

这会生成一个名为 ref.dict 的 SAM 样式头文件，描述 FASTA 文件的内容。

代码语言：javascript复制
@HD VN:1.5
@SQ SN:20   LN:63025520 M5:0dec9660ec1efaaf33281c0d5ea2560f UR:file:/Users/vdauwera/Desktop/germline_mini/ref/ref.fasta

创建FASTA索引文件

我们使用 Samtools 中的 faidx 命令来准备 FASTA 索引文件。该文件描述了 FASTA 文件中每个重叠群的字节偏移量，使我们能够准确计算在 FASTA 文件中的特定基因组坐标处找到特定参考碱基的位置。

代码语言：javascript复制
samtools faidx ref.fasta
# 环境中应自己安装samtools，该函数未集成于GATK

这会生成一个名为 ref.fasta.fai 的文本文件，其中每个 FASTA 重叠群每行一条记录。每条记录都包含重叠群、大小、位置、basesPerLine 和 bytesPerLine。上面生成的索引文件如下所示：

代码语言：javascript复制
20  63025520    4   60  61

这表明我们的 FASTA 文件包含 20 号染色体，长度为 63025520 个碱基，然后是文件中的坐标。

宿主和微生物 BWA 索引镜像

BWA 索引镜像只能通过BwaMemIndexImageCreator创建

注意！！！该步骤极慢，因为GATK代码导致该镜像构建采取了单核心 大内存的运行方式，对于15G的fasta文件，单核运行约20h，占用内存约100G

代码语言：javascript复制
gatk BwaMemIndexImageCreator -I host.fasta
gatk BwaMemIndexImageCreator -I microbe.fasta

生成宿主的k-mer库文件

PathSeqBuildKmers 工具根据宿主的参考 FASTA 文件创建 k-mers 库。使用以下命令创建宿主参考基因组中所有 k-mers 的哈希集：

代码语言：javascript复制
gatk PathSeqBuildKmers 
    --reference host.fasta 
    -O host.hss

best-practice推荐命令

代码语言：javascript复制
gatk PathSeqBuildKmers 
    --referencePath pathseq_host.fa 
    --bloomFalsePositiveProbability 0.001 
    --kSize 31 
    --kmerMask 15 
    -O pathseq_host.bfi
#The k-mer length is 31 and there is a masked base at position 16 in each k-mer. 
#The set is represented using a Bloom filter with false positive probability 1.08e-4.

构建微生物的分类文件

下载最新的 RefSeq 入藏目录 RefSeq-releaseXX.catalog.gz，其中 XX 是最新的 RefSeq 版本号，网址为：ftp://ftp.ncbi.nlm.nih.gov/refseq/release/release-catalog/

下载 NCBI 分类数据文件转储（无需解压存档）：ftp://ftp.ncbi.nlm.nih.gov/pub/taxonomy/taxdump.tar.gz

使用PathSeqBuildReferenceTaxonomy构建分类文件:

示例命令：

代码语言：javascript复制
gatk PathSeqBuildReferenceTaxonomy 
    -R microbe.fasta 
    --refseq-catalog RefSeq-releaseXX.catalog.gz 
    --tax-dump taxdump.tar.gz 
    -O microbe.db

best-practice命令：

代码语言：javascript复制
gatk PathSeqBuildReferenceTaxonomy 
    --reference pathseq_microbe.fa 
    --refseqCatalogPath RefSeq-release81.catalog.gz 
    --taxdumpPath taxdump.tar.gz 
    --minNonVirusContigLength 2000 
    -O pathseq_taxonomy.db

#The minimum non-virus contig length parameter set to 2000 bp to remove low-quality assembly contigs that often contain artifacts.

故障排除

Java 内存不足错误

要使用 --java-options 将内存限制增加到 4GB：

代码语言：javascript复制
gatk --java-options "-Xmx4G" ...

通常应将其设置为大于所有参考文件之和的值。

临时文件/tmp空间不足错误

示例更改了该java命令/tmp文件输出路径：

代码语言：javascript复制
gatk --java-options "-Xmx16G -Djava.io.tmpdir=/store/gatk_tmp/" BwaMemIndexImageCreator -I ./microbes_db.fa --tmp-dir /store/gatk_tmp/

结果解读

PathSeq 输出文件是：

  • output.pathseq.bam：包含与微生物参考对齐的所有高质量非宿主读取。
  • output.pathseq.txt：输入样本微生物组成表，可以将其导入 Excel 查看：
  • 每行提供分类树中单个节点的信息。始终列出与树顶部相对应的“根”节点。分类信息右侧的列是：
  • Score ：根据与该分类单元中对齐的read数量，指示该分类单元存在的证据量。这通过将读数的权重除以每个可能的命中来考虑由于模糊映射读数而导致的不确定性。它还通过基因组长度进行标准化（需要通过参数-divide-by-genome-length true 将原有的score*（1,000,000/bacterial genome size)。
  • Score_normalized ：与得分相同，但在每个kingdom内标准化为总和为 100。
  • reads ：映射读取的数量（模糊 明确）
  • unambigously : 明确映射的读取数
  • Reference_length ：所发现微生物的参考基因组长度

在此示例中，可以看到 PathSeq 检测到 189,580 个读数，这些读数映射到大肠杆菌 K-12 MG1655 的菌株参考。该读取计数沿树（种、属、科等）向上传播到根节点。如果存在其他物种，它们的读取计数将被列出并添加到其相应的祖先分类类别中。

教程参考链接：https://gatk.broadinstitute.org/hc/en-us/articles/360035889911--How-to-Run-the-Pathseq-pipeline

延伸思考

  • 从wgs或者wes这样的DNA测序数据里面更容易找到微生物呢，还是从转录组数据更容易呢？
  • 单细胞转录组数据如何找微生物呢？
  • 加入肿瘤的测序数据里面分析到了微生物，可靠性该如何评判呢？
工具 镜像 内存 数据 索引

0 人点赞