sentieon使用

  • gvcf_joincalling

vi test.sh
ls *.gz >input.setieon.list
export SENTIEON_LICENSE=192.168.102.80:8990
cat input.setieon.list | /home/yzhou/soft/sentieon-genomics-202112.07/bin/sentieon driver -r /home/yzhou/refs/pig/Sus_scrofa.Sscrofa11.1-106/Sus_scrofa.Sscrofa11.1.dna.toplevel.fa -t 20 --algo GVCFtyper -d /home/yzhou/refs/pig/Sus_scrofa.Sscrofa11.1-106/sus_scrofa.vcf.gz /data1/yzhou/ZZQ2003/typedvcf/combZZQ2003_chr13_typed.vcf.gz -
  • 单个个体SNPcalling

cd /data1/data_public/WGS/pig/XinJiang382/2.5G/rawdata/cleandata
mkdir DNAseq
vi  test01.sh
#!/bin/sh
# Copyright (c) 2016-2020 Sentieon Inc. All rights reserved
# *******************************************
# Script to perform DNA seq variant calling
# using a single sample with fastq files
# named 1.fastq.gz and 2.fastq.gz
# *******************************************
set -eu
# Update with the fullpath location of your sample fastq
SM="test01" #sample name
RGID="rg_$SM" #read group ID
PL="ILLUMINA" #or other sequencing platform
FASTQ_FOLDER="/data1/data_public/WGS/pig/XinJiang382/2.5G/rawdata/cleandata"
FASTQ_1="$FASTQ_FOLDER/${SM}_1_val_1.fq.gz"
FASTQ_2="$FASTQ_FOLDER/${SM}_2_val_2.fq.gz"  #If using Illumina paired data
# Update with the location of the reference data files
FASTA_DIR="/home/jychu/home/refs/pig/Sus_scrofa.Sscrofa11.1-106"
FASTA="$FASTA_DIR/Sus_scrofa.Sscrofa11.1.dna.toplevel.fa"
KNOWN_DBSNP="$FASTA_DIR/sus_scrofa.vcf.gz"
# Update with the location of the Sentieon software package and license file
SENTIEON_INSTALL_DIR=/home/yzhou/soft/sentieon-genomics-202112.07
export SENTIEON_LICENSE=192.168.102.80:8990  #or using licsrvr: c1n11.sentieon.com:5443
# Other settings
NT=1  #number of threads to use in computation, set to number of cores in the server
START_DIR="$FASTQ_FOLDER/DNAseq" #Determine where the output files will be stored
# You do not need to modify any of the lines below unless you want to tweak the pipeline
# ************************************************************************************************************************************************************************
# ******************************************
# 0. Setup
# ******************************************
WORKDIR="$START_DIR/${SM}" 
mkdir -p $WORKDIR
LOGFILE=$WORKDIR/run.log
exec >$LOGFILE 2>&1
cd $WORKDIR
# ******************************************
# 1. Mapping reads with BWA-MEM, sorting
# ******************************************
#The results of this call are dependent on the number of threads used. To have number of threads independent results, add chunk size option -K 10000000 
( $SENTIEON_INSTALL_DIR/bin/sentieon bwa mem -R "@RG\tID:$RGID\tSM:$SM\tPL:$PL" \
    -t $NT -K 10000000 $FASTA $FASTQ_1 $FASTQ_2 || { echo -n 'BWA error'; exit 1; } ) | \
    $SENTIEON_INSTALL_DIR/bin/sentieon util sort -r $FASTA -o sorted.bam -t $NT \
    --sam2bam -i - || { echo "Alignment failed"; exit 1; }
# ******************************************
# 2. Metrics
# ******************************************
$SENTIEON_INSTALL_DIR/bin/sentieon driver -r $FASTA -t $NT -i sorted.bam \
    --algo MeanQualityByCycle mq_metrics.txt --algo QualDistribution qd_metrics.txt \
    --algo GCBias --summary gc_summary.txt gc_metrics.txt --algo AlignmentStat \
    --adapter_seq '' aln_metrics.txt --algo InsertSizeMetricAlgo is_metrics.txt || \
    { echo "Metrics failed"; exit 1; }

