开箱即用版本 满分室间质评之GATK Somatic SNV+Indel+CNV+SV（20240430更新）

已于 2024-05-03 21:11:43 修改
阅读量330 收藏 1
点赞数
分类专栏： 生信分析流程 生信分析 图形化通用生信分析系统/平台开发 文章标签： 生物信息学 docker 生信分析流程 GATK python
于 2023-10-07 09:57:40 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/weixin_39900139/article/details/133632711
版权

本次更新可选参考基因组版本和项目bed文件，可以根据变量自行适配

最近准备为sliverworkspace 图形化生信平台开发报告设计器，需要一个较为复杂的pipeline作为测试数据，就想起来把之前的 满分室间质评之GATK Somatic SNV+Indel+CNV+SV（下）性能优化翻出来用一下。跑了一遍发现还是各种问题，于是想把pipeline改造成免部署、首次运行初始化环境的版本，以便需要时候能够直接运行起来，于是有了本文。

一句话描述就是：开箱即用的pipeline，能够根据样本version_reference自动选择参考基因组版本，根据project_bed文件选择项目bed，自动初始化环境、安装所需软件、下载ref文件和数据库的版本

为了让pipeline能够直接运行，无需部署，这里使用docker容器保证运行环境的一致性：见文章：基于docker的生信基础环境镜像构建，这里采用的方案是带ssh服务的docker+conda环境，整个pipeline在一个通用容器中运行。

本文代码较长，可能略复杂，想直接运行的可以下载 workflow文件，导入sliverworkspace图形化生信平台直接运行，获取并查看图形化分析报告。

相关代码和workflow文件可以访问笔者github项目地址

导入操作

在这里插入图片描述

分析流程整体概览

在这里插入图片描述

docker 镜像拉取及部署配置

# 拉取docker镜像
docker     pull     doujiangbaozi/sliverworkspace:1.10

docker-compose.yml配置文件

version: "3"
services:
  GATK:
    image: doujiangbaozi/sliverworkspace:1.10
    container_name: GATK
    volumes:
      - /home/sliver/Data/data:/opt/data:rw                                #挂载输入数据目录
      - /home/sliver/Manufacture/gatk/envs:/root/mambaforge-pypy3/envs     #挂载envs目录
      - /home/sliver/Manufacture/sliver/ref:/opt/ref:rw                    #挂载reference目录
      - /home/sliver/Manufacture/gatk/result:/opt/result:rw                #挂载中间文件和结果目录
    environment:
      - TZ=Asia/Shanghai                                                   #设置时区
      - PS=20191124                                                        #设置ssh密码
      - PT=9024                                                            #设置ssh连接端口

docker 镜像部署运行

# 在docker-compose.yml文件同级目录下运行
docker-compose up -d

# 或者docker-compose -f docker-compose.yml所在目录
docker-compose -f somedir/docker-compose.yml up -d

# 可以通过ssh连接到docker 运行pipeline命令了,连接端口和密码见docker-compose.yml配置文件相关字段
ssh root@127.0.0.1 -p9024

测试数据

测试数据来自2017年卫计委室间质评提供的bed文件（pipeline会自动下载）和测试数据，修改命名以匹配pipeline输入端，也可以替换为自己的数据文件，因为室间质评目前参考基因组还停留在hg19版本，所以本流程仍然使用hg19（GRCH37），如果要切换到hg38，可以将version_reference变量值设置为hg38，project_bed设置为Illumina_pt2_hg38.bed。pipeline会使用hg38（GRCH38）版本和对应的bed文件。

文件名（按照需要有调整）文件大小MD5
B1701_R1.fq.gz4.85G07d3cdccee41dbb3adf5d2e04ab28e5b
B1701_R2.fq.gz4.77Gc2aa4a8ab784c77423e821b9f7fb00a7
B1701NC_R1.fq.gz3.04G4fc21ad05f9ca8dc93d2749b8369891b
B1701NC_R2.fq.gz3.11Gbc64784f2591a27ceede1727136888b9

变量名称

# 变量初始化赋值
  sn=1701                               #样本编号
  pn=GS03                               #pipeline 代号
  project_bed=Illumina_pt2.bed          #参考基因组hg19下的bed文件
  version_reference=hg19                #人参考基因组版本hg19或者hg38
  version_fastqc=0.12.1                 #fastqc 版本
  version_multiqc=1.21                  #multiqc 版本
  version_cnvkit=0.9.10                 #cnvkit 版本
  version_manta=1.6.0                   #manta 版本
  version_gatk=4.3.0.0                  #gatk 版本                                 
  version_sambamba=1.0.1                #sambamba 版本                             
  version_samtools=1.17                 #samtools 版本                             
  version_bwa=0.7.17                    #bwa 版本                                  
  version_fastp=0.23.2                  #fastp 版本                                
  version_vep=108                       #vep 版本                                  
  envs=/root/mambaforge-pypy3/envs	    #mamba envs 目录                           
  threads=32                         	#最大线程数                                   
  memory=32G                        	#内存占用                                    
  scatter=8                          	#Mutect2 分拆并行运行interval list 份数          
  event=2                          	    #gatk FilterMutectCalls --max-events-in-region 数值
  snv_tlod=16.00                      	#snv 过滤参数 tload 值                        
  snv_vaf=0.01                       	#snv 过滤参数 丰度/突变频率                        
  snv_depth=500                        	#snv 过滤参数 支持reads数/depth 测序深度            
  cnv_dep=1000                       	#cnv 过滤参数 支持reads数/depth 测序深度            
  cnv_min=-0.5                       	#cnv 过滤参数 log2最小值                        
  cnv_max=0.5                        	#cnv 过滤参数 log2 最大值                       
  sv_score=200                        	#sv  过滤参数 score                           

## 以上变量个可以具体根据需求调整

表格：

变量名变量值备注
sn1701样本编号
pnGS03pipeline 代号 GATK Somatic 03版本
project_bedIllumina_pt2.bed参考基因组hg19下的bed文件
version_referencehg19人参考基因组版本hg19 / hg38
version_fastqc0.12.1fastqc 版本
version_multiqc1.21multiqc 版本
version_cnvkit0.9.10cnvkit 版本
version_manta1.6.0manta版本
version_gatk4.3.0.0gatk 版本
version_sambamba1.0.1sambamba 版本
version_samtools1.17samtools 版本
version_bwa0.7.17bwa 版本
version_fastp0.23.2fastp 版本
version_vep108vep 版本
envs/root/mambaforge-pypy3/envsmamba envs 目录
threads32最大线程数
memory32G内存占用
scatter8Mutect2 分拆并行运行interval list 份数
event2gatk FilterMutectCalls --max-events-in-region 数值
snv_tlod16.00snv 过滤参数 tload 值
snv_vaf0.01snv 过滤参数 丰度/突变频率
snv_depth500snv 过滤参数 支持reads数/depth 测序深度
cnv_dep1000cnv 过滤参数 支持reads数/depth 测序深度
cnv_min-0.5cnv 过滤参数 log2最小值
cnv_max0.5cnv 过滤参数 log2 最大值
sv_score200sv 过滤参数 score

Pipeline/workflow 具体步骤：

1. fastp 默认参数过滤，看下原始数据质量，clean data

#conda检测环境是否存在，首次运行不存在创建该环境并安装软件
if [ ! -d "${envs}/${pn}.fastp" ]; then
    echo "Creating the environment ${pn}.fastp"
    mamba create -n ${pn}.fastp -y fastp=${version_fastp} fastqc=${version_fastqc} multiqc=${version_multiqc}
fi

mamba	activate ${pn}.fastp

mkdir -p ${result}/${sn}/trimmed
mkdir -p ${result}/${sn}/qc

fastqc -t ${threads}\
	${data}/GS03/${sn}_tumor_R1.fq.gz \
	${data}/GS03/${sn}_tumor_R2.fq.gz \
	-o ${result}/${sn}/qc &

fastqc -t ${threads}\
	${data}/GS03/${sn}_normal_R1.fq.gz \
	${data}/GS03/${sn}_normal_R2.fq.gz \
	-o ${result}/${sn}/qc &

fastp -w 16 \
    -i ${data}/GS03/${sn}_tumor_R1.fq.gz  \
    -I ${data}/GS03/${sn}_tumor_R2.fq.gz  \
    -o ${result}/${sn}/trimmed/${sn}_tumor_R1_trimmed.fq.gz \
    -O ${result}/${sn}/trimmed/${sn}_tumor_R2_trimmed.fq.gz \
    -h ${result}/${sn}/trimmed/${sn}_tumor_fastp.html \
    -j ${result}/${sn}/trimmed/${sn}_tumor_fastp.json &

fastp -w 16 \
    -i ${data}/GS03/${sn}_normal_R1.fq.gz  \
    -I ${data}/GS03/${sn}_normal_R2.fq.gz  \
    -o ${result}/${sn}/trimmed/${sn}_normal_R1_trimmed.fq.gz \
    -O ${result}/${sn}/trimmed/${sn}_normal_R2_trimmed.fq.gz \
    -h ${result}/${sn}/trimmed/${sn}_normal_fastp.html \
    -j ${result}/${sn}/trimmed/${sn}_normal_fastp.json &
wait

fastqc -t ${threads}\
	${result}/${sn}/trimmed/${sn}_tumor_R1_trimmed.fq.gz \
	${result}/${sn}/trimmed/${sn}_tumor_R2_trimmed.fq.gz \
	-o ${result}/${sn}/qc &

