不知覺間，距離我寫下第一篇關(guān)于 WGS 數(shù)據(jù)分析系列的文章已經(jīng)過去了三年多（WGS 系列文章），時(shí)間真的快啊。
很多朋友從當(dāng)時(shí)的文章中得到了啟發(fā)，對(duì)此我也很開心。在后來的日子里，我又合作完成了多個(gè)大規(guī)模的人類基因組學(xué)科研項(xiàng)目，在這個(gè)過程中關(guān)于大規(guī)模的 WGS 數(shù)據(jù)分析（數(shù)量從數(shù)千到十萬、乃至百萬級(jí)別）已經(jīng)是家常便飯。
我自己也積累了更多屬于大數(shù)據(jù)基因組學(xué)研究的分析策略、方法論和最佳實(shí)踐，還在整理之中，現(xiàn)在把其中一部分做成 Python 工具包分享給大家。這個(gè)工具一年前寫了第一個(gè)版本（當(dāng)時(shí)只在我的知識(shí)星球里分享過一次），過程中也調(diào)試了碰到的問題，日趨完善，后續(xù)還將迭代更新。
這個(gè)工具我將其命名為 ilus (/i:l?s/)，這是我看過的一部美劇——《無垠的太空》中通過星環(huán)抵達(dá)的第一個(gè)系外類地行星的名字。它是一個(gè)全面的、輕量的、可拓展且易用的半自動(dòng)化全基因組（Whole genome sequencing, WGS）和全外顯子（Whole exom sequencing，WES）分析流程生成器，是以前我 這篇文章 所提供代碼的高級(jí)版本。
簡介
ilus 的用途就是生成完整的 WGS/WES 分析流程，但ilus不執(zhí)行流程的具體步驟。你需要自己手動(dòng)投遞任務(wù)，不過執(zhí)行過程不再依賴ilus，這也是為何我稱之為半自動(dòng)化的原因（這其實(shí)是一個(gè)優(yōu)點(diǎn)，下文我會(huì)說到）。雖然如此，但 ilus 會(huì)幫你將最重要的流程和分析步驟生成出來，你只要按步驟投遞就可以了。
目前 ilus 含有三個(gè)功能模塊，分別是：
第一、
WGS全基因組數(shù)據(jù)分析流程模塊。這個(gè)模塊基于 GATK 最佳實(shí)踐，采用bwa-mem + GATK，這和我之前的 WGS 系列教程是相互呼應(yīng)的，但在 ilus 中被我做的更好用了。這個(gè)模塊集成了比對(duì)、排序、同一個(gè)樣本多l(xiāng)ane 數(shù)據(jù)合并、標(biāo)記重復(fù)序列（Markduplicates）、HaplotypeCaller gvcf 生成、多樣本聯(lián)合變異檢測（Joint-calling）和變異質(zhì)控（Variant quality score recalibrator, VQSR） 這 5 個(gè)完整的過程。而且這個(gè)模塊同樣適用于WES數(shù)據(jù)的分析，只需要將配置文件variant_calling_interval設(shè)置為WES 的外顯子捕獲區(qū)間即可（下文詳述）。第二、
genotype-joint-calling聯(lián)合變異檢測模塊。這個(gè)模塊是從 ilus WGS 中分離出來的，目的是讓我們可以從樣本的 gvcf 文件 （注意 gvcf 文件和 VCF 文件是完全不同的）開始進(jìn)行變異檢測，這個(gè)功能很有用。特別是在碰到需要分多批次完成 WGS/WES 數(shù)據(jù)分析的時(shí)候（在大型基因組學(xué)項(xiàng)目中，這是很常見的），我們可以分批次生成 gvcf，最后再整理一個(gè)總的 gvcf 文件列表，然后用該模塊完成后續(xù)步驟，這增加了分析流程的靈活度。第三、
VQSR變異質(zhì)控模塊，同樣是從 ilus WGS 中分離出來的，目的是方便我們可以對(duì)變異結(jié)果進(jìn)行獨(dú)立的VQSR質(zhì)控。
需要注意的是 ilus 不包含對(duì)原始 fastq 數(shù)據(jù)的質(zhì)控，ilus
流程默認(rèn)你輸入的測序數(shù)據(jù)都是清洗好的 clean data。
ilus 不直接運(yùn)行任務(wù)有如下兩個(gè)方面的考慮：
- 第一，避免在 ilus 的核心功能之外做關(guān)于任務(wù)調(diào)度方面的優(yōu)化。不同的計(jì)算集群（本地和云上），作業(yè)被調(diào)度的方式是多種多樣的，如果將這些情況都一一考慮進(jìn)去，ilus 會(huì)變得臃腫復(fù)雜，并且還不一定能夠符合真實(shí)的需要，反而會(huì)導(dǎo)致一部分人無法有效使用 ilus，也容易在跟多系統(tǒng)任務(wù)管理纏斗的過程中丟失 ilus 真正要解決的問題。因此，作為一個(gè)輕量級(jí)的工具，我在設(shè)計(jì) ilus 的時(shí)候從一開始就沒將自動(dòng)投遞和運(yùn)行任務(wù)的功能考慮在內(nèi)。我更希望它作為一個(gè)框架程序，嚴(yán)格依據(jù)你的輸入數(shù)據(jù)和配置文件的信息，生成符合你分析需求的流程腳本。
如果要實(shí)現(xiàn)作業(yè)的自動(dòng)投遞和監(jiān)控，ilus 希望將來可以通過外部插件的形式來實(shí)現(xiàn)，但目前需要手動(dòng)去投遞這些按次序生成的腳本。
- 第二，增加流程靈活性和可操控性。ilus 所生成的流程都有一個(gè)特點(diǎn)，就是它們都將完全獨(dú)立于 ilus，不會(huì)再依賴 ilus 的任何功能，這個(gè)時(shí)候即使你將 ilus 徹底從系統(tǒng)中卸載刪除掉都沒關(guān)系（這個(gè)特性是我有意而為的）。另外，ilus 所生成的可執(zhí)行腳本(shell) 每一行都可以獨(dú)立運(yùn)行，彼此之間互不影響。
這樣做的好處是，你可以按照計(jì)算機(jī)集群的特點(diǎn)將腳本拆分成若干個(gè)子腳本（如果你的樣本不多或者集群資源不充足，也可以不拆分），然后再分別投遞任務(wù)，這樣可以大大加快任務(wù)的完成速度，這也是目前并行投遞 ilus 任務(wù)的唯一方式，我也非常建議你應(yīng)該這樣去做（我的項(xiàng)目基本都是通過這樣的方式來完成的）。
至于每一步要拆成多少個(gè)子腳本取決于你自己的具體情況。舉個(gè)例子，比如你一共有 100 個(gè)樣本，第一步的 xxx.step1.bwa.sh 比對(duì)腳本中將一共有 100 個(gè)比對(duì)命令，每一行都是一個(gè)樣本的比對(duì)指令。由于這 100 個(gè)命令彼此獨(dú)立互不依賴，因此你可以放心地將該步驟拆分為 100（或者任意小于 100 ）個(gè)子腳本，然后再手動(dòng)投遞這些任務(wù)。
至于如何將一個(gè)完整的執(zhí)行腳本拆分為多個(gè)，你既可以自己寫程序完成，也可以使用我在 ilus 中提供的 yhbatch_slurm_jobs.py 程序來完成，但要注意，我提供的這個(gè)程序是基于 slurm 任務(wù)調(diào)度系統(tǒng)的（尚未測試過其它系統(tǒng)），不一定符合你的集群配置，但它的作用和意義我剛剛也做了說明，如果你覺得這個(gè)程序不能直接滿足你的需求，你可以參照和修改。雖然 ilus 不做任務(wù)調(diào)度和執(zhí)行，但是這個(gè)方法卻能夠增加流程控制的靈活性和操控性。
此外，當(dāng)你有成千上萬的樣本需要分析時(shí)，實(shí)現(xiàn)對(duì)任務(wù)完成狀態(tài)的監(jiān)控是極其重要的一項(xiàng)任務(wù)。這個(gè)過程如果你打算手動(dòng)來檢查，那么一定是低效、易出錯(cuò)且極易讓人發(fā)狂的。
我在 ilus 中充分考慮到了這一點(diǎn)，因此在生成流程的時(shí)候會(huì)為每個(gè)任務(wù)添加一個(gè)可識(shí)別的結(jié)束標(biāo)記，我們只需要查看這個(gè)標(biāo)記就行了（參考下文 WGS 的例子）。
有結(jié)束標(biāo)志其實(shí)還不夠，一旦我們的任務(wù)數(shù)量龐大，假如每次都要手動(dòng)打開某些文件檢查這個(gè)標(biāo)志的話，那未免太過于麻煩。因此，我在 ilus 中還同時(shí)實(shí)現(xiàn)了一個(gè)專門用來檢查任務(wù)作業(yè)完成狀態(tài)的程序，具體用法同樣參考下文 WGS 例子。
如何安裝
ilus 是基于 Python 編寫的，同時(shí)支持 Python3.7+ 和 Python2.7+，穩(wěn)定版本的代碼發(fā)布至 PyPI。因此要使用 ilus, 直接通過 pip 這個(gè) Python 包管理工具進(jìn)行安裝即可，非常方便：
$ pip install ilus
該命令除了主程序 ilus 之外，還會(huì)自動(dòng)將 ilus 所依賴的其它 Python 包自動(dòng)裝上。安裝完成之后，在命令行中執(zhí)行 ilus，如果能看到類似如下的內(nèi)容，那么就說明安裝成功了。
$ ilus
usage: ilus [-h] {WGS,genotype-joint-calling,VQSR} ...
ilus: error: too few arguments
ilus 的源代碼托管在 github 中：https://github.com/ShujiaHuang/ilus
快速開始
執(zhí)行 ilus --help 可以看到 WGS, genotype-joint-calling 和 VQSR 這三個(gè)功能模塊。
$ ilus --help
usage: ilus [-h] {WGS,genotype-joint-calling,VQSR} ...
ilus: A WGS analysis pipeline.
optional arguments:
-h, --help show this help message and exit
ilus commands:
{WGS,genotype-joint-calling,VQSR}
WGS Creating pipeline for WGS(from fastq to genotype VCF)
genotype-joint-calling Genotype from GVCFs.
VQSR VQSR
下面，我通過例子逐一介紹如何使用這三個(gè)模塊。
全基因組和全外顯子數(shù)據(jù)分析
全基因組數(shù)據(jù)分析流程（WGS）的運(yùn)行腳本通過 ilus WGS 來生成，用法如下：
$ ilus WGS --help
usage: ilus WGS [-h] -C SYSCONF -L FASTQLIST [-P WGS_PROCESSES]
[-n PROJECT_NAME] [-f] [-c] -O OUTDIR
optional arguments:
-h, --help show this help message and exit
-C SYSCONF, --conf SYSCONF
YAML configuration file specifying details about
system.
-L FASTQLIST, --fastqlist FASTQLIST
Alignment FASTQ Index File.
-O OUTDIR, --outdir OUTDIR
A directory for output results.
-n PROJECT_NAME, --name PROJECT_NAME
Name of the project. Default value: test
-P WGS_PROCESSES, --Process WGS_PROCESSES
Specific one or more processes (separated by comma) of
WGS pipeline. Defualt value:
align,markdup,BQSR,gvcf,genotype,VQSR. Possible
values: {align,markdup,BQSR,gvcf,genotype,VQSR}
-f, --force_overwrite
Force overwrite existing shell scripts and folders.
-c, --cram Covert BAM to CRAM after BQSR and save alignment file storage.
其中，-C, -L 和 -O 是三個(gè)最重要且 必須的參數(shù)，其余的按照你的實(shí)際需要做選擇。 -O 參數(shù)比較簡單，就是一個(gè)輸出目錄，該目錄如果不存在，無需擔(dān)心，ilus 會(huì)自動(dòng)創(chuàng)建。
最重要的是 -C 和 -L 這兩個(gè)參數(shù)。 -C 參數(shù)是 ilus 的配置文件，如果沒有這個(gè)文件 ilus 就無法正確生成分析流程了，因此它十分重要；-L 是輸入文件，這個(gè)文件的格式有固定要求，一共 5 列，每一列都是流程所必須的信息。
下面，我分別對(duì)這兩個(gè)文件的格式展開說明。
首先是 -C 配置文件，你需要在文件中填寫好分析流程所需的所有程序路徑、GATK bundle 文件路徑、參考基因組 fasta 文件的路徑以及各個(gè)關(guān)鍵步驟所對(duì)應(yīng)的參數(shù)。
需要注意的是 bwa mem 的比對(duì)索引文件前綴要與配置文件的 resources.reference 的名字相同，并放在同一個(gè)文件夾里。如下：
/path/human_reference/GRCh38/
|-- human_GRCh38.fa
|-- human_GRCh38.dict
|-- human_GRCh38.fa.amb
|-- human_GRCh38.fa.ann
|-- human_GRCh38.fa.bwt
|-- human_GRCh38.fa.fai
|-- human_GRCh38.fa.pac
`-- human_GRCh38.fa.sa
human_GRCh38.fa 是 resources.reference，后面 6 份以 human_GRCh38.fa 為前綴的是 bwa mem 所需的比對(duì)索引文件，human_GRCh38.dict 也是一份必須的文件。配置文件要使用 Yaml 語法 來編寫，這里我提供一份 配置文件的模板，大家可以直接參考：
aligner:
bwa: /path/to/BioSoftware/local/bin/bwa
bwamem_options: [-Y -M -t 8]
samtools:
samtools: /path/to/BioSoftware/local/bin/samtools
sort_options: ["-@ 8"]
merge_options: ["-@ 8 -f"]
stats_options: ["-@ 8"]
bcftools:
bcftools: /path/to/BioSoftware/local/bin/bcftools
concat_options: ["-a --rm-dups all"]
bedtools:
bedtools: /path/to/BioSoftware/local/bin/bedtools
genomecov_options: ["-bga -split"]
sambamba:
sambamba: /path/to/BioSoftware/local/bin/sambamba
sort_options: ["-t 8"]
merge_options: ["-t 8"]
markdup_options: []
verifyBamID2:
# Link of VerifyBamID2: https://github.com/Griffan/VerifyBamID
verifyBamID2: /path/to/BioSoftware/local/bin/verifyBamID2
options: [
# Download from: https://github.com/Griffan/VerifyBamID/tree/master/resource
"--SVDPrefix /path/to/BioSoftware/verifyBamID2/1.0.6/resource/1000g.phase3.10k.b38.vcf.gz.dat"
]
bgzip: /path/to/BioSoftware/local/bin/bgzip
tabix: /path/to/BioSoftware/local/bin/tabix
gatk:
gatk: /path/to/BioSoftware/gatk/4.1.4.1/gatk
markdup_java_options: ["-Xmx10G", "-Djava.io.tmpdir=/your_path/cache"]
bqsr_java_options: ["-Xmx8G", "-Djava.io.tmpdir=/your_path/cache"]
hc_gvcf_java_options: ["-Xmx4G"]
genotype_java_options: ["-Xmx8G"]
vqsr_java_options: ["-Xmx10G"]
CollectAlignmentSummaryMetrics_jave_options: ["-Xmx10G"]
# Adapter sequencing of BGISEQ-500. If you use illumina(or other) sequencing system you should
# change the value of this parameter.
CollectAlignmentSummaryMetrics_options: [
"--ADAPTER_SEQUENCE AAGTCGGAGGCCAAGCGGTCTTAGGAAGACAA",
"--ADAPTER_SEQUENCE AAGTCGGATCGTAGCCATGTCGTTCTGTGAGCCAAGGAGTTG"
]
genomicsDBImport_options: ["--reader-threads 12"]
use_genomicsDBImport: false # Do not use genomicsDBImport to combine GVCFs by default
vqsr_options: [
"-an DP -an QD -an FS -an SOR -an ReadPosRankSum -an MQRankSum -an InbreedingCoeff",
"-tranche 100.0 -tranche 99.9 -tranche 99.5 -tranche 99.0 -tranche 95.0 -tranche 90.0",
"--max-gaussians 6"
]
# ``interval`` for create gvcf. The value could be a interval region file in bed format or could be a list here
interval: ["chr1", "chr2", "chr3", "chr4", "chr5", "chr6", "chr7", "chr8", "chr9",
"chr10", "chr11", "chr12", "chr13", "chr14", "chr15", "chr16", "chr17",
"chr18", "chr19", "chr20", "chr21", "chr22", "chrX", "chrY", "chrM"]
# Specific variant calling intervals.
# The value could be a file in bed format (I show you a example bellow) or a interval of list.
# Bed format of interval file only contain three columns: ``Sequencing ID``, ``region start`` and ``region end``,e.g.:
# chr1 10001 207666
# chr1 257667 297968
# These invertals could be any regions alone the genome as you wish or just set the same as ``interval`` parameter above.
variant_calling_interval: ["chr1", "chr2", "chr3", "chr4", "chr5", "chr6", "chr7", "chr8", "chr9",
"chr10", "chr11", "chr12", "chr13", "chr14", "chr15", "chr16", "chr17",
"chr18", "chr19", "chr20", "chr21", "chr22", "chrX", "chrY", "chrM"]
# variant_calling_interval: ["./wgs_calling_regions.GRCh38.interval.bed"]
# GATK bundle
bundle:
hapmap: /path/to/BioDatahub/gatk/bundle/hg38/hapmap_3.3.hg38.vcf.gz
omni: /path/to/BioDatahub/gatk/bundle/hg38/1000G_omni2.5.hg38.vcf.gz
1000G: /path/to/BioDatahub/gatk/bundle/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz
mills: /path/to/BioDatahub/gatk/bundle/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz
1000G_known_indel: /path/to/BioDatahub/gatk/bundle/hg38/Homo_sapiens_assembly38.known_indels.vcf.gz
dbsnp: /path/to/BioDatahub/gatk/bundle/hg38/Homo_sapiens_assembly38.dbsnp138.vcf.gz
# Set Reference
resources:
reference: /path/to/BioDatahub/human_reference/GRCh38/human_GRCh38.fa
如果你不知道各個(gè)步驟最佳的參數(shù)是什么的話，不妨使用我配置模板中所提供的。
在配置文件中，bwa、samtools、bcftools、bedtools、gatk、bgzip 和 tabix 都是必須的生信軟件，需要自行安裝，然后再將路徑填入到對(duì)應(yīng)的參數(shù)里。
verifyBamID2 僅用于計(jì)算樣本是否存在污染，它并不是一個(gè)必填的參數(shù)，如果你的配置文件中沒有這個(gè)參數(shù)，則代表流程不對(duì)樣本的污染情況進(jìn)行推算，如果有那么你要自行安裝并下載與之配套的 resource 數(shù)據(jù)，模板里我也告訴你該去哪里下載相關(guān)的數(shù)據(jù)了（沒注意的話，建議回頭去看看上面的配置模板信息）。
還有，配置文件 GATK 參數(shù)中的 Djava.io.tmpdir 緩存文件路徑記得重新設(shè)置為你自己的文件夾路徑。
此外，要注意配置文件中的 variant_calling_interval 參數(shù)。這是一個(gè)專門用來指定變異檢測區(qū)間的參數(shù)。在以上配置例子里，我給出了從 chr1 到 chrM 這 25 條染色體，意思就是告訴流程只對(duì)這 25 條染色體做變異檢測。如果你在參數(shù)里只列出一條染色體，或者僅僅給出一個(gè)染色體區(qū)間（區(qū)間信息按照 chrom:start-end 格式），比如 chr1:1-10000，那么 ilus 將只在你給定的這個(gè)區(qū)間里做變異檢測。
這是一個(gè)非常靈活有用的參數(shù)，因?yàn)? variant_calling_interval 可以任意指定的，除了可以按照我例子給出的賦值方式之外，你還可以將區(qū)間 文件的路徑 賦給這個(gè)參數(shù)。我們知道 WGS 和 WES 有很多步驟是完全相同的，只在變異檢測的區(qū)間上存在差別——WES數(shù)據(jù)沒有必要也不能 在全染色體上做變異檢測，它只在外顯子捕獲區(qū)域里進(jìn)行。
這個(gè)時(shí)候你只需要將外顯子捕獲區(qū)域的文件——注意是文件，這個(gè)文件的內(nèi)容可以是 .interval_list 格式、.list 格式、.intervals 格式或者 .bed 格式。其中，.list 格式和 .intervals文件格式必須如下所示：
chr1:63697-63697
chr1:101158-101158
chr1:103241-103241
chr1:104108-104108
chr1:185336-185336
chr1:261495-261495
chr1:598862-598862
chr1:601606-601606
chr1:700596-700596
chr1:725086-725086
而 .interval_list 格式和 .bed 格式參照 GATK的說明，你不需要手動(dòng)拆分成一個(gè)個(gè)的區(qū)間，只需將這個(gè)文件的路徑賦給這個(gè)參數(shù)就可以了，這時(shí)流程就成了 WES 分析流程。這也是為何 ilus 同時(shí)是一個(gè)WGS 和 WES 分析流程生成器的原因。
另外，ilus 必需的公用數(shù)據(jù)集是：gatk bundle 和基因組參考序列（fasta 文件）。
【注意】如果你項(xiàng)目的樣本量少于 10 那么 GATK 將不計(jì)算
InbreedingCoeff的值，此時(shí)配置文件中vqsr_options不需要設(shè)置-an InbreedingCoeff，可以將其去掉。
接下來是由 -L 參數(shù)指定的輸入文件，文件里包含了 WGS/WES 分析流程所必需的一切信息。這是一個(gè)需要你自己來準(zhǔn)備的文件，文件各列的格式信息如下：
- [1] SAMPLE，樣本名；
- [2] RGID，Read Group，使用
bwa mem時(shí)通過 -R 參數(shù)指定的read group； - [3] FASTQ1，F(xiàn)astq1 文件的路徑；
- [4] FASTQ2，F(xiàn)astq2 文件路徑，如果是Single End測序，沒有fastq2的話，該列用空格代替；
- [5] LANE，fastq 的 lane 編號(hào)。
這五個(gè)信息中
RGID要求嚴(yán)格，最容易出錯(cuò)。RGID一定要設(shè)置正確（正確的編寫方式參考下面的例子或者這篇文章），否則你流程會(huì)出錯(cuò)。
另外，假如某些樣本有多個(gè) lane 的測序數(shù)據(jù)，或者同一個(gè) lane 的數(shù)據(jù)被拆分成了很多個(gè)子文件，這個(gè)時(shí)候也不需要手動(dòng)合并這些 fastq 數(shù)據(jù)，只需要依照 -L 的格式要求編寫在輸入文件里即可。對(duì)于這些同屬一個(gè)樣本的數(shù)據(jù)，ilus 都會(huì)根據(jù)樣本 ID 在生成的流程中自動(dòng)幫你合并。
下面我給出一個(gè) -L 輸入文件的例子，其中樣本 HG002, HG003和 HG004 的數(shù)據(jù)就有分拆的情況（哪怕碎成一萬份也沒問題）：
#SAMPLE RGID FASTQ1 FASTQ2 LANE
HG002 "@RG\tID:CL100076190_L01\tPL:COMPLETE\tPU:CL100076190_L01_HG002\tLB:CL100076190_L01\tSM:HG002" /path/to/NA24385_CL100076190_L01_read_1.clean.fq.gz /path/to/NA24385_CL100076190_L01_read_2.clean.fq.gz CL100076190_L01
HG002 "@RG\tID:CL100076190_L02\tPL:COMPLETE\tPU:CL100076190_L02_HG002\tLB:CL100076190_L02\tSM:HG002" /path/to/NA24385_CL100076190_L02_read_1.clean.fq.gz /path/to/NA24385_CL100076190_L02_read_2.clean.fq.gz CL100076190_L02
HG003 "@RG\tID:CL100076246_L01\tPL:COMPLETE\tPU:CL100076246_L01_HG003\tLB:CL100076246_L01\tSM:HG003" /path/to/NA24149_CL100076246_L01_read_1.clean.fq.gz /path/to/NA24149_CL100076246_L01_read_2.clean.fq.gz CL100076246_L01
HG003 "@RG\tID:CL100076246_L02\tPL:COMPLETE\tPU:CL100076246_L02_HG003\tLB:CL100076246_L02\tSM:HG003" /path/to/NA24149_CL100076246_L02_read_1.clean.fq.gz /path/to/NA24149_CL100076246_L02_read_2.clean.fq.gz CL100076246_L02
HG004 "@RG\tID:CL100076266_L01\tPL:COMPLETE\tPU:CL100076266_L01_HG004\tLB:CL100076266_L01\tSM:HG004" /path/to/NA24143_CL100076266_L01_read_1.clean.fq.gz /path/to/NA24143_CL100076266_L01_read_2.clean.fq.gz CL100076266_L01
HG004 "@RG\tID:CL100076266_L02\tPL:COMPLETE\tPU:CL100076266_L02_HG004\tLB:CL100076266_L02\tSM:HG004" /path/to/NA24143_CL100076266_L02_read_1.clean.fq.gz /path/to/NA24143_CL100076266_L02_read_2.clean.fq.gz CL100076266_L02
HG005 "@RG\tID:CL100076244_L01\tPL:COMPLETE\tPU:CL100076244_L01_HG005\tLB:CL100076244_L01\tSM:HG005" /path/to/NA24631_CL100076244_L01_read_1.clean.fq.gz /path/to/NA24631_CL100076244_L01_read_2.clean.fq.gz CL100076244_L01
接下來，舉例說明 ilus WGS 的使用和整個(gè)流程的結(jié)構(gòu)。
例子1：從頭開始生成 WGS 分析流程
$ ilus WGS -c -n my_wgs_project -C ilus_sys.yaml -L input.list -O output/
這個(gè)命令的意思是，項(xiàng)目 (-n) my_wgs_project 依據(jù)配置文件 (-C) ilus_sys.yaml
和輸入數(shù)據(jù) (-L )input.list 在輸出目錄 output 中生成一個(gè) WGS
分析流程。同時(shí)流程在完成分析之后將 BAM 比對(duì)文件自動(dòng)轉(zhuǎn)為 (-c) CRAM。CRAM 和 BAM 格式和可被操作的方式幾乎相同，但可以節(jié)省大約三分之一的空間，如果不設(shè)置 -c 參數(shù)，則保留原來的 BAM 文件。
以上命令順利運(yùn)行之后，輸出目錄 output 里會(huì)生成如下 4 個(gè)文件夾：
00.shell/
01.alignment/
02.gvcf/
03.genotype/
它們分別用于存放流程產(chǎn)生的各類不同數(shù)據(jù)，其中：
-
00.shell流程shell腳本的匯集目錄； -
01.alignment以樣本為單位存放比對(duì)結(jié)果； -
02.gvcf存放各個(gè)樣本的gvcf結(jié)果； -
03.genotype存放最后變異檢測的結(jié)果。
00.shell 目錄里有分析流程執(zhí)行腳本和日志目錄，我們要投遞和執(zhí)行腳本都在這：
/00.shell
├── loginfo
│ ├── 01.alignment
│ ├── 01.alignment.e.log.list
│ ├── 01.alignment.o.log.list
│ ├── 02.markdup
│ ├── 02.markdup.e.log.list
│ ├── 02.markdup.o.log.list
│ ├── 03.BQSR
│ ├── 03.BQSR.e.log.list
│ ├── 03.BQSR.o.log.list
│ ├── 04.gvcf
│ ├── 04.gvcf.e.log.list
│ ├── 04.gvcf.o.log.list
│ ├── 05.genotype
│ ├── 05.genotype.e.log.list
│ ├── 05.genotype.o.log.list
│ ├── 06.VQSR
│ ├── 06.VQSR.e.log.list
│ └── 06.VQSR.o.log.list
├── my_wgs_project.step1.bwa.sh
├── my_wgs_project.step2.markdup.sh
├── my_wgs_project.step3.bqsr.sh
├── my_wgs_project.step4.gvcf.sh
├── my_wgs_project.step5.genotype.sh
└── my_wgs_project.step6.VQSR.sh
投遞任務(wù)運(yùn)行流程時(shí)，按順序從 step1 依次執(zhí)行到 step6 即可。loginfo/ 文件夾下記錄了各個(gè)樣本所有步驟的運(yùn)行狀態(tài)，你可以通過檢查各個(gè)任務(wù)的 .o.log.list 日志文件，獲得每個(gè)樣本是否都有成功結(jié)束的標(biāo)記。
如果成功了，你將在日志文件的末尾看到一句類似于 [xxxx] xxxx done 的標(biāo)記。用我在 ilus 中提供的程序 check_jobs_status.py，你就可以很方便地知道哪些已經(jīng)順利完成，哪些還沒有。這個(gè)腳本會(huì)幫你將那些未完成的任務(wù)全部輸出，方便檢查問題或重新執(zhí)行這部分未完成的任務(wù)。
check_jobs_status 的用法如下：
$ python check_jobs_status.py loginfo/01.alignment.o.log.list > bwa.unfinish.list
如果這個(gè) bwa.unfinish.list 文件是空的，并輸出了** All Jobs done **，那么就代表你 step1 的所有任務(wù)都成功結(jié)束了，以此類推。如果有未完成的，那么你要提取出來，重新投遞。
例子2：只生成 WGS 流程中的某個(gè)/某些步驟
有時(shí)，我們并不打算（或者沒有必要）從頭到尾完整地將 WGS 流程執(zhí)行下去，比如我們只想完成從 fastq 比對(duì)到生成 gvcf 結(jié)束，暫時(shí)不想執(zhí)行 genotype 和 VQSR——這情況在大型基因組項(xiàng)目中很普遍，這時(shí)該怎么辦呢？ilus 的 -P 參數(shù)就可以實(shí)現(xiàn)這個(gè)目的：
$ ilus WGS -c -n my_wgs_project -C ilus_sys.yaml -L input.list -P align,markdup,BQSR,gvcf -O output/
這樣就只生成從 bwa 到 gvcf 的執(zhí)行腳本，這對(duì)于需要分批次完成分析的項(xiàng)目來說是很有用的。而且 ilus 所輸出的結(jié)果是以樣本為文件夾作區(qū)分的，因此在相同的輸出目錄下，只要樣本編號(hào)是不同的，那么不同批次的數(shù)據(jù)就不會(huì)存在相互覆蓋的問題。
除此之外，-P 參數(shù)還有一個(gè)用途，那就是萬一某個(gè) WGS 步驟跑錯(cuò)了，需要重跑，那你就可以用 -P重跑特定的步驟。比如我需要重新生成 BQSR 這個(gè)步驟的運(yùn)行腳本，那么就可以這樣做：
$ ilus WGS -c -n my_wgs_project -C ilus_sys.yaml -L input.list -P BQSR -O ./output
不過，要特別注意，ilus 為了節(jié)省項(xiàng)目對(duì)存儲(chǔ)空間的消耗，只會(huì)為每一個(gè)樣本保留 BQSR 之后的總 BAM/CRAM 文件。因此，如果你要重新跑 BQSR 那就需要先確保 BQSR 的前一步（即，markdup）的 BAM 文件還沒有被被刪除。如果你在項(xiàng)目中一直使用的都是 ilus 那么是不用擔(dān)心這個(gè)問題的，因?yàn)?ilus 執(zhí)行任務(wù)時(shí)具有 原子屬性，也就是說只有當(dāng)步驟中所有過程都成功結(jié)束了才會(huì)將那些完全不需要的文件刪除掉。所以，如果 ilus 的 BQSR 沒有正常結(jié)束，那么前一步 markdup 的 BAM 文件是會(huì)被完整保留住的。
目前
-P參數(shù)能指定的分析模塊必須屬于「align,markdup,BQSR,gvcf,genotype,VQSR」中的一個(gè)或多個(gè)，并用英文逗號(hào)隔開，中間不可以有空格。
genotype-joint-calling
如果我們已經(jīng)有了各個(gè)樣本的 gvcf 數(shù)據(jù)，現(xiàn)在要用這些 gvcf 完成多樣本的聯(lián)合變異檢測（Joint-calling），那么就要用 genotype-joint-calling 了。具體用法如下：
$ ilus genotype-joint-calling --help
usage: ilus genotype-joint-calling [-h] -C SYSCONF -L GVCFLIST
[-n PROJECT_NAME] [--as_pipe_shell_order]
[-f] -O OUTDIR
optional arguments:
-h, --help show this help message and exit
-C SYSCONF, --conf SYSCONF
YAML configuration file specifying details about
system.
-L GVCFLIST, --gvcflist GVCFLIST
GVCFs file list. One gvcf_file per-row and the format
should looks like: [interval gvcf_file_path]. Column
[1] is a symbol which could represent the genome
region of the gvcf_file and column [2] should be the
path.
-O OUTDIR, --outdir OUTDIR
A directory for output results.
-n PROJECT_NAME, --name PROJECT_NAME
Name of the project. [test]
--as_pipe_shell_order
Keep the shell name as the order of `WGS`.
-f, --force Force overwrite existing shell scripts and folders.
-L 是 genotype-joint-calling 的輸入?yún)?shù)，它要一個(gè) gvcf list 文件，這個(gè)文件很簡單，由兩列構(gòu)成：第一列是每個(gè) gvcf 文件所對(duì)應(yīng)的區(qū)間或者染色體編號(hào)；第二列是這份 gvcf 文件的路徑。目前 ilus 要求各個(gè)樣本的 gvcf 都按照主要染色體（1-22、X、Y、M）分開，舉個(gè)例子：
$ ilus genotype-joint-calling -n my_project -C ilus_sys.yaml -L gvcf.list -O genotype --as_pipe_shell_order
其中 gvcf.list 的格式如下：
chr1 /path/sample1.chr1.g.vcf.gz
chr1 /paht/sample2.chr1.g.vcf.gz
chr2 /path/sample1.chr2.g.vcf.gz
chr2 /path/sample2.chr2.g.vcf.gz
...
chrM /path/sample1.chrM.g.vcf.gz
chrM /path/sample2.chrM.g.vcf.gz
這個(gè)例子里 gvcf.list 只有兩個(gè)樣本（文件順序隨意）。參數(shù) --as_pipe_shell_order 可加也可不加（默認(rèn)是不加，但我建議加），它唯一的作用就是按照 ilus WGS 流程的方式輸出執(zhí)行腳本的名字，維持和 WGS 流程一樣的次序和相同的輸出結(jié)構(gòu)。
VQSR
這個(gè)模塊只用來生成基于 GATK VQSR 的變異質(zhì)控執(zhí)行腳本。VQSR 的大致原理我以前在我的知識(shí)星球 (達(dá)爾文生信星球) 中簡要回答過。
我們?nèi)绻呀?jīng)有了最終的變異檢測（VCF格式）結(jié)果，現(xiàn)在只想用 GATK VQSR 對(duì)這個(gè)變異數(shù)據(jù)做質(zhì)控，那么就可以使用這個(gè)模塊了，用法和 genotype-joint-calling 大同小異，如下：
$ ilus VQSR --help
usage: ilus VQSR [-h] -C SYSCONF -L VCFLIST [-n PROJECT_NAME]
[--as_pipe_shell_order] [-f] -O OUTDIR
optional arguments:
-h, --help show this help message and exit
-C SYSCONF, --conf SYSCONF
YAML configuration file specifying details about
system.
-L VCFLIST, --vcflist VCFLIST
VCFs file list. One vcf_file per-row and the format
should looks like: [interval vcf_file_path]. Column
[1] is a symbol which could represent the genome
region of the vcf_file and column [2] should be the
path.
-O OUTDIR, --outdir OUTDIR
A directory for output results.
-n PROJECT_NAME, --name PROJECT_NAME
Name of the project. [test]
--as_pipe_shell_order
Keep the shell name as the order of `WGS`.
-f, --force Force overwrite existing shell scripts and folders.
跟 genotype-joint-calling 相比不同的是，ilus VQSR 的輸入文件是 VCF 文件列表，并且 每行只有一個(gè) VCF 文件路徑，例子如下：
/path/chr1.vcf.gz
/path/chr2.vcf.gz
...
/path/chrM.vcf.gz
其它參數(shù)與 genotype-joint-calling 相同。文件列表中的 vcf 不需要事先合并，ilus VQSR 會(huì)幫你合并（如你已經(jīng)合并成一份了，也一樣使用），但你不要將兩個(gè)或者多個(gè)不同項(xiàng)目的 VCF 混在一起，那完全不是一回事，必須得是同一個(gè)項(xiàng)目的。以下提供一個(gè)完整的例子：
$ ilus VQSR -C ilus_sys.yaml -L vcf.list -O genotype --as_pipe_shell_order
其它
關(guān)于文中我提到的 yhbatch_slurm_jobs.py 和 check_jobs_status.py 都在 ilus github 代碼倉庫的 scripts 目錄里:
以上，就是關(guān)于 WGS/WES 流程生成器 ilus 的完整介紹，希望在有需要的時(shí)候它能成為你的一把小刀子，輔助你解決問題，使用過程中有任何問題都可以告知我或者在 github 上提給我。
訂閱
首發(fā)于個(gè)人公眾號(hào)：helixminer（堿基礦工）