$SENTIEON_INSTALL_DIR/bin/sentieon plot GCBias -o gc-report.pdf gc_metrics.txt
$SENTIEON_INSTALL_DIR/bin/sentieon plot QualDistribution -o qd-report.pdf qd_metrics.txt
$SENTIEON_INSTALL_DIR/bin/sentieon plot MeanQualityByCycle -o mq-report.pdf mq_metrics.txt
$SENTIEON_INSTALL_DIR/bin/sentieon plot InsertSizeMetricAlgo -o is-report.pdf is_metrics.txt

# ******************************************
# 3. Remove Duplicate Reads. It is possible
# to remove instead of mark duplicates
# by adding the --rmdup option in Dedup
# ******************************************
$SENTIEON_INSTALL_DIR/bin/sentieon driver -t $NT -i sorted.bam --algo LocusCollector \
    --fun score_info score.txt || { echo "LocusCollector failed"; exit 1; }
$SENTIEON_INSTALL_DIR/bin/sentieon driver -t $NT -i sorted.bam --algo Dedup \
    --score_info score.txt --metrics dedup_metrics.txt deduped.bam || \
    { echo "Dedup failed"; exit 1; }
# ******************************************
# 5. Base recalibration
# ******************************************
$SENTIEON_INSTALL_DIR/bin/sentieon driver -r $FASTA -t $NT -i deduped.bam --algo QualCal \
    -k $KNOWN_DBSNP recal_data.table
$SENTIEON_INSTALL_DIR/bin/sentieon driver -r $FASTA -t $NT -i deduped.bam -q recal_data.table \
    --algo QualCal -k $KNOWN_DBSNP  recal_data.table.post
$SENTIEON_INSTALL_DIR/bin/sentieon driver -t $NT --algo QualCal --plot \
    --before recal_data.table --after recal_data.table.post recal.csv
$SENTIEON_INSTALL_DIR/bin/sentieon plot QualCal -o recal_plots.pdf recal.csv
# ******************************************
# 6b. HC Variant caller
# ******************************************
$SENTIEON_INSTALL_DIR/bin/sentieon driver -r $FASTA -t $NT -i deduped.bam -q recal_data.table \
    --algo Haplotyper -d $KNOWN_DBSNP ${SM}-hc.vcf.gz || \
    { echo "Haplotyper failed"; exit 1; }
# ******************************************
# 5b. ReadWriter to output recalibrated bam
# This stage is optional as variant callers
# can perform the recalibration on the fly
# using the before recalibration bam plus
# the recalibration table
# ******************************************
# $SENTIEON_INSTALL_DIR/bin/sentieon driver -r $FASTA -t $NT -i deduped.bam -q recal_data.table  --algo ReadWriter recaled.bam || { echo "ReadWriter failed"; exit 1; }
for i in `ls *_1_*|tr "_" "\t"|cut -f1`;do sed "s/test01/$i/g"  test01.sh >SNP_CALL$i.sh;done
for i in `ls *_1_*|tr "_" "\t"|cut -f1`;do  chmod a+x SNP_CALL$i.sh;done
  • 大群体joincalling