fastqc -t ${threads}\
	${result}/${sn}/trimmed/${sn}_normal_R1_trimmed.fq.gz \
	${result}/${sn}/trimmed/${sn}_normal_R2_trimmed.fq.gz \
	-o ${result}/${sn}/qc &

wait

mamba	deactivate

2. normal文件fastq比对到参考基因组，sort 排序，mark duplicate 得到 marked.bam

#conda检测环境是否存在，首次运行不存在创建该环境并安装软件
if [ ! -d "${envs}/${pn}.align" ]; then
    mamba create -n ${pn}.align -y bwa=${version_bwa} samtools=${version_samtools} 
fi

#从github下载sambamba static 比 mamba 安装的版本速度快1倍以上.这是个很诡异的地方
if [ ! -f "${envs}/${pn}.align/bin/sambamba" ]; then
    aria2c https://github.com/biod/sambamba/releases/download/v${version_sambamba}/sambamba-${version_sambamba}-linux-amd64-static.gz -d ${envs}/${pn}.align/bin
    gzip -cdf ${envs}/${pn}.align/bin/sambamba-${version_sambamba}-linux-amd64-static.gz  >  ${envs}/${pn}.align/bin/sambamba 
    chmod a+x ${envs}/${pn}.align/bin/sambamba
fi

mamba	activate ${pn}.align

if [ "${version_reference}" == hg19 ]; then
    echo	"USE reference Version : ${version_reference}"
    mkdir	-p /opt/ref/hg19
    #如果没有检测到参考基因组序列，则下载序列并使用bwa创建索引
    if [ ! -f "/opt/ref/hg19/hg19.fasta" ]; then
    	aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg19/ucsc.hg19.fasta.gz -d /opt/ref/hg19 -o hg19.fasta.gz
    	cd /opt/ref/hg19 && gzip -d /opt/ref/hg19/hg19.fasta.gz
	else
        if [ -f "/opt/ref/hg19/ucsc.hg19.fasta.gz.aria2" ]; then
            aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg19/ucsc.hg19.fasta.gz -c -d /opt/ref/hg19 -o hg19.fasta.gz
            cd /opt/ref/hg19 && gzip -d /opt/ref/hg19/hg19.fasta.gz
        fi
	fi
	
    if  [ ! -f /opt/ref/hg19/hg19.fasta.amb ] ||
        [ ! -f /opt/ref/hg19/hg19.fasta.ann ] ||
        [ ! -f /opt/ref/hg19/hg19.fasta.bwt ] ||
        [ ! -f /opt/ref/hg19/hg19.fasta.pac ] ||
        [ ! -f /opt/ref/hg19/hg19.fasta.sa ]; then
        bwa index /opt/ref/hg19/hg19.fasta
    fi
    #检测samtools索引是否存在，如不存在则使用samtools创建参考基因组索引
    if [ ! -f "/opt/ref/hg19/hg19.fasta.fai" ]; then
    	samtools faidx /opt/ref/hg19/hg19.fasta
    fi
elif [ "${version_reference}" == "hg38" ]; then
    echo	"USE reference Version : ${version_reference}"
    mkdir	-p /opt/ref/hg38
#如果没有检测到参考基因组序列，则下载序列并使用bwa创建索引
    if [ ! -f "/opt/ref/hg38/hg38.fasta" ]; then
    	aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Homo_sapiens_assembly38.fasta.gz -d /opt/ref/hg38 -o hg38.fasta.gz
        cd /opt/ref/hg38 && gzip -d /opt/ref/hg38/hg38.fasta.gz
	else
        if [ -f "/opt/ref/hg38/hg38.fasta.gz.aria2" ]; then
            aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Homo_sapiens_assembly38.fasta.gz -d /opt/ref/hg38 -d /opt/ref/hg38 -o hg38.fasta.gz
            cd /opt/ref/hg38 && gzip -d /opt/ref/hg38/hg38.fasta.gz
        fi
    fi
    
    if  [ ! -f /opt/ref/hg38/hg38.fasta.amb ] ||
        [ ! -f /opt/ref/hg38/hg38.fasta.ann ] ||
        [ ! -f /opt/ref/hg38/hg38.fasta.bwt ] ||
        [ ! -f /opt/ref/hg38/hg38.fasta.pac ] ||
        [ ! -f /opt/ref/hg38/hg38.fasta.sa ]; then
        bwa index /opt/ref/hg38/hg38.fasta
    fi
    #检测samtools索引是否存在，如不存在则使用samtools创建参考基因组索引
    if [ ! -f "/opt/ref/hg38/hg38.fasta.fai" ]; then
    	samtools faidx /opt/ref/hg38/hg38.fasta
    fi
fi




mkdir -p ${result}/${sn}/aligned
#比对基因组管道输出成bam文件，管道输出排序
bwa mem \
    -t ${threads} -M \
    -R "@RG\\tID:${sn}_normal\\tLB:${sn}_normal\\tPL:Illumina\\tPU:Miseq\\tSM:${sn}_normal" \
    /opt/ref/${version_reference}/${version_reference}.fasta  ${result}/${sn}/trimmed/${sn}_normal_R1_trimmed.fq.gz ${result}/${sn}/trimmed/${sn}_normal_R2_trimmed.fq.gz \
    | sambamba view -S -f bam -l 0 /dev/stdin \
    | sambamba sort -t ${threads} -m 2G --tmpdir=${result}/${sn}/aligned -o ${result}/${sn}/aligned/${sn}_normal_sorted.bam /dev/stdin

#防止linux打开文件句柄数超过限制，报错退出
ulimit -n 10240

#使用sambamba对sorted bam文件标记重复
sambamba markdup \
    --tmpdir ${result}/${sn}/aligned \
    -t ${threads} ${result}/${sn}/aligned/${sn}_normal_sorted.bam ${result}/${sn}/aligned/${sn}_normal_marked.bam 

#修改marked bam文件索引名，gatk和sambamba索引文件名需要保持一致
mv  ${result}/${sn}/aligned/${sn}_normal_marked.bam.bai ${result}/${sn}/aligned/${sn}_normal_marked.bai
#删除sorted bam文件
rm -f ${result}/${sn}/aligned/${sn}_normal_sorted.bam*

mamba	deactivate

3. tumor文件fastq比对到参考基因组，排序，mark duplicate 得到 marked.bam，与第2步类似

if [ ! -d "${envs}/${pn}.align" ]; then
	mamba create -n ${pn}.align -y bwa=${version_bwa} samtools=${version_samtools} 
fi

#从github下载sambamba static 比 mamba 安装的版本速度快1倍以上.
if [ ! -f "${envs}/${pn}.align/bin/sambamba" ]; then
    aria2c https://github.com/biod/sambamba/releases/download/v${version_sambamba}/sambamba-${version_sambamba}-linux-amd64-static.gz -d ${envs}/${pn}.align/bin
    gzip -cdf ${envs}/${pn}.align/bin/sambamba-${version_sambamba}-linux-amd64-static.gz  >  ${envs}/${pn}.align/bin/sambamba 
    chmod a+x ${envs}/${pn}.align/bin/sambamba
fi

mamba	activate ${pn}.align

if [ "${version_reference}" == hg19 ]; then
    echo	"USE reference Version : ${version_reference}"
    mkdir	-p /opt/ref/hg19
	#如果没有检测到参考基因组序列，则下载序列并使用bwa创建索引
	if [ ! -f "/opt/ref/hg19/hg19.fasta" ]; then
        aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg19/ucsc.hg19.fasta.gz -d /opt/ref/hg19 -o hg19.fasta.gz
        cd /opt/ref/hg19 && gzip -d /opt/ref/hg19/hg19.fasta.gz
    else
        if [ -f "/opt/ref/hg19/ucsc.hg19.fasta.gz.aria2" ]; then
            aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg19/ucsc.hg19.fasta.gz -c -d /opt/ref/hg19 -o hg19.fasta.gz
            cd /opt/ref/hg19 && gzip -d /opt/ref/hg19/hg19.fasta.gz
        fi
	fi
	
    if  [ ! -f /opt/ref/hg19/hg19.fasta.amb ] ||
        [ ! -f /opt/ref/hg19/hg19.fasta.ann ] ||
        [ ! -f /opt/ref/hg19/hg19.fasta.bwt ] ||
        [ ! -f /opt/ref/hg19/hg19.fasta.pac ] ||
        [ ! -f /opt/ref/hg19/hg19.fasta.sa ]; then
        bwa index /opt/ref/hg19/hg19.fasta
    fi
    #检测samtools索引是否存在，如不存在则使用samtools创建参考基因组索引
    if [ ! -f "/opt/ref/hg19/hg19.fasta.fai" ]; then
    	samtools faidx /opt/ref/hg19/hg19.fasta
    fi
elif [ "${version_reference}" == "hg38" ]; then
    echo	"USE reference Version : ${version_reference}"
    mkdir	-p /opt/ref/hg38
    #如果没有检测到参考基因组序列，则下载序列并使用bwa创建索引
    if [ ! -f "/opt/ref/hg38/hg38.fasta" ]; then
    	aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Homo_sapiens_assembly38.fasta.gz -d /opt/ref/hg38 -o hg38.fasta.gz
        cd /opt/ref/hg38 && gzip -d /opt/ref/hg38/hg38.fasta.gz
	else
        if [ -f "/opt/ref/hg38/hg38.fasta.gz.aria2" ]; then
            aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Homo_sapiens_assembly38.fasta.gz -d /opt/ref/hg38 -d /opt/ref/hg38 -o hg38.fasta.gz
            cd /opt/ref/hg38 && gzip -d /opt/ref/hg38/hg38.fasta.gz
        fi
    fi
    
    if  [ ! -f /opt/ref/hg38/hg38.fasta.amb ] ||
        [ ! -f /opt/ref/hg38/hg38.fasta.ann ] ||
        [ ! -f /opt/ref/hg38/hg38.fasta.bwt ] ||
        [ ! -f /opt/ref/hg38/hg38.fasta.pac ] ||
        [ ! -f /opt/ref/hg38/hg38.fasta.sa ]; then
        bwa index /opt/ref/hg38/hg38.fasta
    fi
    #检测samtools索引是否存在，如不存在则使用samtools创建参考基因组索引
    if [ ! -f "/opt/ref/hg38/hg38.fasta.fai" ]; then
    	samtools faidx /opt/ref/hg38/hg38.fasta
    fi
