illumina 肿瘤分析流程

本地 Win 电脑部署了一套 ILLUMINA 的肿瘤分析流程，有 DNA 和 RNA 两个流程。需要将这套分析流程部署到 Linux 服务器上。

通过查看日志，发现 DNA 的 Workflow 主要有以下几个部分：

数据拆分
- 数据拆分：通过 bcl2fastq 拆分数据（现有流程已有拆分步骤，可忽略）；
bam 文件处理
- bwa比对：将 fastq 文件比对到基因组；
- TranslateBam：作用未知，似乎是通过 target 文件对 bam 文件区域做了限定；
- indel重排序：对 bam 文件中的 indel 进行重新对齐，不清楚是左对齐还是右对齐；
- 合并reads：合并 overlap paired-end read；
VariantCalling
- call 变异：使用 Pisces 进行 call 变异；
- 突变质量重矫正：对结果 vcf 结果进行重矫正
结果注释及统计
- 结果注释：对 vcf 结果进行数据库注释
- bam文件统计：对 bam 文件进行统计，但是使用的 puma.exe 无法在 Linux 上运行。也找不到相关资料
- 突变结果统计：对最终 vcf 结果进行统计，但并没有结果文件生成

软件依赖

流程脚本基本是 dll 脚本，可以通过安装.NET Core来跨平台运行。

有一些脚本，在网上找不到下载地址，是直接从 Win 电脑上复制过来的，需要修改相应的XXX.runtimeconfig.json，将framework的version参数从1.0.0改成自己安装的.NET Core版本，比如我现在装的是2.1.7。

另外针对所有的dll脚本，需要在XXX.runtimeconfig.json添加一项configProperties配置，具体如下：

