Sentieon | 鸡全基因组（WGS）分析流程

#!/bin/bash
echo $0 \$SAMPLEID \$WORKDIR \$FASTQ_1 \$FASTQ_2 \$FASTA \$SUFFIX \$DATATYPE \$KEEP_CLEAN \$KEEP_BAM \$PLOIDY
set -euxo pipefail

export SAMPLEID=$1
export WORKDIR=$2
export FASTQ_1=$3
export FASTQ_2=$4
export FASTA=$5
BSUFFIX=$6
TYPE=${7:-"raw"}
KEEP_CLEAN=${8:-keep}
KEEP_BAM=${9:-keep}
PLOIDY=${10:-2}
LOGFILE=$SAMPLEID.run.log
export root=/Path/sentieon/202503/sentieon-genomics-202503/bin/

if [ "$BSUFFIX" = "bam" ] || [ "$BSUFFIX" = "cram" ] || [ "$BSUFFIX" = " " ]; then
    echo "$BSUFFIX check"
else
    die "Error: check 6th blank, BSUFFIX must be 'bam' or 'cram' or space"
fi

if [ -e $LOGFILE ];then
    export PLATFORM=$(awk '/Platform/{print $NF;exit}'$LOGFILE)
else
    export count=$(zcat $FASTQ_1|head -n 1|awk '{print NF}')
    if [ $count -eq 1 ];then
        export PLATFORM="DNBSEQ"
    elif [ $count == 2 ];then
        export PLATFORM="ILLUMINA"
    elif [ $count == 3 ];then
        export PLATFORM="ILLUMINA"
    else
        echo"Unrecognized platform"
        export PLATFORM="ILLUMINA"
    fi
fi

FAI=$FASTA.fai
if [ "$BSUFFIX" = "bam" ] || [ "$BSUFFIX" = "cram" ];then
export SUFFIX=$(awk -v preset=$BSUFFIX 'BEGIN{max=0}{if($2>max)max=$2}END{if(max>536870911){print "cram"}else{print preset}}' $FAI)
else
export SUFFIX=$(awk 'BEGIN{max=0}{if($2>max)max=$2}END{if(max>536870911){print "cram"}else{print "bam"}}' $FAI)
fi

export TMPDIR=$WORKDIR
export THREADS=$(nproc)
[ -e $WORKDIR ]||mkdir -p $WORKDIR
cd $WORKDIR
exec >>$LOGFILE 2>&1
echo "SampleID:" $SAMPLEID
echo "DataType:" $TYPE
……

timer(){
    start_time=$(date +%s)
    eval $2 && touch $3
    end_time=$(date +%s)
    cost_time=$[ $end_time-$start_time ]
    echo -e "TIMER: $1\t$(($cost_time/60)) min $(($cost_time%60)) s"
}

raw2clean(){
    cmd="fastp -w 16 -i $FASTQ_1 -I $FASTQ_2 -o $clean1 -O $clean2 -j $SAMPLEID.qc.json -h $SAMPLEID.qc.html&&rm $FASTQ_1 $FASTQ_2"
    timer raw2clean "$cmd" qc.ok
}

raw2clean(){
    cmd="fastp -w 16 -i $FASTQ_1 -o $clean1 -j $SAMPLEID.qc.json -h $SAMPLEID.qc.html&&rm $FASTQ_1"
    timer raw2clean "$cmd" qc.ok
}

alignment(){
    tag="@RG\tID:rg_$SAMPLEID\tSM:$SAMPLEID\tPL:$PLATFORM"
    cmd="sentieon bwa mem -R \"$tag\" -t $THREADS -K 10000000 -x $ML_MODEL/bwa.model $FASTA $clean1 $clean2|sentieon util sort --temp_dir $TMPDIR -r $FASTA -o $SAMPLEID.sorted.$SUFFIX -t $THREADS --sam2bam -i -"
    echo $cmd
    timer alignment "$cmd" align.ok
}

alignment(){
    tag="@RG\tID:rg_$SAMPLEID\tSM:$SAMPLEID\tPL:$PLATFORM"
    cmd="$root/sentieon bwa mem -R \"$tag\" -t $THREADS -K 10000000 $FASTA -x $ML_MODEL/bwa.model $clean1 |$root/sentieon util sort --temp_dir $TMPDIR -r $FASTA -o $SAMPLEID.sorted.$SUFFIX -t $THREADS --sam2bam -i -"
    echo $cmd
    timer alignment "$cmd" align.ok
}