fi

mkdir	-p ${result}/${sn}/aligned

bwa	mem \
    -t ${threads} -M \
    -R "@RG\\tID:${sn}_tumor\\tLB:${sn}_tumor\\tPL:Illumina\\tPU:Miseq\\tSM:${sn}_tumor" \
    /opt/ref/${version_reference}/${version_reference}.fasta  ${result}/${sn}/trimmed/${sn}_tumor_R1_trimmed.fq.gz ${result}/${sn}/trimmed/${sn}_tumor_R2_trimmed.fq.gz \
    | sambamba view -S -f bam -l 0 /dev/stdin \
    | sambamba sort -t ${threads} -m 2G --tmpdir=${result}/${sn}/aligned -o ${result}/${sn}/aligned/${sn}_tumor_sorted.bam /dev/stdin
ulimit -n 10240
sambamba  markdup \
    --tmpdir ${result}/${sn}/aligned \
    -t ${threads} ${result}/${sn}/aligned/${sn}_tumor_sorted.bam ${result}/${sn}/aligned/${sn}_tumor_marked.bam
    mv  ${result}/${sn}/aligned/${sn}_tumor_marked.bam.bai ${result}/${sn}/aligned/${sn}_tumor_marked.bai
    rm  -f ${result}/${sn}/aligned/${sn}_tumor_sorted.bam*

mamba	deactivate

​

4. 对上述bam文件生成重新校准表，为后续BQSR使用；Generates recalibration table for Base Quality Score Recalibration (BQSR)

#conda检测环境是否存在，首次运行不存在创建该环境并安装软件
if [ ! -f "${envs}/${pn}.gatk/bin/gatk" ]; then
    mamba create -n ${pn}.gatk -y gatk4=${version_gatk} \
        r-base=3.6.2 \
        r-data.table=1.12.8 \
        r-dplyr=0.8.5 \
        r-getopt=1.20.3 \
        r-ggplot2=3.3.0 \
        r-gplots=3.0.3 \
        r-gsalib=2.1 \
        r-optparse=1.6.4 \
        r-backports=1.1.10 \
        tabix
fi

mamba activate ${pn}.gatk

if [ "${version_reference}" == hg19 ]; then
    echo	"USE reference Version : ${version_reference}"
    #这里有个巨坑，broadinstitute ftp站点bundle目录提供的hg19版本参考文件，默认格式运行会报错，提示没有索引，使用gatk创建索引仍然报错，其实是gz格式需要使用bgzip重新压缩并且使用tabix创建索引才行
    if [ ! -f "/opt/ref/hg19/dbsnp_138.hg19.vcf.gz" ]; then
        aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg19/dbsnp_138.hg19.vcf.gz -d /opt/ref/hg19
    else 
        if [ -f "/opt/ref/hg19/dbsnp_138.hg19.vcf.gz.aria2" ]; then
            echo 'download continue...'
            aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg19/dbsnp_138.hg19.vcf.gz -c -d /opt/ref/hg19
        fi
    fi
    
    if [ ! -f "/opt/ref/hg19/dbsnp_138.hg19.vcf.gz.tbi" ]; then
        gzip -k -f -d /opt/ref/hg19/dbsnp_138.hg19.vcf.gz > /opt/ref/hg19/dbsnp_138.hg19.vcf
        bgzip -f --threads ${threads} /opt/ref/hg19/dbsnp_138.hg19.vcf > /opt/ref/hg19/dbsnp_138.hg19.vcf.gz
        tabix -f /opt/ref/hg19/dbsnp_138.hg19.vcf.gz
    fi

    if [ ! -f "/opt/ref/hg19/Mills_and_1000G_gold_standard.indels.hg19.sites.vcf.gz" ]; then
        aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg19/Mills_and_1000G_gold_standard.indels.hg19.sites.vcf.gz -d /opt/ref/hg19
    else 
        if [ -f "/opt/ref/hg19/Mills_and_1000G_gold_standard.indels.hg19.sites.vcf.gz.aria2" ]; then
            echo 'download continue...'
            aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg19/Mills_and_1000G_gold_standard.indels.hg19.sites.vcf.gz -c -d /opt/ref/hg19
        fi
    fi
    
    if [ ! -f "/opt/ref/hg19/Mills_and_1000G_gold_standard.indels.hg19.sites.vcf.gz.tbi" ]; then
        gzip -k -f -d /opt/ref/hg19/Mills_and_1000G_gold_standard.indels.hg19.sites.vcf.gz > /opt/ref/hg19/Mills_and_1000G_gold_standard.indels.hg19.sites.vcf
        bgzip -f --threads ${threads} /opt/ref/hg19/Mills_and_1000G_gold_standard.indels.hg19.sites.vcf > /opt/ref/hg19/Mills_and_1000G_gold_standard.indels.hg19.sites.vcf.gz
        tabix -f /opt/ref/hg19/Mills_and_1000G_gold_standard.indels.hg19.sites.vcf.gz
    fi

    if [ ! -f "/opt/ref/hg19/1000G_phase1.snps.high_confidence.hg19.sites.vcf.gz" ]; then
        aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg19/1000G_phase1.snps.high_confidence.hg19.sites.vcf.gz -d /opt/ref/hg19
    else 
        if [ -f "/opt/ref/hg19/1000G_phase1.snps.high_confidence.hg19.sites.vcf.gz.aria2" ]; then
            echo 'download continue...'
            aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg19/1000G_phase1.snps.high_confidence.hg19.sites.vcf.gz -c -d /opt/ref/hg19
        fi
    fi
    
    if [ ! -f "/opt/ref/hg19/1000G_phase1.snps.high_confidence.hg19.sites.vcf.gz.tbi" ]; then
        gzip -k -f -d /opt/ref/hg19/1000G_phase1.snps.high_confidence.hg19.sites.vcf.gz > /opt/ref/hg19/1000G_phase1.snps.high_confidence.hg19.sites.vcf
        bgzip -f --threads ${threads} /opt/ref/hg19/1000G_phase1.snps.high_confidence.hg19.sites.vcf > /opt/ref/hg19/1000G_phase1.snps.high_confidence.hg19.sites.vcf.gz
        tabix -f /opt/ref/hg19/1000G_phase1.snps.high_confidence.hg19.sites.vcf.gz
    fi

    if [ ! -f "/opt/ref/hg19/1000G_phase1.indels.hg19.sites.vcf.gz" ]; then
        aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg19/1000G_phase1.indels.hg19.sites.vcf.gz -d /opt/ref/
    else 
        if [ -f "/opt/ref/hg19/1000G_phase1.indels.hg19.sites.vcf.gz.aria2" ]; then
            echo 'download continue...'
            aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg19/1000G_phase1.indels.hg19.sites.vcf.gz -c -d /opt/ref/hg19
        fi
    fi
    
    if [ ! -f "/opt/ref/hg19/1000G_phase1.indels.hg19.sites.vcf.gz.tbi" ]; then
        gzip -k -f -d /opt/ref/hg19/1000G_phase1.indels.hg19.sites.vcf.gz > /opt/ref/hg19/1000G_phase1.indels.hg19.sites.vcf
        bgzip -f --threads ${threads} /opt/ref/hg19/1000G_phase1.indels.hg19.sites.vcf > /opt/ref/hg19/1000G_phase1.indels.hg19.sites.vcf.gz
        tabix -f /opt/ref/hg19/1000G_phase1.indels.hg19.sites.vcf.gz
    fi
    #创建参考序列hg19的dict字典文件
    if [ ! -f "/opt/ref/hg19/hg19.dict" ]; then
        gatk CreateSequenceDictionary -R /opt/ref/hg19/hg19.fasta -O /opt/ref/hg19/hg19.dict
    fi
    
    if [ ! -f "/opt/ref/${version_reference}/${project_bed}" ]; then
        #根据下载的Illumina_pt2.bed 文件创建interval list文件，坐标转换，起始坐标0修改为1
        #sed 's/chr//; s/\t/ /g' /opt/ref/hg19/Illumina_pt2.bed > /opt/ref/hg19/Illumina_pt2.processed.bed
        mkdir  -p /opt/ref/${version_reference}
		aria2c https://raw.githubusercontent.com/doujiangbaozi/sliverworkspace-util/main/somatic/projects/${project_bed} -d /opt/ref/hg19
        if [ ! -f "/opt/ref/hg19/${project_bed}.interval_list" ]; then
            gatk BedToIntervalList \
                -I	/opt/ref/hg19/${project_bed} \
                -SD	/opt/ref/hg19/hg19.dict \
                -O	/opt/ref/hg19/${project_bed}.interval_list
        fi
    fi
    


    mkdir -p ${result}/${sn}/recal

    gatk BaseRecalibrator \
            --known-sites /opt/ref/hg19/dbsnp_138.hg19.vcf.gz \
            --known-sites /opt/ref/hg19/Mills_and_1000G_gold_standard.indels.hg19.sites.vcf.gz \
            --known-sites /opt/ref/hg19/1000G_phase1.indels.hg19.sites.vcf.gz \
            -L /opt/ref/hg19/${project_bed}.interval_list \
            -R /opt/ref/hg19/hg19.fasta \
            -I ${result}/${sn}/aligned/${sn}_tumor_marked.bam \
            -O ${result}/${sn}/recal/${sn}_tumor_recal.table &

    gatk BaseRecalibrator \
            --known-sites /opt/ref/hg19/dbsnp_138.hg19.vcf.gz \
            --known-sites /opt/ref/hg19/Mills_and_1000G_gold_standard.indels.hg19.sites.vcf.gz \
            --known-sites /opt/ref/hg19/1000G_phase1.indels.hg19.sites.vcf.gz \
            -L /opt/ref/hg19/${project_bed}.interval_list \
            -R /opt/ref/hg19/hg19.fasta \
            -I ${result}/${sn}/aligned/${sn}_normal_marked.bam \
            -O ${result}/${sn}/recal/${sn}_normal_recal.table &

    wait