cd /data1/data_public/WGS/pig/XinJiang382/2.5G/rawdata/cleandata
mkdir DNAseq
vi  test01.sh
#!/bin/sh
# Copyright (c) 2016-2020 Sentieon Inc. All rights reserved
# *******************************************
# Script to perform DNA seq variant calling
# using a single sample with fastq files
# named 1.fastq.gz and 2.fastq.gz
# *******************************************
cd  /data1/data_public/WGS/pig/XinJiang382/2.5G/rawdata/cleandata
set -eu
# Update with the fullpath location of your sample fastq
SM="test01" #sample name
RGID="rg_$SM" #read group ID
PL="ILLUMINA" #or other sequencing platform
FASTQ_FOLDER="/data1/data_public/WGS/pig/XinJiang382/2.5G/rawdata/cleandata"
FASTQ_1="$FASTQ_FOLDER/${SM}_1_val_1.fq.gz"
FASTQ_2="$FASTQ_FOLDER/${SM}_2_val_2.fq.gz"  #If using Illumina paired data
# Update with the location of the reference data files
FASTA_DIR="/home/jychu/home/refs/pig/Sus_scrofa.Sscrofa11.1-106"
FASTA="$FASTA_DIR/Sus_scrofa.Sscrofa11.1.dna.toplevel.fa"
KNOWN_DBSNP="$FASTA_DIR/sus_scrofa.vcf.gz"
# Update with the location of the Sentieon software package and license file
SENTIEON_INSTALL_DIR=/home/yzhou/soft/sentieon-genomics-202112.07
export SENTIEON_LICENSE=192.168.102.80:8990  #or using licsrvr: c1n11.sentieon.com:5443
# Other settings
NT=1  #number of threads to use in computation, set to number of cores in the server
START_DIR="$FASTQ_FOLDER/DNAseq" #Determine where the output files will be stored
# You do not need to modify any of the lines below unless you want to tweak the pipeline
# ************************************************************************************************************************************************************************
# ******************************************
# 0. Setup
# ******************************************
WORKDIR="$START_DIR/${SM}" 
mkdir -p $WORKDIR
LOGFILE=$WORKDIR/run.log
exec >$LOGFILE 2>&1
cd $WORKDIR
# ******************************************
# 1. Mapping reads with BWA-MEM, sorting
# ******************************************
#The results of this call are dependent on the number of threads used. To have number of threads independent results, add chunk size option -K 10000000 
( $SENTIEON_INSTALL_DIR/bin/sentieon bwa mem -R "@RG\tID:$RGID\tSM:$SM\tPL:$PL" \
    -t $NT -K 10000000 $FASTA $FASTQ_1 $FASTQ_2 || { echo -n 'BWA error'; exit 1; } ) | \
    $SENTIEON_INSTALL_DIR/bin/sentieon util sort -r $FASTA -o sorted.bam -t $NT \
    --sam2bam -i - || { echo "Alignment failed"; exit 1; }
# ******************************************
# 2. Metrics
# ******************************************
$SENTIEON_INSTALL_DIR/bin/sentieon driver -r $FASTA -t $NT -i sorted.bam \
    --algo MeanQualityByCycle mq_metrics.txt --algo QualDistribution qd_metrics.txt \
    --algo GCBias --summary gc_summary.txt gc_metrics.txt --algo AlignmentStat \
    --adapter_seq '' aln_metrics.txt --algo InsertSizeMetricAlgo is_metrics.txt || \
    { echo "Metrics failed"; exit 1; }

$SENTIEON_INSTALL_DIR/bin/sentieon plot GCBias -o gc-report.pdf gc_metrics.txt
$SENTIEON_INSTALL_DIR/bin/sentieon plot QualDistribution -o qd-report.pdf qd_metrics.txt
$SENTIEON_INSTALL_DIR/bin/sentieon plot MeanQualityByCycle -o mq-report.pdf mq_metrics.txt
$SENTIEON_INSTALL_DIR/bin/sentieon plot InsertSizeMetricAlgo -o is-report.pdf is_metrics.txt

# ******************************************
# 3. Remove Duplicate Reads. It is possible
# to remove instead of mark duplicates
# by adding the --rmdup option in Dedup
# ******************************************
$SENTIEON_INSTALL_DIR/bin/sentieon driver -t $NT -i sorted.bam --algo LocusCollector \
    --fun score_info score.txt || { echo "LocusCollector failed"; exit 1; }
$SENTIEON_INSTALL_DIR/bin/sentieon driver -t $NT -i sorted.bam --algo Dedup \
    --score_info score.txt --metrics dedup_metrics.txt deduped.bam || \
    { echo "Dedup failed"; exit 1; }
# ******************************************
# 5. Base recalibration
# ******************************************
$SENTIEON_INSTALL_DIR/bin/sentieon driver -r $FASTA -t $NT -i deduped.bam --algo QualCal \
    -k $KNOWN_DBSNP recal_data.table
$SENTIEON_INSTALL_DIR/bin/sentieon driver -r $FASTA -t $NT -i deduped.bam -q recal_data.table \
    --algo QualCal -k $KNOWN_DBSNP  recal_data.table.post
