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 ./
gatk BwaMemIndexImageCreator -I host.fasta
gatk BwaMemIndexImageCreator -I microbe.fasta
gatk PathSeqBuildKmers \
--reference host.fasta \
-O host.hss
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
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 # 输入样本的微生物组成表
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会根据实际进行调整
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 |
每行提供分类树中单个节点的信息。始终列出与树顶部相对应的“根”节点。分类信息右侧的列是:
原创声明:本文系作者授权腾讯云开发者社区发表,未经许可,不得转载。
如有侵权,请联系 cloudcommunity@tencent.com 删除。
原创声明:本文系作者授权腾讯云开发者社区发表,未经许可,不得转载。
如有侵权,请联系 cloudcommunity@tencent.com 删除。