elif [ "${version_reference}" == "hg38" ]; then
    echo	"USE reference Version : ${version_reference}"
    if [ ! -f "/opt/ref/hg38/dbsnp_146.hg38.vcf.gz" ]; then
        aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/dbsnp_146.hg38.vcf.gz -d /opt/ref/hg38
    else 
        if [ -f "/opt/ref/hg38/dbsnp_146.hg38.vcf.gz.aria2" ]; then
            echo 'download continue...'
            aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/dbsnp_146.hg38.vcf.gz -c -d /opt/ref/hg38
        fi
    fi
    if [ ! -f "/opt/ref/hg38/dbsnp_146.hg38.vcf.gz.tbi" ]; then
        aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/dbsnp_146.hg38.vcf.gz.tbi -d /opt/ref/hg38
    else 
        if [ -f "/opt/ref/hg38/dbsnp_146.hg38.vcf.gz.tbi.aria2" ]; then
            echo 'download continue...'
            aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/dbsnp_146.hg38.vcf.gz.tbi -c -d /opt/ref/hg38
        fi
    fi

    if [ ! -f "/opt/ref/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz" ]; then
        aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz -d /opt/ref/hg38
    else 
        if [ -f "/opt/ref/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz.aria2" ]; then
            echo 'download continue...'
            aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz -c -d /opt/ref/hg38
        fi
    fi
    if [ ! -f "/opt/ref/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz.tbi" ]; then
        aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz.tbi -d /opt/ref/hg38
    else 
        if [ -f "/opt/ref/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz.tbi.aria2" ]; then
            echo 'download continue...'
            aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz.tbi -c -d /opt/ref/hg38
        fi
    fi

    if [ ! -f "/opt/ref/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz" ]; then
        aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz -d /opt/ref/hg38
    else 
        if [ -f "/opt/ref/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz.aria2" ]; then
            echo 'download continue...'
            aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz -c -d /opt/ref/hg38
        fi
    fi
    if [ ! -f "/opt/ref/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz.tbi" ]; then
        aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz.tbi -d /opt/ref/hg38
    else 
        if [ -f "/opt/ref/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz.tbi.aria2" ]; then
            echo 'download continue...'
            aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz.tbi -c -d /opt/ref/hg38
        fi
    fi
    #创建参考序列hg38的dict字典文件
    if [ ! -f "/opt/ref/hg38/hg38.dict" ]; then
        gatk CreateSequenceDictionary -R /opt/ref/hg38/hg38.fasta -O /opt/ref/hg38/hg38.dict
    fi
    
    if [ ! -f "/opt/ref/${version_reference}/${project_bed}" ]; then
        #根据下载的Illumina_pt2.bed 文件创建interval list文件，坐标转换，起始坐标0修改为1
        #sed 's/chr//; s/\t/ /g' /opt/ref/hg38/Illumina_pt2.bed > /opt/ref/hg38/Illumina_pt2.processed.bed
        mkdir  -p /opt/ref/hg38
        aria2c https://raw.githubusercontent.com/doujiangbaozi/sliverworkspace-util/main/somatic/projects/${project_bed} -d /opt/ref/hg38 -o ${project_bed}
        if [ ! -f "/opt/ref/hg38/${project_bed}.interval_list" ]; then
            gatk BedToIntervalList \
                -I	/opt/ref/hg38/${project_bed} \
                -SD	/opt/ref/hg38/hg38.dict \
                -O	/opt/ref/hg38/${project_bed}.interval_list
        fi
    fi
    
    mkdir -p ${result}/${sn}/recal

    gatk BaseRecalibrator \
            --known-sites /opt/ref/hg38/dbsnp_146.hg38.vcf.gz \
            --known-sites /opt/ref/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz \
            --known-sites /opt/ref/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz \
            -L /opt/ref/hg38/${project_bed}.interval_list \
            -R /opt/ref/hg38/hg38.fasta \
            -I ${result}/${sn}/aligned/${sn}_tumor_marked.bam \
            -O ${result}/${sn}/recal/${sn}_tumor_recal.table &

    gatk BaseRecalibrator \
            --known-sites /opt/ref/hg38/dbsnp_146.hg38.vcf.gz \
            --known-sites /opt/ref/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz \
            --known-sites /opt/ref/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz \
            -L /opt/ref/hg38/${project_bed}.interval_list \
            -R /opt/ref/hg38/hg38.fasta \
            -I ${result}/${sn}/aligned/${sn}_normal_marked.bam \
            -O ${result}/${sn}/recal/${sn}_normal_recal.table &

    wait
fi


​    
mamba deactivate

​

5. 使用校准表对bam碱基质量校准，因为这一步gatk效率感人，所以同时计算insertsize，拆分interval list（后续mutect2并行运行需要），运行cnvkit batch，运行samtools depth计算测序深度，samtools flagstat 统计mapping比例及质量

mkdir -p ${result}/${sn}/bqsr
mkdir -p ${result}/${sn}/stat
mkdir -p ${result}/${sn}/cnv
mkdir -p ${result}/${sn}/interval

mamba activate ${pn}.gatk

gatk ApplyBQSR \
    --bqsr-recal-file ${result}/${sn}/recal/${sn}_tumor_recal.table \
    -L /opt/ref/${version_reference}/${project_bed}.interval_list \
    -R /opt/ref/${version_reference}/${version_reference}.fasta \
    -I ${result}/${sn}/aligned/${sn}_tumor_marked.bam \
    -O ${result}/${sn}/bqsr/${sn}_tumor_bqsr.bam &


​    
gatk ApplyBQSR \
    --bqsr-recal-file ${result}/${sn}/recal/${sn}_normal_recal.table \
    -L /opt/ref/${version_reference}/${project_bed}.interval_list \
    -R /opt/ref/${version_reference}/${version_reference}.fasta \
    -I ${result}/${sn}/aligned/${sn}_normal_marked.bam \
    -O ${result}/${sn}/bqsr/${sn}_normal_bqsr.bam &


​    
gatk CollectInsertSizeMetrics \
    -I ${result}/${sn}/aligned/${sn}_tumor_marked.bam \
    -O ${result}/${sn}/stat/${sn}_tumor_insertsize_metrics.txt \
    -H ${result}/${sn}/stat/${sn}_tumor_insertsize_histogram.pdf &


​    
gatk CollectInsertSizeMetrics \
    -I ${result}/${sn}/aligned/${sn}_normal_marked.bam \
    -O ${result}/${sn}/stat/${sn}_normal_insertsize_metrics.txt \
    -H ${result}/${sn}/stat/${sn}_normal_insertsize_histogram.pdf &

