pexels-photo-262577.jpeg

不知覺間棋枕，距離我寫下第一篇關(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ā)狂的旅挤。

img

我在 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á)爾文生信星球) 中簡要回答過穗慕。

vqsr_method

我們?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 目錄里:

• https://github.com/ShujiaHuang/ilus/blob/master/scripts/yhbatch_slurm_jobs.py
• https://github.com/ShujiaHuang/ilus/blob/master/scripts/check_jobs_status.py

以上倔毙，就是關(guān)于 WGS/WES 流程生成器 ilus 的完整介紹埃仪，希望在有需要的時(shí)候它能成為你的一把小刀子，輔助你解決問題陕赃，使用過程中有任何問題都可以告知我或者在 github 上提給我卵蛉。

訂閱

首發(fā)于個(gè)人公眾號(hào)：helixminer（堿基礦工）