metrics(){
    cmd="sentieon driver --temp_dir $TMPDIR -r $FASTA -t $THREADS -i $SAMPLEID.sorted.$SUFFIX --algo WgsMetricsAlgo $SAMPLEID.WGS_METRICS.txt --algo MeanQualityByCycle $SAMPLEID.mq_metrics.txt --algo QualDistribution $SAMPLEID.qd_metrics.txt --algo GCBias --summary $SAMPLEID.gc_summary.txt $SAMPLEID.gc_metrics.txt --algo AlignmentStat $SAMPLEID.aln_metrics.txt --algo BaseDistributionByCycle $SAMPLEID.bd_metrics.txt --algo QualityYield $SAMPLEID.qy_metrics.txt --algo InsertSizeMetricAlgo $SAMPLEID.is_metrics.txt"
    timer metrics "$cmd" metrics.ok
}

dedup(){
    cmd="sentieon driver -r $FASTA --temp_dir $TMPDIR -t $THREADS -i $SAMPLEID.sorted.$SUFFIX --algo LocusCollector --fun score_info $SAMPLEID.score.txt&&sentieon driver -r $FASTA -t $THREADS -i $SAMPLEID.sorted.$SUFFIX --algo Dedup --score_info $SAMPLEID.score.txt --metrics $SAMPLEID.dedup_metrics.txt $SAMPLEID.deduped.$SUFFIX&&rm $SAMPLEID.sorted.$SUFFIX* $SAMPLEID.score.txt*"
    timer markdup "$cmd" markdup.ok
}

bqsr(){
    sentieon driver -t $THREADS -r $FASTA -i $SAMPLEID.deduped.$SUFFIX --algo QualCal -k $KNOWN_SITES $SAMPLEID.RECAL_DATA.TABLE
    sentieon driver -t $THREADS -r $FASTA -i $SAMPLEID.deduped.$SUFFIX -q $SAMPLEID.RECAL_DATA.TABLE --algo QualCal -k $KNOWN_SITES $SAMPLEID.RECAL_DATA.TABLE.POST --algo ReadWriter $SAMPLEID.deduped.RECALIBRATED.$SUFFIX
    timer markdup "$cmd" markdup.ok
}

dnaseq(){
    cmd="sentieon driver --temp_dir $TMPDIR -r $FASTA -t $THREADS -i $SAMPLEID.deduped.$SUFFIX --algo Haplotyper --emit_conf=30 --call_conf=30 --emit_mode gvcf --ploidy $PLOIDY $SAMPLEID.hc.gvcf.gz && md5sum $SAMPLEID.hc.gvcf.gz > $SAMPLEID.hc.gvcf.gz.md5"
    timer Haplotyper "$cmd" hc.ok
}

#!/bin/bash
[ $# -eq 0 ]&&echo Usage: $(basename $0) \$FASTA \$GVCF_LIST  \$NUM \$DATADIR&&exit
start_time=$(date +%s)

set -euo pipefail

DATADIR=$4
FASTA=$1 #"$FASTA_DIR/genomeEN_split.fa"
GVCF_LIST=$2
NUM=$3

NT=$(nproc) #number of threads to use in computation, set to number of cores in the server
WORKDIR="$DATADIR/JointCall-${NUM}"
[ -e $WORKDIR ]||mkdir -p $WORKDIR
#[ -e $file ]&&exit 0
cd $WORKDIR
LOGFILE=$WORKDIR/joint-call${NUM}_run.log
exec >$LOGFILE 2>&1

root=/APP/u22/x86_com/sentieon/202503/sentieon-genomics-202503
cat $GVCF_LIST|$root/bin/sentieon driver -r $FASTA  --algo GVCFtyper   \
      $WORKDIR/output${NUM}-jc.vcf.gz - || { echo echo "GVCFtyper failed"; exit 1; }

end_time=$(date +%s)
cost_time=$[ $end_time-$start_time ]
echo "joint calling time is $(($cost_time/60))min $(($cost_time%60))s"
echo "the Joint-calling done at `date +%H:%M:%S` !!!"

wget -c 
https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/016/699/485/GCF_016699485.2_bGalGal1.mat.broiler.GRCg7b/GCF_016699485.2_bGalGal1.mat.broiler.GRCg7b_genomic.fna.gz

curl -C - -O --progress-bar 
https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/016/699/485/GCF_016699485.2_bGalGal1.mat.broiler.GRCg7b/GCF_016699485.2_bGalGal1.mat.broiler.GRCg7b_genomic.fna.gz

gunzip GCF_016699485.2_bGalGal1.mat.broiler.GRCg7b_genomic.fna.gz