rm -f ${result}/${sn}/interval/*.interval_list
gatk SplitIntervals \
    -L /opt/ref/${version_reference}/${project_bed}.interval_list \
    -R /opt/ref/${version_reference}/${version_reference}.fasta \
    -O ${result}/${sn}/interval \
    --scatter-count ${scatter} &

mamba deactivate

if [ ! -d "${envs}/${pn}.cnvkit" ]; then
    mamba create -n ${pn}.cnvkit -y cnvkit=${version_cnvkit}
fi

if [ ! -f "/opt/ref/${version_reference}/refFlat.txt" ]; then
    aria2c -x 16 -s 16 http://hgdownload.soe.ucsc.edu/goldenPath/${version_reference}/database/refFlat.txt.gz -d /opt/ref/${version_reference}
    cd /opt/ref/${version_reference} && gzip -d refFlat.txt.gz
fi

mamba activate ${pn}.cnvkit

rm -f ${result}/${sn}/cnv/${sn}_reference.cnn

cnvkit.py batch \
    ${result}/${sn}/aligned/${sn}_tumor_marked.bam \
    --normal ${result}/${sn}/aligned/${sn}_normal_marked.bam \
    --method hybrid \
    --fasta  /opt/ref/${version_reference}/${version_reference}.fasta \
    --targets /opt/ref/${version_reference}/${project_bed} \
    --annotate /opt/ref/${version_reference}/refFlat.txt \
    --output-reference ${result}/${sn}/cnv/${sn}_reference.cnn \
    --output-dir ${result}/${sn}/cnv/ \
    --diagram \
    -p ${threads} &

mamba deactivate

mamba activate ${pn}.align

samtools depth -a -b /opt/ref/${version_reference}/Illumina_pt2.bed  --threads ${threads} \
    ${result}/${sn}/aligned/${sn}_tumor_marked.bam > \
    ${result}/${sn}/stat/${sn}_tumor_marked.depth  &

samtools depth -a -b /opt/ref/${version_reference}/Illumina_pt2.bed  --threads ${threads} \
    ${result}/${sn}/aligned/${sn}_normal_marked.bam > \
    ${result}/${sn}/stat/${sn}_normal_marked.depth &

samtools flagstat --threads ${threads} \
    ${result}/${sn}/aligned/${sn}_tumor_marked.bam  > \
    ${result}/${sn}/stat/${sn}_tumor_marked.flagstat   &

samtools flagstat --threads ${threads} \
    ${result}/${sn}/aligned/${sn}_normal_marked.bam > \
    ${result}/${sn}/stat/${sn}_normal_marked.flagstat &
    
mamba deactivate

wait

​

6. 计算堆叠数据( pileup metrics )以便后续评估污染，也可以根据拆分的interval list并行处理，处理之后合并。

#官方巨坑,默认提供的small_exac_common_3_b37.vcf.gz默认染色体坐标不是以chr开头而是数字
mamba activate ${pn}.gatk

if [ "${version_reference}" == hg19 ]; then
    echo	"USE reference Version : ${version_reference}"
    #这里有个巨坑，从broadinstitute ftp 站点bundle Mutect2目录下载的参考文件，与同样下载的参考序列基因组坐标系不一致，参考基因组参考序列是chr1这种格式，这个af-only-gnomad是1，2，3这种格式，需要编写脚本处理
    if [ ! -f "/opt/ref/hg19/small_exac_common_3_b37.processed.vcf.gz" ]; then
        aria2c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/Mutect2/GetPileupSummaries/small_exac_common_3_b37.vcf.gz -d /opt/ref/hg19
    else
        if [ -f "/opt/ref/hg19/small_exac_common_3_b37.vcf.gz.aria2" ]; then
            aria2c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/Mutect2/GetPileupSummaries/small_exac_common_3_b37.vcf.gz -c -d /opt/ref/hg19
        fi
    fi
    if [ ! -f "/opt/ref/${version_reference}/small_exac_common_3_b37.processed.vcf.gz" ]; then
		if [ ! -f "${envs}/VcfProcessUtil.py" ]; then
			aria2c https://raw.fgit.cf/doujiangbaozi/sliverworkspace-util/main/somatic/VcfProcessUtil.py -d ${envs}/
			#aria2c https://raw.githubusercontent.com/doujiangbaozi/sliverworkspace-util/main/somatic/VcfProcessUtil.py -d ${envs}/
			chmod a+x ${envs}/VcfProcessUtil.py
		fi
		mamba activate ${pn}.cnvkit
		${envs}/VcfProcessUtil.py \
			-f /opt/ref/${version_reference}/small_exac_common_3_b37.vcf.gz \
			-o /opt/ref/${version_reference}/small_exac_common_3_b37.processed.vcf
		mamba deactivate
		mamba activate ${pn}.gatk
		cd /opt/ref/${version_reference}
		bgzip -f --threads ${threads} small_exac_common_3_b37.processed.vcf
		tabix -f small_exac_common_3_b37.processed.vcf.gz
		mamba deactivate
    fi
    mamba activate ${pn}.gatk
    for i in `ls ${result}/${sn}/interval/*.interval_list`;
    do
        echo $i
        gatk GetPileupSummaries \
            -R /opt/ref/${version_reference}/${version_reference}.fasta \
            -I ${result}/${sn}/bqsr/${sn}_tumor_bqsr.bam \
            -O ${i%.*}-tumor-pileups.table \
            -V /opt/ref/${version_reference}/small_exac_common_3_b37.processed.vcf.gz \
            -L $i \
            --interval-set-rule INTERSECTION &

        gatk GetPileupSummaries \
            -R /opt/ref/${version_reference}/${version_reference}.fasta \
            -I ${result}/${sn}/bqsr/${sn}_normal_bqsr.bam \
            -O ${i%.*}-normal-pileups.table \
            -V /opt/ref/${version_reference}/small_exac_common_3_b37.processed.vcf.gz \
            -L $i \
            --interval-set-rule INTERSECTION &

    done
    wait
    mamba deactivate
elif [ "${version_reference}" == "hg38" ]; then
    echo	"USE reference Version : ${version_reference}"
    if [ ! -f "/opt/ref/${version_reference}/small_exac_common_3.hg38.vcf.gz" ]; then
        aria2c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/Mutect2/GetPileupSummaries/small_exac_common_3.hg38.vcf.gz -d /opt/ref/${version_reference}
    else
        if [ -f "/opt/ref/hg19/small_exac_common_3_b37.vcf.gz.aria2" ]; then
            aria2c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/Mutect2/GetPileupSummaries/small_exac_common_3.hg38.vcf.gz -c -d /opt/ref/${version_reference}
        fi
    fi
    if [ ! -f "/opt/ref/${version_reference}/small_exac_common_3.hg38.vcf.gz.tbi" ]; then
        aria2c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/Mutect2/GetPileupSummaries/small_exac_common_3.hg38.vcf.gz.tbi -d /opt/ref/${version_reference}
    else
        if [ -f "/opt/ref/hg19/small_exac_common_3_b37.vcf.gz.tbi.aria2" ]; then
            aria2c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/Mutect2/GetPileupSummaries/small_exac_common_3.hg38.vcf.gz.tbi -c -d /opt/ref/${version_reference}
        fi
    fi
    mamba activate ${pn}.gatk
    for i in `ls ${result}/${sn}/interval/*.interval_list`;
    do
        echo $i
        gatk GetPileupSummaries \
            -R /opt/ref/${version_reference}/${version_reference}.fasta \
            -I ${result}/${sn}/bqsr/${sn}_tumor_bqsr.bam \
            -O ${i%.*}-tumor-pileups.table \
            -V /opt/ref/${version_reference}/small_exac_common_3.hg38.vcf.gz \
            -L $i \
            --interval-set-rule INTERSECTION &

        gatk GetPileupSummaries \
            -R /opt/ref/${version_reference}/${version_reference}.fasta \
            -I ${result}/${sn}/bqsr/${sn}_normal_bqsr.bam \
            -O ${i%.*}-normal-pileups.table \
            -V /opt/ref/${version_reference}/small_exac_common_3.hg38.vcf.gz \
            -L $i \
            --interval-set-rule INTERSECTION &

    done
    wait
    mamba deactivate
fi


mamba activate ${pn}.gatk
tables=
for i in `ls ${result}/${sn}/interval/*-tumor-pileups.table`;
do
    tables="$tables -I $i"
done

gatk GatherPileupSummaries \
    --sequence-dictionary /opt/ref/${version_reference}/${version_reference}.dict \
    $tables \
    -O ${result}/${sn}/stat/${sn}_tumor_pileups.table

nctables=
for i in `ls ${result}/${sn}/interval/*-normal-pileups.table`;
do
    nctables="$nctables -I $i"
done

gatk GatherPileupSummaries \
    --sequence-dictionary /opt/ref/${version_reference}/${version_reference}.dict \
    $nctables \
    -O ${result}/${sn}/stat/${sn}_normal_pileups.table
    
mamba deactivate

​

7. 使用GetPileupSummaries计算结果评估跨样本污染，结果用于后面 FilterMutectCall 过滤Mutect2输出结果

mamba activate ${pn}.gatk

gatk CalculateContamination \
    -matched ${result}/${sn}/stat/${sn}_normal_pileups.table \
    -I ${result}/${sn}/stat/${sn}_tumor_pileups.table \
    -O ${result}/${sn}/stat/${sn}_contamination.table
    
mamba deactivate

​

8. Mutect2 call 突变，使用拆分的interval list，结束后将结果合并；同时并行运行manta call sv突变

mkdir -p ${result}/${sn}/sv
mkdir -p ${result}/${sn}/snv

mamba activate ${pn}.gatk

if [ "${version_reference}" == hg19 ]; then
    echo	"USE reference Version : ${version_reference}"
    #这里有个巨坑，从broadinstitute ftp 站点bundle Mutect2目录下载的参考文件，与同样下载的参考序列基因组坐标系不一致，参考基因组参考序列是chr1这种格式，这个af-only-gnomad是1，2，3这种格式，需要编写脚本处理;hg38貌似没有这个问题，hg19的数据都不维护了么？
    if [ ! -f "/opt/ref/hg19/af-only-gnomad.raw.sites.b37.vcf.gz" ]; then
        aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/Mutect2/af-only-gnomad.raw.sites.b37.vcf.gz -d /opt/ref/hg19
    else
        if [ -f "/opt/ref/hg19/af-only-gnomad.raw.sites.b37.vcf.gz.aria2" ]; then
            aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/Mutect2/af-only-gnomad.raw.sites.b37.vcf.gz -c -d /opt/ref/hg19
        fi
        if [ ! -f "${envs}/VcfProcessUtil.py" ]; then
            aria2c https://raw.fgit.cf/doujiangbaozi/sliverworkspace-util/main/somatic/VcfProcessUtil.py -d ${envs}/
            #aria2c https://raw.githubusercontent.com/doujiangbaozi/sliverworkspace-util/main/somatic/VcfProcessUtil.py -d ${envs}/
            chmod a+x ${envs}/VcfProcessUtil.py
        fi
        mamba activate ${pn}.cnvkit
        ${envs}/VcfProcessUtil.py \
            -f /opt/ref/hg19/af-only-gnomad.raw.sites.b37.vcf.gz \
            -o /opt/ref/hg19/af-only-gnomad.raw.sites.b37.processed.vcf
        mamba deactivate
        mamba activate ${pn}.gatk
        cd /opt/ref/hg19
        bgzip -f --threads ${threads} af-only-gnomad.raw.sites.b37.processed.vcf
        tabix -f af-only-gnomad.raw.sites.b37.processed.vcf.gz
        mamba deactivate
    fi
elif [ "${version_reference}" == "hg38" ]; then
    echo	"USE reference Version : ${version_reference}"
    if [ ! -f "/opt/ref/${version_reference}/af-only-gnomad.hg38.vcf.gz" ]; then
        aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/Mutect2/af-only-gnomad.hg38.vcf.gz -d /opt/ref/${version_reference}
    else
        if [ -f "/opt/ref/${version_reference}/af-only-gnomad.hg38.vcf.gz.aria2" ]; then
            aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/Mutect2/af-only-gnomad.hg38.vcf.gz -c -d /opt/ref/${version_reference}
        fi
    fi
    if [ ! -f "/opt/ref/${version_reference}/af-only-gnomad.hg38.vcf.gz.tbi" ]; then
        aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/Mutect2/af-only-gnomad.hg38.vcf.gz.tbi -d /opt/ref/${version_reference}
    else
        if [ -f "/opt/ref/${version_reference}/af-only-gnomad.hg38.vcf.gz.tbi.aria2" ]; then
            aria2c -x 16 -s 32 ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/Mutect2/af-only-gnomad.hg38.vcf.gz.tbi -c -d /opt/ref/${version_reference}
        fi
    fi
fi


​    
mamba activate ${pn}.gatk
if [ ! -f "/opt/ref/${version_reference}/${project_bed}.gz" ]; then
    bgzip -f -c /opt/ref/${version_reference}/${project_bed} > /opt/ref/hg19/${project_bed}.gz
    tabix -f -p bed /opt/ref/${version_reference}/${project_bed}.gz
else
    if [ ! -f "/opt/ref/${version_reference}/${project_bed}.gz.tbi" ]; then
        tabix -f -p bed /opt/ref/${version_reference}/${project_bed}.gz
    fi
fi
mamba deactivate

if [ ! -d "${envs}/${pn}.manta" ]; then
    mamba create -n ${pn}.manta -y manta=1.6.0
fi

mamba activate ${pn}.manta

rm -f ${result}/${sn}/sv/runWorkflow.py*
configManta.py  \
    --normalBam ${result}/${sn}/bqsr/${sn}_normal_bqsr.bam \
    --tumorBam  ${result}/${sn}/bqsr/${sn}_tumor_bqsr.bam \
    --referenceFasta /opt/ref/${version_reference}/${version_reference}.fasta \
    --exome \
    --callRegions /opt/ref/${version_reference}/${project_bed}.gz \
    --runDir ${result}/${sn}/sv

rm -rf ${result}/${sn}/sv/workspace
python ${result}/${sn}/sv/runWorkflow.py -m local -j ${threads} &

mamba deactivate

mamba activate ${pn}.gatk

if [ "${version_reference}" == hg19 ]; then
    echo	"USE reference Version : ${version_reference}"
    rm -f ${result}/${sn}/snv/vcf-file.list
    touch ${result}/${sn}/snv/vcf-file.list
    for i in `ls ${result}/${sn}/interval/*.interval_list`;
    do
        rm -f ${i%.*}_bqsr.vcf.gz
        gatk Mutect2 \
            -R /opt/ref/${version_reference}/${version_reference}.fasta \
            -I ${result}/${sn}/bqsr/${sn}_tumor_bqsr.bam  -tumor  ${sn}_tumor  \
            -I ${result}/${sn}/bqsr/${sn}_normal_bqsr.bam -normal ${sn}_normal \
            -L $i \
            -O ${i%.*}_bqsr.vcf.gz \
            --max-mnp-distance 10 \
            --germline-resource /opt/ref/${version_reference}/af-only-gnomad.raw.sites.b37.processed.vcf.gz \
            --native-pair-hmm-threads ${threads} &
        echo ${i%.*}_bqsr.vcf.gz >> ${result}/${sn}/snv/vcf-file.list
    done
elif [ "${version_reference}" == "hg38" ]; then
    echo	"USE reference Version : ${version_reference}"
    rm -f ${result}/${sn}/snv/vcf-file.list
    touch ${result}/${sn}/snv/vcf-file.list
    for i in `ls ${result}/${sn}/interval/*.interval_list`;
    do
        rm -f ${i%.*}_bqsr.vcf.gz
        gatk Mutect2 \
            -R /opt/ref/${version_reference}/${version_reference}.fasta \
            -I ${result}/${sn}/bqsr/${sn}_tumor_bqsr.bam  -tumor  ${sn}_tumor  \
            -I ${result}/${sn}/bqsr/${sn}_normal_bqsr.bam -normal ${sn}_normal \
            -L $i \
            -O ${i%.*}_bqsr.vcf.gz \
            --max-mnp-distance 10 \
            --germline-resource /opt/ref/${version_reference}/af-only-gnomad.hg38.vcf.gz \
            --native-pair-hmm-threads ${threads} &
        echo ${i%.*}_bqsr.vcf.gz >> ${result}/${sn}/snv/vcf-file.list
    done
fi
wait

​    
rm -f ${result}/${sn}/snv/${sn}_bqsr.vcf.gz.stats
stats=
for z in `ls ${result}/${sn}/interval/*_bqsr.vcf.gz.stats`;
do
    stats="$stats -stats $z"
done

gatk MergeMutectStats $stats \
    -O ${result}/${sn}/snv/${sn}_bqsr.vcf.gz.stats

gatk MergeVcfs \
    -I ${result}/${sn}/snv/vcf-file.list \
    -O ${result}/${sn}/snv/${sn}_bqsr.vcf.gz
    
mamba deactivate

​

9. FilterMutectCalls 对Mutect结果突变过滤

mamba activate ${pn}.gatk

gatk FilterMutectCalls \
    --max-events-in-region ${event} \
    --contamination-table ${result}/${sn}/stat/${sn}_contamination.table \
    -R /opt/ref/${version_reference}/${version_reference}.fasta \
    -V ${result}/${sn}/snv/${sn}_bqsr.vcf.gz \
    -O ${result}/${sn}/snv/${sn}_filtered.vcf.gz
    
mamba deactivate

​

10. 使用Vep注释过滤结果

#conda检测环境是否存在，首次运行不存在创建该环境并安装软件
if [ ! -d "${envs}/${pn}.vep" ]; then
    echo "Creating the environment ${pn}.vep"
    mamba create -n ${pn}.vep -y ensembl-vep=${version_vep}
fi

mkdir -p /opt/result/${sn}/vcf

if [ "${version_reference}" == hg19 ]; then
echo	"USE reference Version : ${version_reference}"
    #检测vep注释数据库是否存在如果不存在则先下载
    if [ ! -d "/opt/ref/vep-cache/homo_sapiens/${version_vep}_GRCh37" ]; then
        aria2c -x 16 -s 48 https://ftp.ensembl.org/pub/release-${version_vep}/variation/indexed_vep_cache/homo_sapiens_vep_${version_vep}_GRCh37.tar.gz -d /opt/ref/
        tar -zxvf /opt/ref/homo_sapiens_vep_${version_vep}_GRCh37.tar.gz -C /opt/ref/vep-cache/
    elif [ -f "/opt/ref/homo_sapiens_vep_${version_vep}_GRCh37.tar.gz.aria2" ]; then
        aria2c -x 16 -s 48 https://ftp.ensembl.org/pub/release-${version_vep}/variation/indexed_vep_cache/homo_sapiens_vep_${version_vep}_GRCh37.tar.gz -c -d /opt/ref/
        tar -zxvf /opt/ref/homo_sapiens_vep_${version_vep}_GRCh37.tar.gz -C /opt/ref/vep-cache/
    fi

    if [ ! -d "/opt/ref/vep-cache/homo_sapiens_refseq/${version_vep}_GRCh37" ]; then
        aria2c -x 16 -s 48 http://ftp.ensembl.org/pub/release-${version_vep}/variation/vep/homo_sapiens_refseq_vep_${version_vep}_GRCh37.tar.gz -d /opt/ref/
        tar -zxvf /opt/ref/homo_sapiens_refseq_vep_${version_vep}_GRCh37.tar.gz -C /opt/ref/vep-cache/
    elif [ -f "/opt/ref/homo_sapiens_refseq_vep_${version_vep}_GRCh37.tar.gz.aria2" ]; then
        aria2c -x 16 -s 48 http://ftp.ensembl.org/pub/release-${version_vep}/variation/vep/homo_sapiens_refseq_vep_${version_vep}_GRCh37.tar.gz -c -d /opt/ref/
        tar -zxvf /opt/ref/homo_sapiens_refseq_vep_${version_vep}_GRCh37.tar.gz -C /opt/ref/vep-cache/
    fi
    
    mamba activate ${pn}.vep

    mkdir -p ${result}/${sn}/annotation
    vep \
        -i ${result}/${sn}/snv/${sn}_filtered.vcf.gz  \
        -o ${result}/${sn}/annotation/${sn}_filtered_vep.tsv \
        --offline \
        --cache \
        --cache_version ${version_vep} \
        --everything \
        --dir_cache /opt/ref/vep-cache/ \
        --dir_plugins /opt/ref/vep-cache/Plugins \
        --species homo_sapiens \
        --assembly GRCh37 \
        --fasta /opt/ref/${version_reference}/${version_reference}.fasta \
        --refseq \
        --force_overwrite \
        --format vcf \
        --tab \
        --shift_3prime 1  \
        --offline

    mamba deactivate
    
elif [ "${version_reference}" == "hg38" ]; then
    echo	"USE reference Version : ${version_reference}"
    #检测vep注释数据库是否存在如果不存在则先下载
    if [ ! -d "/opt/ref/vep-cache/homo_sapiens/${version_vep}_GRCh38" ]; then
        aria2c -x 16 -s 48 https://ftp.ensembl.org/pub/release-${version_vep}/variation/indexed_vep_cache/homo_sapiens_vep_${version_vep}_GRCh38.tar.gz -d /opt/ref/
        tar -zxvf /opt/ref/homo_sapiens_vep_${version_vep}_GRCh38.tar.gz -C /opt/ref/vep-cache/
    elif [ -f "/opt/ref/homo_sapiens_vep_${version_vep}_GRCh38.tar.gz.aria2" ]; then
        aria2c -x 16 -s 48 https://ftp.ensembl.org/pub/release-${version_vep}/variation/indexed_vep_cache/homo_sapiens_vep_${version_vep}_GRCh38.tar.gz -c -d /opt/ref/
        tar -zxvf /opt/ref/homo_sapiens_vep_${version_vep}_GRCh38.tar.gz -C /opt/ref/vep-cache/
    fi

    if [ ! -d "/opt/ref/vep-cache/homo_sapiens_refseq/${version_vep}_GRCh38" ]; then
        aria2c -x 16 -s 48 http://ftp.ensembl.org/pub/release-${version_vep}/variation/vep/homo_sapiens_refseq_vep_${version_vep}_GRCh38.tar.gz -d /opt/ref/
        tar -zxvf /opt/ref/homo_sapiens_refseq_vep_${version_vep}_GRCh38.tar.gz -C /opt/ref/vep-cache/
    elif [ -f "/opt/ref/homo_sapiens_refseq_vep_${version_vep}_GRCh38.tar.gz.aria2" ]; then
        aria2c -x 16 -s 48 http://ftp.ensembl.org/pub/release-${version_vep}/variation/vep/homo_sapiens_refseq_vep_${version_vep}_GRCh38.tar.gz -c -d /opt/ref/
        tar -zxvf /opt/ref/homo_sapiens_refseq_vep_${version_vep}_GRCh38.tar.gz -C /opt/ref/vep-cache/
    fi
    
    mamba activate ${pn}.vep

    mkdir -p ${result}/${sn}/annotation
    vep \
        -i ${result}/${sn}/snv/${sn}_filtered.vcf.gz  \
        -o ${result}/${sn}/annotation/${sn}_filtered_vep.tsv \
        --offline \
        --cache \
        --cache_version ${version_vep} \
        --everything \
        --dir_cache /opt/ref/vep-cache/ \
        --dir_plugins /opt/ref/vep-cache/Plugins \
        --species homo_sapiens \
        --assembly GRCh37 \
        --fasta /opt/ref/${version_reference}/${version_reference}.fasta \
        --refseq \
        --force_overwrite \
        --format vcf \
        --tab \
        --shift_3prime 1  \
        --offline

    mamba deactivate
fi

11. 使用脚本处理注释结果和过滤vcf结果，输出和室间质评要求格式的数据表格

mamba activate ${pn}.cnvkit

if [ ! -f "${envs}/MatchedSnvVepAnnotationFilter.py" ]; then
    aria2c https://raw.githubusercontent.com/doujiangbaozi/sliverworkspace-util/main/somatic/MatchedSnvVepAnnotationFilter.py -d ${envs}/
    chmod a+x ${envs}/MatchedSnvVepAnnotationFilter.py
fi

${envs}/MatchedSnvVepAnnotationFilter.py \
    -e normal_artifact   \
    -e germline   \
    -i strand_bias   \
    -i clustered_events   \
    --min-vaf=${snv_vaf} \
    --min-tlod=${snv_tlod} \
    --min-depth=${snv_depth} \
    -v ${result}/${sn}/snv/${sn}_filtered.vcf.gz   \
    -a ${result}/${sn}/annotation/${sn}_filtered_vep.tsv   \
    -o ${result}/${sn}/annotation/${sn}.result.SNV.tsv

mamba deactivate

​

12. 使用cnvkit提供工具输出分布图和热图

mamba activate ${pn}.cnvkit
    
    cnvkit.py scatter ${result}/${sn}/cnv/${sn}_tumor_marked.cnr \
    -s ${result}/${sn}/cnv/${sn}_tumor_marked.cns \
    -i ' ' \
    -n ${sn}_normal \
    -o ${result}/${sn}/cnv/${sn}_cnv_scatter.png -t  &&

cnvkit.py heatmap ${result}/${sn}/cnv/${sn}_tumor_marked.cns \
-o ${result}/${sn}/cnv/${sn}_cnv_heatmap.png

mamba deactivate

​

13. 使用cnvkit call 根据cnvkit batch输出结果推算拷贝数

mamba activate ${pn}.cnvkit

cnvkit.py call ${result}/${sn}/cnv/${sn}_tumor_marked.cns \
	-o ${result}/${sn}/cnv/${sn}_tumor_marked.call.cns

mamba deactivate

​

14. 编写脚本处理cnvkit输出，计算cnv基因，exon位置，gain/lost，cn数

mamba activate ${pn}.cnvkit

if [ ! -f "${envs}/CnvAnnotationFilter.py" ]; then
    aria2c https://raw.githubusercontent.com/doujiangbaozi/sliverworkspace-util/main/somatic/CnvAnnotationFilter.py -d ${envs}/
    chmod a+x ${envs}/CnvAnnotationFilter.py
fi

if [ ! -f "/opt/ref/${version_reference}/refGene.txt" ]; then
    aria2c -x 16 -s 16 http://hgdownload.cse.ucsc.edu/goldenPath/${version_reference}/database/refGene.txt.gz -d /opt/ref/${version_reference} -o refGene.txt.gz
    cd /opt/ref/${version_reference} && gzip -d refGene.txt.gz
fi

python ${envs}/CnvAnnotationFilter.py  \
    -r /opt/ref/${version_reference}/refGene.txt \
    -i ${cnv_min} \
    -x ${cnv_max} \
    -D ${cnv_dep} \
    -f ${result}/${sn}/cnv/${sn}_tumor_marked.call.cns \
    -o ${result}/${sn}/cnv/${sn}.result.CNV.tsv

mamba deactivate

​

15. 编写脚本处理manta的输出，获取最终sv输出结果，起始位置，基因、频率等

mamba activate ${pn}.cnvkit

if [ ! -f "${envs}/SvAnnotationFilter.py" ]; then
    aria2c https://raw.githubusercontent.com/doujiangbaozi/sliverworkspace-util/main/somatic/SvAnnotationFilter.py -d ${envs}/
    chmod a+x ${envs}/SvAnnotationFilter.py
fi

if [ ! -f "/opt/ref/${version_reference}/refGene.txt" ]; then
    aria2c -x 16 -s 16 http://hgdownload.cse.ucsc.edu/goldenPath/${version_reference}/database/refGene.txt.gz -d /opt/ref/${version_reference} -o refGene.txt.gz
    cd /opt/ref/${version_reference} && gzip -d refGene.txt.gz
fi

${envs}/SvAnnotationFilter.py \
    -r /opt/ref/${version_reference}/refGene.txt \
    -s ${sv_score} \
    -f ${result}/${sn}/sv/results/variants/somaticSV.vcf.gz \
    -o ${result}/${sn}/sv/${sn}.result.SV.tsv

mamba deactivate

​

16. 根据之前fastp,samtools depth,samtools flagstat,gatk CollectInsertSizeMetrics等输出，给出综合 QC数据

mamba activate ${pn}.cnvkit

if [ ! -f "${envs}/MatchedQcProcessor.py" ]; then
    aria2c https://raw.githubusercontent.com/doujiangbaozi/sliverworkspace-util/main/somatic/MatchedQcProcessor.py -d ${envs}/
    chmod a+x ${envs}/MatchedQcProcessor.py
fi

${envs}/MatchedQcProcessor.py  --bed /opt/ref/${version_reference}/${project_bed} \
    --out ${result}/${sn}/stat/${sn}.result.QC.tsv \
    --sample-fastp=${result}/${sn}/trimmed/${sn}_tumor_fastp.json \
    --sample-depth=${result}/${sn}/stat/${sn}_tumor_marked.depth \
    --sample-flagstat=${result}/${sn}/stat/${sn}_tumor_marked.flagstat \
    --sample-insertsize=${result}/${sn}/stat/${sn}_tumor_insertsize_metrics.txt \
    --normal-fastp=${result}/${sn}/trimmed/${sn}_normal_fastp.json \
    --normal-depth=${result}/${sn}/stat/${sn}_normal_marked.depth \
    --normal-flagstat=${result}/${sn}/stat/${sn}_normal_marked.flagstat  \
    --normal-insertsize=${result}/${sn}/stat/${sn}_normal_insertsize_metrics.txt

mamba deactivate

mamba activate ${pn}.fastp

multiqc ${result}/${sn}/ -f -o ${result}/${sn}/qc

mamba deactivate

​

最终输出

文件名备注
report.pdf由sliverworkspace设计的图形化分析报告打印输出为PDF
1701/qc/multiqc_report.htmlmultiqc聚合报告
1701/1701.result.SNV.tsvSNV最终突变结果
1701/1701/cnv/1701_cnv_heatmap.pngCNV结果热图
1701/cnv/1701_cnv_scatter.pngCNV结果分布图
1701/cnv/1701.result.CNV.tsvCNV最终结果
1701.result.SV.tsvSV最终结果
1701.result.QC.tsv最终质控结果
1701/qc/multiqc_report.htmlmultiqc报告
weixin_39900139
关注
  • 0
    点赞
  • 踩
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
专栏目录
ExAC:将 ExAC 数据导入管道
06-26
ExAC 来自 master 的变化 来自上游的变化，用于 saito test 创建分支基因组中心 ExAC 概述 该项目是为将 ExAC 数据导入管道而开发的。 检查已在 ExAC 的 VCF 文件中注册的每个样本的 SNV 位置。 关于 ExAC 数据库 麻省理工学院布罗德研究所。 链接： : 2014 年 10 月发布的 61,486 人（包括 4,378 名亚洲人）外显子组数据（0.2 版） 9,462,741 个变异位置 过滤器：DP >= 10 & GQ >= 20 2.9GB (gzip) aaaaaaaaaaaaa VCF 规范 (v4.1) 示例程序 (getExAC.pl) 读取每行的 VCFfile 导入 chr, pos, alt, AC, AC_EAS, AN, AN_EAS 到 hash-table (key: chr, pos) 从每个
GATK-HaplotypeCaller 变异检测详解
qq_28723681的博客
03-18 9343
1.GATK-HaplotypeCaller简介 基因组变异检测是基因组学领域一个非常重要的问题，是遗传性疾病溯源，物种进化等分析的前提。而目前最主流、使用最广泛的变异检测软件当属 Broad Institute 开发的 GATK(Genome Analysis ToolKit) 组件。GATK 设计之初是用于分析人类的全外显子和全基因组数据，随着不断发展，现在也可以应用于其他的物种。GATK官网提供了一整套完整的变异检测分析流程：GATK Best Practices。如下图示： 其中，Haplotyp
[更新中] NGS常见软件和数据库（肿瘤/遗传病）
Rys_Ben_Blogs
02-27 1691
变异分析
满分室间质评之GATK Somatic SNV+Indels+CNV+SV（上）
SliverWorkspace
07-31 3241
卫计委在2017年，2019年，2020年（还没有答案）提供标准数据用于肿瘤生信分析的室间质评。这样预知结果的数据自然是不能放过了，本文尝试参考GATK Best Practice：Somatic SNVs + Indels ，Cnvkit，Manta的pipeline来完成满分流程分析，也可以使用标准数据反向判断GATK Mutect2的实际准确度,算法优劣。 注：本文仅用于学习，距离真正的临床应用还有相当大距离，欢迎大佬批评指正 **1. 分析流程概览如下： 2. 本文用到的分析系统及分析流程文件
满分室间质评之GATK Somatic SNV+Indel+CNV+SV（下）性能优化
SliverWorkspace
08-04 2942
我们接上文：满分室间质评之GATK Somatic SNV+Indel+CNV+SV一文中实现了对于卫计委室间质评数据分析以及与满分结果的匹配。本文将着重解决，保证最终结果一致的情况下，如何优化分析性能（并行化），如何将分析时间从 3h 59m 53s缩短至 1h 10m 38s。 优化的方向：实际运行GATK4.X的工具如Mutect2时，发现其运行效率相当低，从CPU占用率，内存占用，硬盘I/O都占用很低，起初自己DIY时候，将要分析的bed/interval_list文件按照染色体编号拆分（不太确定
GenomicsDB：C ++中的高性能数据存储，用于导入，查询和转换带有CC ++ JavaSpark绑定的变体数据。 在gatk4中使用
02-05
变异数据本质上是稀疏的（相对于整个基因组稀疏），并且使用稀疏数组数据存储非常适合存储此类数据。 GenomicsDB是用C ++编写的高性能可扩展数据存储，用于导入，查询和转换基因组变异数据。 支持的平台：Linux和...
gatk-sv
03-19
GATK-SV 用于Illumina短读全基因组测序（WGS）数据的结构变异发现管道。 目录 原始呼叫者和证据收集 批量质量控制 -gCNV模型创建 批量证据合并，BAF生成和深度调用者 网站集群 网站指标 过滤 跨批次站点合并 基因...
gatk:GATK版本4及更高版本的官方代码存储库
05-11
GATK4旨在在简化的框架下将来自和代码库的完善工具集中在一起，并使选定的工具能够使用在本地集群或云中大规模并行地运行。 它还包含许多新开发的工具，该工具包的早期版本中没有这些工具。 目录 将JVM选项传递给...
gatk_recalibrate_and_realign:GATK Indel校正和质量校正
05-07
此应用程序对Mappings对象应用了许多重新计算，以提高变量调用的质量。 该过程遵循Broad Institute关于变式调用最佳实践的建议。 映射在已知插入缺失的位点周围进行重复数据删除和重新排列，并通过查看质量指标与...
gatk, GATK版本 4和更高版本的官方代码库.zip
09-18
gatk, GATK版本 4和更高版本的官方代码库 注意：请查看网站，你可以下载预编译可执行文件，阅读文档，提问，并接收技术支持。 GATK 4这个库包含了下一代基因组分析工具包( GATK ) 。 这个存储库的内容是 100%开放...
图形化开放式生信分析系统开发 - 9 Illumina测序仪测序数据自动拆分
SliverWorkspace
12-05 2164
前文链接： 图形化开放式生信分析系统开发 - 1 需求分析及技术实现 图形化开放式生信分析系统开发 - 2 样本信息处理 图形化开放式生信分析系统开发 - 3 生信分析流程的进化 图形化开放式生信分析系统开发 - 4 生信分析流程的图形化 图形化开放式生信分析系统开发 - 5 生信分析流程服务器端运行 图形化开放式生信分析系统开发 - 6 生信分析流程批量运行与过程控制 图形化开放式生信分析系统...
新冠病毒分型和突变分析（SARS-CoV2_ARTIC_Nanopore）
SliverWorkspace
11-30 1199
SARS-Cov2 新冠病毒分析 Nanopore Artic 生信分析流程 新冠病毒分型 新冠病毒突变分析
使用docker完成生信分析环境搭建
SliverWorkspace
12-12 1179
生信开发人员最头疼的问题，可能就是平台搭建和软件安装了。部署和迁移上要费很大力气。本文讲述使用docker制作一个镜像，后续通过导入自己定制的镜像，复制文件完成分析流程的部署和迁移。 如何使用docker，推荐阅读 Docker — 从入门到实践 一、准备工作 首先我们需要对手头的pipeline做一个分类。为了保持docker镜像尽可能的小，不能将所有文件全部放在docker镜像里面，需要外...
图形化开放式生信分析系统开发 - 3 生信分析流程的进化
SliverWorkspace
12-02 1171
接上两篇内容，本文主要讲述工作中NGS从科研进入医学临床领域，工作中接触到生信流程，以及最终在自动图形化开放式生信分析系统开发中生信workflow设计实现的过程。 接触二代测序，生信分析，那真是打开了一个新世界的大门，各种名次术语满天飞，搞的头晕脑胀。什么“什么是高通量测序/NGS”、Sanger法测序（一代测序）、外显子测序（whole exon sequencing）、mRNA测序 （...
图形化开放式生信分析系统开发 - 4 生信分析流程图形化设计
SliverWorkspace
12-02 1093
前文链接： 自动化图形生物信息分析系统开发 - 1 需求分析及技术实现 自动化图形生物信息分析系统开发 - 2 样本信息处理 自动图形化开放式生信分析系统开发 - 3 生信分析pipeline的进化 在上文自动图形化开放式生信分析系统开发 - 3 生信分析pipeline的进化 讨论了生信分析pipeline的进化，从手动到自动，但仍然停留在终端命令行阶段，为了让更多非生信专业的人能够使用...
新冠病毒分型和突变分析（SARS-CoV2_ARTIC_Illumina）
SliverWorkspace
11-09 1017
SARS-CoV2 SARS-CoV-2 Illumina 新冠病毒分型和突变分析
使用宏基因组的方法快速鉴定新冠病毒SARS-CoV2
SliverWorkspace
12-07 899
使用宏基因组的方法快速鉴定新冠病毒SARS-CoV2
图形化开放式生信分析系统开发 - 2样本信息处理
SliverWorkspace
12-02 726
一、基于生信分析云平台的需求，在下图中可以看出样本信息在整个软件中的作用 样本信息用于数据拆分、运行运行流程、分析报告，是整个系统的基础数据之一。 下面来具体归纳一下样本信息在软件系统的作用并列出详细的数据字段： 1、用于标记分析运行状态的 序号 字段名称 数据类型 作用 1 ...
spark gatk
最新发布
01-24
Spark GATK是指在Spark集群上运行GATK（Genome Analysis Toolkit）工具。GATK是一个广泛使用的工具集，用于分析基因组数据。它提供了一系列的工具，用于变异检测、变异过滤、变异注释等。 在Spark集群上运行GATK可以显著加速分析流程的执行速度，因为Spark具有分布式计算的能力，可以并行处理大规模的基因组数据。通过将GATK与Spark集成，可以充分利用集群中的多个计算节点，提高数据处理的效率。 要在Spark集群上运行GATK，您可以按照以下步骤进行操作： 1. 首先，确保您已经安装了Spark和GATK，并且集群已经正确配置。 2. 将GATK工具和相关的输入数据上传到集群中。 3. 使用Spark的API或命令行工具来提交GATK任务。您可以使用Spark的分布式计算能力来处理大规模的数据集。 4. 监控任务的执行情况，并根据需要进行调优和优化。 请注意，并非所有的GATK工具都使用Spark。具体使用哪些工具取决于您的需求和数据集的规模。您可以参考GATK的文档和技术支持资源，了解更多关于在Spark集群上运行GATK的详细信息。

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值