From 67d6780120902ec6bfdcbb0913ac1c910d24af61 Mon Sep 17 00:00:00 2001 From: Isaac Turner Date: Wed, 20 Jan 2016 14:58:57 +0000 Subject: [PATCH] Update klebsiella pneumoniae analysis [thesis][skip CI] --- .../kleb_pneumoniae/assembly/get-max-covg.sh | 10 +++ .../assembly/make-dot-plots.sh | 80 +++++++++++++++++++ .../kleb_pneumoniae/cortex/analysis.sh | 2 +- .../kleb_pneumoniae/freebayes/analysis.sh | 2 +- .../kleb_pneumoniae/indels/kleb-plot-hist.R | 5 +- .../kleb_pneumoniae/mcrun/analysis.sh | 2 +- .../kleb_pneumoniae/platypus/analysis.sh | 2 +- 7 files changed, 97 insertions(+), 6 deletions(-) create mode 100755 results/klebsiella/kleb_pneumoniae/assembly/get-max-covg.sh create mode 100755 results/klebsiella/kleb_pneumoniae/assembly/make-dot-plots.sh diff --git a/results/klebsiella/kleb_pneumoniae/assembly/get-max-covg.sh b/results/klebsiella/kleb_pneumoniae/assembly/get-max-covg.sh new file mode 100755 index 00000000..eabafcb8 --- /dev/null +++ b/results/klebsiella/kleb_pneumoniae/assembly/get-max-covg.sh @@ -0,0 +1,10 @@ +#!/bin/bash + +# usage: get-max-covg.sh +# output: + +genome=$1 +seqfile=$2 + +dnacat -L "$seqfile" | cut -f2 | sort -rn | \ + awk '{if(x+$1>'$genome'){exit;} x+=$1; n+=1; l=$1; if(!f){f=$1}} END{print n,x,f,l}' diff --git a/results/klebsiella/kleb_pneumoniae/assembly/make-dot-plots.sh b/results/klebsiella/kleb_pneumoniae/assembly/make-dot-plots.sh new file mode 100755 index 00000000..80ee70cf --- /dev/null +++ b/results/klebsiella/kleb_pneumoniae/assembly/make-dot-plots.sh @@ -0,0 +1,80 @@ +#!/bin/bash + +set -ou pipefail +set -o xtrace + +# +# Make MUMMER dot plot +# Requires: mummer (nucmer, mummerplot), gnuplot, dnacat, ghostscript (ps2pdf) +# +# Clean up: +# rm -rf links plain truthref + +LINKCONTIGS=k31_links/klebPneu.rmdup.fa +PLAINCONTIGS=k31_plain/klebPneu.rmdup.fa +TRUTH=../truth/CAV1016.fa +REF=../ref/GCF_000016305.1_ASM1630v1_genomic.fa + +CONTIGSTATS=~/mccortex/libs/bioinf-perl/fastn_scripts/contig_stats.pl +ASSEM_NCONTIGS=./get-max-covg.sh + +REFLEN=$(dnacat -P $REF | awk '{x+=length($1)} END{print x}') + +# Fix issues, re-run gnuplot +function fix_mummer_gp { + local prefix=$1 + sed -i '.tmp' 's/^set mouse clipboardformat/\#set mouse clipboardformat/g' $prefix.gp + gnuplot $prefix.gp +} + +# make plot $prefix'_cov.pdf' +# make plot $prefix'_dot.pdf' +function mk_mummer_plots { + local prefix=$1; shift + # local otheropts=$2 + # coverage plot + mummerplot -t postscript -c -p $prefix"_cov" $prefix.delta + fix_mummer_gp $prefix"_cov" + ps2pdf $prefix"_cov.ps" $prefix"_cov.pdf" + # dot plot + # removed: --breaklen 10000 + mummerplot -t postscript -l "$@" -p $prefix"_dot" $prefix.delta + fix_mummer_gp $prefix"_dot" + ps2pdf $prefix"_dot.ps" $prefix"_dot.pdf" +} + +function get_longest_contigs { + local seqin=$1 + local seqout=$2 + local reflen=$3 + local ncontigs=$($ASSEM_NCONTIGS $reflen $seqin | awk '{print $1}') + set +o xtrace + dnacat -P $seqin | awk '{print length($1),$1}' | sort -k1rn,2 | \ + awk '{print $2}' | head -$ncontigs | \ + dnacat -F -M <(i=0; while true; do echo "r$i"; ((i++)); done) > $seqout + set -o xtrace +} + +# mkdir -p links +# get_longest_contigs $LINKCONTIGS links/contigs.fa $REFLEN +# nucmer -mum -p links/links $TRUTH links/contigs.fa +# mk_mummer_plots links/links +# nucmer -mum -p links/linksswap links/contigs.fa $TRUTH +# mk_mummer_plots links/linksswap + +# mkdir -p plain +# get_longest_contigs $PLAINCONTIGS plain/contigs.fa $REFLEN +# nucmer -mum -p plain/plain $TRUTH plain/contigs.fa +# mk_mummer_plots plain/plain +# nucmer -mum -p plain/plainswap plain/contigs.fa $TRUTH +# mk_mummer_plots plain/plainswap + +# # Make plot between ref and truth +# mkdir -p truthref +# nucmer -mum -p truthref/truthref $TRUTH $REF +# mk_mummer_plots truthref/truthref +# nucmer -mum -p truthref/truthrefswap $REF $TRUTH +# mk_mummer_plots truthref/truthrefswap + +$CONTIGSTATS links/contigs.fa > links/contig_stats.txt +$CONTIGSTATS plain/contigs.fa > plain/contig_stats.txt diff --git a/results/klebsiella/kleb_pneumoniae/cortex/analysis.sh b/results/klebsiella/kleb_pneumoniae/cortex/analysis.sh index afc1d2ec..40714952 100755 --- a/results/klebsiella/kleb_pneumoniae/cortex/analysis.sh +++ b/results/klebsiella/kleb_pneumoniae/cortex/analysis.sh @@ -6,7 +6,7 @@ set -o xtrace CTXDIR=../../../../ REF=../ref/GCF_000016305.1_ASM1630v1_genomic.fna.gz -MUMMER=../mummer/ref_qry.good.nodupes.vcf.gz +MUMMER=../mummer/mummer.vcf.gz TRUTH=../truth/CAV1016.fa MAPPING_TEST=$CTXDIR/scripts/analysis/mapping-vars-test.sh MUMMER_ISEC=$CTXDIR/scripts/analysis/mummer-vcf-isec.sh diff --git a/results/klebsiella/kleb_pneumoniae/freebayes/analysis.sh b/results/klebsiella/kleb_pneumoniae/freebayes/analysis.sh index c01172a0..7f9c391b 100755 --- a/results/klebsiella/kleb_pneumoniae/freebayes/analysis.sh +++ b/results/klebsiella/kleb_pneumoniae/freebayes/analysis.sh @@ -6,7 +6,7 @@ set -o xtrace CTXDIR=../../../../ REF=../ref/GCF_000016305.1_ASM1630v1_genomic.fna.gz -MUMMER=../mummer/ref_qry.good.nodupes.vcf.gz +MUMMER=../mummer/mummer.vcf.gz TRUTH=../truth/CAV1016.fa MAPPING_TEST=$CTXDIR/scripts/analysis/mapping-vars-test.sh MUMMER_ISEC=$CTXDIR/scripts/analysis/mummer-vcf-isec.sh diff --git a/results/klebsiella/kleb_pneumoniae/indels/kleb-plot-hist.R b/results/klebsiella/kleb_pneumoniae/indels/kleb-plot-hist.R index dc1ac42d..d75716e7 100755 --- a/results/klebsiella/kleb_pneumoniae/indels/kleb-plot-hist.R +++ b/results/klebsiella/kleb_pneumoniae/indels/kleb-plot-hist.R @@ -18,6 +18,8 @@ colours <- c('red', 'green', 'blue', 'black', 'purple') # points <- c('o','o','o', 'o', 'o') points <- rep(20, 5) +title <- expression(paste(italic('klebsiella pneumoniae'),' indels')) + ncalls <- length(hists) lim <- 100 m <- matrix(0, nrow=2*lim+1, ncol=ncalls+1) @@ -44,8 +46,7 @@ pdf(file="kleb_indels.pdf",width=6,height=6) plot(jitter(x,factor=jf), m[,2], log="y", col=colours[1],pch=points[1], - main="Klebsiella Pneumoniae Indels", - xlab="indel size", ylab="count (log)") + main=title, xlab="indel size", ylab="count (log)") for(i in 2:length(names)) { points(jitter(x,factor=jf),m[,i+1],col=colours[i],pch=points[i]) # ,type='b' diff --git a/results/klebsiella/kleb_pneumoniae/mcrun/analysis.sh b/results/klebsiella/kleb_pneumoniae/mcrun/analysis.sh index db005548..fc0fd2bd 100755 --- a/results/klebsiella/kleb_pneumoniae/mcrun/analysis.sh +++ b/results/klebsiella/kleb_pneumoniae/mcrun/analysis.sh @@ -10,7 +10,7 @@ function myreadlink() { } REF=$(myreadlink ../ref/GCF_000016305.1_ASM1630v1_genomic.fna.gz) -MUMMER=$(myreadlink ../mummer/ref_qry.good.nodupes.vcf.gz) +MUMMER=$(myreadlink ../mummer/mummer.vcf.gz) TRUTH=$(myreadlink ../truth/CAV1016.fa) BCFTOOLS=$(myreadlink $CTXDIR/libs/bcftools/bcftools) MAPPING_TEST=$(myreadlink $CTXDIR/scripts/analysis/mapping-vars-test.sh) diff --git a/results/klebsiella/kleb_pneumoniae/platypus/analysis.sh b/results/klebsiella/kleb_pneumoniae/platypus/analysis.sh index 6e2f5b15..f18b1d1e 100755 --- a/results/klebsiella/kleb_pneumoniae/platypus/analysis.sh +++ b/results/klebsiella/kleb_pneumoniae/platypus/analysis.sh @@ -6,7 +6,7 @@ set -o xtrace CTXDIR=../../../../ REF=../ref/GCF_000016305.1_ASM1630v1_genomic.fna.gz -MUMMER=../mummer/ref_qry.good.nodupes.vcf.gz +MUMMER=../mummer/mummer.vcf.gz TRUTH=../truth/CAV1016.fa MAPPING_TEST=$CTXDIR/scripts/analysis/mapping-vars-test.sh MUMMER_ISEC=$CTXDIR/scripts/analysis/mummer-vcf-isec.sh