{"runtimeOptions": {"tfm": "netcoreapp2.0","framework": {"name": "Microsoft.NETCore.App","version": "2.0.0" },"configProperties": {                    #这部分配置"System.Globalization.Invariant": true #需自己添加}}
}

否则运行的时候，会报一个关于 ICU package的错误，并退出。

基因组，注释数据库等相关资源文件，都是从 Win 电脑上复制的，版本较旧。也可以从 ILLUMINA 网站重新下载。

.NET Core下载地址

bam 文件处理

通过 bwa 比对 fastq 文件到基因组，产生 bam 文件，并对 bam 文件相应的处理。

bwa 比对

在本地 Win 电脑上部署的流程，拆分的时候会将一个样本的数据拆成 4 个lane。所有一个样本会有 4 对 R1R2 文件，导致比对时生成 4 个 bam 文件，最后通过 samtool cat命令合并成一个 bam 文件。

在 Linux 服务器上，拆分后只生成 1 对 R1R2 文件，直接使用 bwa 比对即可，不需要后续合并步骤。

这部分命令如下：

bwa mem -M -t 4 -R "@RG\tID:LC10-tz2020001DNA\tPL:ILLUMINA\tSM:LC10-tz2020001DNA" -L 50 -Y ../hg19/Sequence/WholeGenomeFasta/genome.fa LC10-tz2020001DNA_S1_R1_001.fastq.gz LC10-tz2020001DNA_S1_R2_001.fastq.gz | samtools view -Sb - > LC10-tz2020002DNA.bam

参数按照 Win 版流程的参数。

TranslateBam

脚本： TranslateBam.dll ，该脚本从 Win 电脑复制，网上无法下载；
作用：未知，可能是对 bam 文件进行限定。不跑这一步，下一步运行时间会非常久；
命令：

## TranslateBam
~/dotnet/dotnet ~/TranslateBam/TranslateBam.dll ../Focus.dna_manifest.20171207.txt.targets.txt ../hg19/Sequence/WholeGenomeFasta/GenomeSize.xml LC10-tz2020001DNA.bam LC10-tz2020001DNA_translate.bam LC10-tz2020001DNA.TargetHits.tsv 0.04 0.1 30 1## bam文件排序
samtools sort -@ 8 LC10-tz2020001DNA_translate.bam -o LC10-tz2020001DNA_translate.sorted.bam# 建立索引
samtools index LC10-tz2020001DNA_translate.sorted.bam

其中的Focus.dna_manifest.20171207.txt.targets.txt文件是由原始的Focus.dna_manifest.20171207.txt文件转化而来。
其中TargetStart和TargetEnd，是Focus.dna_manifest.20171207.txt中是Start、End减去prob sequence的长度。
这个文件可以作为固定的文件。

Indel重排序

脚本：Hygea.dll，从 Github 下载最新版；
作用：Indel重排序；
命令：

~/dotnet/dotnet ~/Hygea_5.2.10.49/Hygea.dll -b Alignment/LC10-tz2020001DNA_translate.sorted.bam -genomeFolders hg19/Sequence/WholeGenomeFasta/ -useAlignmentScorer True -softclipCoefficient 0 -outFolder hygea

合并reads

脚本：Stitcher.dll，从 Github 下载最新版；
作用：合并reads；
命令：

~/dotnet/dotnet ~/Stitcher_5.2.10.49/Stitcher.dll --bam hygea/LC10-tz2020001DNA_translate.sorted.bam --OutFolder stitcher --MinBaseCallQuality 20 --minmapquality 1 --FilterDuplicates True --FilterForProperPairs False --FilterUnstitchablePairs False --StitchGappedPairs False --UseSoftClippedBases True --NifyUnstitchablePairs True --NumThreads 1## 排序
samtools sort -@ 8 LC10-tz2020001DNA_translate.sorted.stitched.bam -o LC10-tz2020001DNA.stitched.sorted.bam## 建立索引
samtools index LC10-tz2020001DNA.stitched.sorted.bam

NumThreads 这个线程参数，好像得设置为 1 。测试的时候，设置为 8 ，但是程序一直不结束。当设置为 1 时，可以较快结束。未过多测试，不知道是否是 BUG。

VariantCalling

使用Pisces软件进行call变异，适用与多重PCR的数据，并推荐用在tumor-only的流程中。

call 变异

命令：Pisces.dll，从 Github 下载最新版；
作用：call 变异；
命令：

~/dotnet/dotnet ~/Pisces_5.2.10.49/Pisces.dll -b stitcher/LC10-tz2020001DNA.stitched.sorted.bam -g hg19/Sequence/WholeGenomeFasta -gVCF True -OutFolder Pisces -i interval/intervals.picard -CallMNVs True -MaxMNVLength 3 -MaxGapBetweenMNV 1 -MinimumFrequency 0.01 -MinBaseCallQuality 20 -MaxVariantQScore 100 -MinVariantQScore 20 -CrushVcf False -MinDP 10 -ReportNoCalls True -Ploidy somatic -EnableSingleStrandFilter False -MinDPFilter 10 -MaxAcceptableStrandBiasFilter 0.5 -VariantQualityFilter 30 -MinVariantFrequencyFilter 0.02 -RMxNFilter 3,6,0.20 -MaxNumThreads 1

其中-i参数，其实就是平时流程中的bed文件。这个文件也可以作为固定文件，不需要每次生成。

质量重矫正

命令：VariantQualityRecalibration.dll，从 Github 下载最新版；
作用：结果 vcf 文件质量重矫正；
命令：

~/dotnet/dotnet ~/VariantQualityRecalibration_5.2.10.49/VariantQualityRecalibration.dll -vcf Pisces/LC10-tz2020001DNA.stitched.sorted.genome.vcf -o VariantQualityRecalibration -b 20 -z 3 -Q 100 -f 30## 压缩vcf文件
/work/bin/bgzip -f Pisces/LC10-tz2020001DNA.stitched.sorted.genome.vcf## 建立vcf文件索引
/work/bin/tabix -p vcf -f Pisces/LC10-tz2020001DNA.stitched.sorted.genome.vcf.gz

质量重矫正没有新文件生成，似乎没什么作用？应该与每个样本的结果有关。

注释及统计

使用Nirvana对突变结果进行注释。
Illumina 的官方 Github 仓库提供一键安装的脚本，会下载需要的数据库资源文件。
因为数据量非常大(30G)，所以直接使用的是 Win 电脑里的旧版本脚本和数据。

结果注释

脚本：Nirvana.dll，本地 Win 电脑复制；
作用：对突变进行注释；
命令：

~/dotnet/dotnet ~/Nirvana/Nirvana.dll --in Pisces/LC10-tz2020001DNA.stitched.sorted.genome.vcf.gz --sd hg19/Annotation/Nirvana/SupplementaryDatabase/38 --out nirvana/LC10-tz2020001DNA --ref hg19/Annotation/Nirvana/Reference/5/Reference.dat --cache hg19/Annotation/Nirvana/Cache/24/RefSeq84 --gvcf --vcf## 建立索引/work/bin/tabix -p vcf -f nirvana/LC10-tz2020001DNA.vcf.gz
/work/bin/tabix -p vcf -f nirvana/LC10-tz2020001DNA.genome.vcf.gz

--out参数，会将路径的最后一部分作为输出文件前缀。如上述命令，会在nirvana路径下生成LC10-tz2020001DNA开头的文件。

bam文件统计

这部分在 Win 电脑上使用的是 puma.exe 软件，无法在 Linux 上运行。并且在网上也找不到相关资料，可以使用 bamdst 软件来代替。

突变结果统计

脚本：VariantStatsGenerator.dll，本地 Win 电脑上复制；
作用：应该是突变统计？但是没有结果文件生成，只有日志文件；
命令：

~/dotnet/dotnet ~/VariantStatsGenerator/VariantStatsGenerator.dll -m samllVar -w single -t nirvana/LC10-tz2020001DNA.vcf.gz StatisticsEvaluation

最后一个参数是输出目录。

总结

如果要部署 ILLUMINA 的流程，可以只到 call变异部分即可。后续的注释部分，可以使用自己的流程。

在 Pisces 的Github仓库上有更新的流程，会用到新的软件脚本，也可以作为参考。

本文首发于公众号：柠檬培养师（ID: yantinger90），欢迎关注！