$SENTIEON_INSTALL_DIR/bin/sentieon driver -t $NT --algo QualCal --plot \
    --before recal_data.table --after recal_data.table.post recal.csv
$SENTIEON_INSTALL_DIR/bin/sentieon plot QualCal -o recal_plots.pdf recal.csv
 # ******************************************
 # 6b. HC Variant caller for GVCF
 # ******************************************
 $SENTIEON_INSTALL_DIR/bin/sentieon driver -r $FASTA -t $NT -i deduped.bam -q recal_data.table \
      --algo Haplotyper -d $KNOWN_DBSNP --emit_conf=30 --call_conf=30 --emit_mode gvcf \
      ${SM}-hc.g.vcf.gz || { echo "Haplotyper failed"; exit 1; }

for i in `ls *_1_*|tr "_" "\t"|cut -f1`;do sed "s/test01/$i/g"  test01.sh >SNP_CALL$i.sh;done
for i in `ls *_1_*|tr "_" "\t"|cut -f1`;do  chmod a+x SNP_CALL$i.sh;done
  • 小群体joincalling

vi joincalling.sh
 #!/bin/sh
# Copyright (c) 2016-2020 Sentieon Inc. All rights reserved
# *******************************************
# Script to perform joint calling on 3 samples
# with fastq files named sample<i>_1.fastq.gz
# and sample<i>_2.fastq.gz
# *******************************************
set -eu
# Update with the fullpath location of your sample fastq
SM_LIST=$(ls *_1_val_1.fq.gz|tr "_" "\t"|cut -f1) # list of sample names
PL="ILLUMINA" #or other sequencing platform
FASTQ_FOLDER="/data1/data_public/WGS/pig/XinJiang382/2.5G/rawdata/cleandata"
FASTQ_1_SUFFIX="1_val_1.fq.gz"
FASTQ_2_SUFFIX="2_val_2.fq.gz"
# Update with the location of the reference data files
FASTA_DIR="/home/jychu/home/refs/pig/Sus_scrofa.Sscrofa11.1-106"
FASTA="$FASTA_DIR/Sus_scrofa.Sscrofa11.1.dna.toplevel.fa"
KNOWN_DBSNP="$FASTA_DIR/sus_scrofa.vcf.gz"
# Update with the location of the Sentieon software package and license file
SENTIEON_INSTALL_DIR=/home/yzhou/soft/sentieon-genomics-202112.07
export SENTIEON_LICENSE=192.168.102.80:8990  #or using licsrvr: c1n11.sentieon.com:5443
# Other settings
NT=40  #number of threads to use in computation, set to number of cores in the server
START_DIR="/data1/data_public/WGS/pig/XinJiang382/2.5G/rawdata/cleandata/DNAseq_joint" #Determine where the output files will be stored
# You do not need to modify any of the lines below unless you want to tweak the pipeline
# ************************************************************************************************************************************************************************
# ******************************************
# 0. Setup
# ******************************************
WORKDIR="$START_DIR/joint_call"
mkdir -p $WORKDIR
LOGFILE=$WORKDIR/run.log
exec >$LOGFILE 2>&1
# ******************************************
# 0. Process all samples independently
# ******************************************
GVCF_INPUTS=""
for SM in $SM_LIST; do
  RGID="rg_$SM"
  mkdir $WORKDIR/$SM
  cd $WORKDIR/$SM
  # ******************************************
  # 1. Mapping reads with BWA-MEM, sorting
  # ******************************************
  ( $SENTIEON_INSTALL_DIR/bin/sentieon bwa mem -R "@RG\tID:$RGID\tSM:$SM\tPL:$PL" \
      -t $NT -K 10000000 $FASTA $FASTQ_FOLDER/${SM}_$FASTQ_1_SUFFIX \
      $FASTQ_FOLDER/${SM}_$FASTQ_2_SUFFIX || { echo -n 'bwa error'; exit 1; } ) | \
      $SENTIEON_INSTALL_DIR/bin/sentieon util sort -r $FASTA -o sorted.bam -t $NT --sam2bam -i - || \
      { echo "Alignment failed"; exit 1; }
  # ******************************************
  # 2. Metrics
  # ******************************************
  $SENTIEON_INSTALL_DIR/bin/sentieon driver -r $FASTA -t $NT -i sorted.bam \
      --algo MeanQualityByCycle mq_metrics.txt --algo QualDistribution qd_metrics.txt \
      --algo GCBias --summary gc_summary.txt gc_metrics.txt --algo AlignmentStat \
      --adapter_seq '' aln_metrics.txt --algo InsertSizeMetricAlgo is_metrics.txt || \
      { echo "Metrics failed"; exit 1; }
  $SENTIEON_INSTALL_DIR/bin/sentieon plot GCBias -o gc-report.pdf gc_metrics.txt
  $SENTIEON_INSTALL_DIR/bin/sentieon plot QualDistribution -o qd-report.pdf qd_metrics.txt
  $SENTIEON_INSTALL_DIR/bin/sentieon plot MeanQualityByCycle -o mq-report.pdf mq_metrics.txt
  $SENTIEON_INSTALL_DIR/bin/sentieon plot InsertSizeMetricAlgo -o is-report.pdf is_metrics.txt
 
  # ******************************************
  # 3. Remove Duplicate Reads. It is possible
  # to remove instead of mark duplicates
  # by adding the --rmdup option in Dedup
  # ******************************************
  $SENTIEON_INSTALL_DIR/bin/sentieon driver -t $NT -i sorted.bam --algo LocusCollector \
      --fun score_info score.txt || { echo "LocusCollector failed"; exit 1; }

  $SENTIEON_INSTALL_DIR/bin/sentieon driver -t $NT -i sorted.bam --algo Dedup \
      --score_info score.txt --metrics dedup_metrics.txt deduped.bam || \
      { echo "Dedup failed"; exit 1; }
  # ******************************************
  # 2a. Coverage metrics
  # ******************************************
  $SENTIEON_INSTALL_DIR/bin/sentieon driver -r $FASTA -t $NT -i deduped.bam \
      --algo CoverageMetrics coverage_metrics || { echo "CoverageMetrics failed"; exit 1; }
  # ******************************************
  # 5. Base recalibration
  # ******************************************
  $SENTIEON_INSTALL_DIR/bin/sentieon driver -r $FASTA -t $NT -i deduped.bam --algo QualCal \
      -k $KNOWN_DBSNP  recal_data.table
  $SENTIEON_INSTALL_DIR/bin/sentieon driver -r $FASTA -t $NT -i deduped.bam -q recal_data.table \
      --algo QualCal -k $KNOWN_DBSNP  recal_data.table.post
  $SENTIEON_INSTALL_DIR/bin/sentieon driver -t $NT --algo QualCal --plot \
      --before recal_data.table --after recal_data.table.post recal.csv
  $SENTIEON_INSTALL_DIR/bin/sentieon plot QualCal -o recal_plots.pdf recal.csv
  # ******************************************
  # 6b. HC Variant caller for GVCF
  # ******************************************
  $SENTIEON_INSTALL_DIR/bin/sentieon driver -r $FASTA -t $NT -i deduped.bam -q recal_data.table \
      --algo Haplotyper -d $KNOWN_DBSNP --emit_conf=30 --call_conf=30 --emit_mode gvcf \
      output-hc.g.vcf.gz || { echo "Haplotyper failed"; exit 1; }
  GVCF_INPUTS="$GVCF_INPUTS -v $WORKDIR/$SM/output-hc.g.vcf.gz"
done
# ******************************************
# Perform the joint calling 
# ******************************************
$SENTIEON_INSTALL_DIR/bin/sentieon driver -r $FASTA -t $NT --algo GVCFtyper $GVCF_INPUTS \
    -d $KNOWN_DBSNP $WORKDIR/output-jc.vcf.gz || { echo "GVCFtyper failed"; exit 1; }
chmod a+x joincalling.sh
nohup ./joincalling.sh &
©著作权归作者所有,转载或内容合作请联系作者
平台声明:文章内容(如有图片或视频亦包括在内)由作者上传并发布,文章内容仅代表作者本人观点,简书系信息发布平台,仅提供信息存储服务。

推荐阅读更多精彩内容