Skip to content
Snippets Groups Projects
Commit 1588fae4 authored by Brandi Cantarel's avatar Brandi Cantarel
Browse files

update union and snsp with mnv->snv

parent 4b2bd485
No related merge requests found
...@@ -33,21 +33,21 @@ if [[ -z $SLURM_CPUS_ON_NODE ]] ...@@ -33,21 +33,21 @@ if [[ -z $SLURM_CPUS_ON_NODE ]]
then then
SLURM_CPUS_ON_NODE=1 SLURM_CPUS_ON_NODE=1
fi fi
if [[ -a "${index_path}/dbSnp.vcf.gz" ]] if [[ -s "${index_path}/dbSnp.vcf.gz" ]]
then then
dbsnp="${index_path}/dbSnp.vcf.gz" dbsnp="${index_path}/dbSnp.vcf.gz"
else else
echo "Missing dbSNP File: ${index_path}/dbSnp.vcf.gz" echo "Missing dbSNP File: ${index_path}/dbSnp.vcf.gz"
usage usage
fi fi
if [[ -a "${index_path}/GoldIndels.vcf.gz" ]] if [[ -s "${index_path}/GoldIndels.vcf.gz" ]]
then then
knownindel="${index_path}/GoldIndels.vcf.gz" knownindel="${index_path}/GoldIndels.vcf.gz"
else else
echo "Missing InDel File: ${index_path}/GoldIndels.vcf.gz" echo "Missing InDel File: ${index_path}/GoldIndels.vcf.gz"
usage usage
fi fi
if [[ -a "${index_path}/genome.fa" ]] if [[ -s "${index_path}/genome.fa" ]]
then then
reffa="${index_path}/genome.fa" reffa="${index_path}/genome.fa"
dict="${index_path}/genome.dict" dict="${index_path}/genome.dict"
...@@ -87,7 +87,10 @@ then ...@@ -87,7 +87,10 @@ then
gatk4_dbsnp=${index_path}/clinseq_prj/dbSnp.gatk4.vcf.gz gatk4_dbsnp=${index_path}/clinseq_prj/dbSnp.gatk4.vcf.gz
user=$USER user=$USER
module load gatk/4.x singularity/2.6.1 module load gatk/4.x singularity/2.6.1
mkdir /tmp/${user} if [[ ! -e "/tmp/${user}" ]]
then
mkdir /tmp/${user}
fi
export TMP_HOME=/tmp/${user} export TMP_HOME=/tmp/${user}
gvcflist='' gvcflist=''
for i in *.bam; do for i in *.bam; do
...@@ -98,7 +101,9 @@ then ...@@ -98,7 +101,9 @@ then
gvcflist="$gvcflist --variant ${prefix}.gatk.g.vcf" gvcflist="$gvcflist --variant ${prefix}.gatk.g.vcf"
done done
singularity exec -H /tmp/$user /project/apps/singularity-images/gatk4/gatk-4.x.simg /gatk/gatk --java-options "-Xmx32g" GenotypeGVCFs $gvcflist -R ${reffa} -D ${gatk4_dbsnp} -O gatk.vcf interval=`cat ${reffa}.fai |cut -f 1 |grep -v decoy |grep -v 'HLA' |grep -v alt |grep -v 'chrUn' |grep -v 'random' |paste -sd "," -`
singularity exec -H /tmp/$user /project/apps/singularity-images/gatk4/gatk-4.x.simg /gatk/gatk --java-options "-Xmx32g" GenomicsDBImport $gvcflist --genomicsdb-workspace-path gendb --intervals $interval
singularity exec -H /tmp/$user /project/apps/singularity-images/gatk4/gatk-4.x.simg /gatk/gatk --java-options "-Xmx32g" GenotypeGVCFs -V gendb://gendb -R ${reffa} -D ${gatk4_dbsnp} -O gatk.vcf
bcftools norm -c s -f ${reffa} -w 10 -O v gatk.vcf | vcf-annotate -n --fill-type gatk.vcf | bgzip > ${pair_id}.gatk.vcf.gz bcftools norm -c s -f ${reffa} -w 10 -O v gatk.vcf | vcf-annotate -n --fill-type gatk.vcf | bgzip > ${pair_id}.gatk.vcf.gz
tabix ${pair_id}.gatk.vcf.gz tabix ${pair_id}.gatk.vcf.gz
elif [[ $algo == 'platypus' ]] elif [[ $algo == 'platypus' ]]
......
...@@ -30,4 +30,5 @@ perl $baseDir\/uniform_vcf_gt.pl $pair_id $vcf ...@@ -30,4 +30,5 @@ perl $baseDir\/uniform_vcf_gt.pl $pair_id $vcf
bgzip -f ${pair_id}.uniform.vcf bgzip -f ${pair_id}.uniform.vcf
j=${pair_id}.uniform.vcf.gz j=${pair_id}.uniform.vcf.gz
tabix -f $j tabix -f $j
bcftools norm -m - -O z -o ${pair_id}.norm.vcf.gz $j bcftools norm -m - -O v $j | /project/shared/bicf_workflow_ref/seqprg/vt/vt decompose_blocksub - -a -o ${pair_id}.norm.vcf
bgzip -f ${pair_id}.norm.vcf
...@@ -125,10 +125,10 @@ foreach $vcf (@vcffiles) { ...@@ -125,10 +125,10 @@ foreach $vcf (@vcffiles) {
my @filts = split(";",$filter); my @filts = split(";",$filter);
my %filterqc = map {$_ => 1} @filts; my %filterqc = map {$_ => 1} @filts;
delete $filterqc{'.'}; delete $filterqc{'.'};
if ($gtfilt{'StrandBias'} || ($hash{FS} && $hash{FS} > 60) || if ($gtfilt{'StrandBias'} || ($hash{FS} && $hash{FS} > 60)) { # ||($hash{SAP} && $hash{SAP} > 20)) {
($hash{SAP} && $hash{SAP} > 20)) {
$filterqc{strandBias} = 1; $filterqc{strandBias} = 1;
}if (scalar(keys %filterqc) > 1) { $annot .= ';strandBias=1';
}if (scalar(keys %filterqc) > 0) {
$filter = join(";",keys %filterqc); $filter = join(";",keys %filterqc);
}else { }else {
$filter = '.'; $filter = '.';
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment