山海之间

illumina 肿瘤分析流程

illumina 肿瘤分析流程
2020-05-09 · 8 min read
生信 Illumina 流程

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

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

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

软件依赖

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

有一些脚本,在网上找不到下载地址,是直接从 Win 电脑上复制过来的,需要修改相应的XXX.runtimeconfig.json,将frameworkversion参数从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文件转化而来。
其中TargetStartTargetEnd,是Focus.dna_manifest.20171207.txt中是StartEnd减去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),欢迎关注!

Powered by Gridea,浙ICP备17039354号-1,© 2019 - 2020🍋