Skip to content

Commit

Permalink
Update klebsiella pneumoniae analysis [thesis][skip CI]
Browse files Browse the repository at this point in the history
  • Loading branch information
noporpoise committed Jan 20, 2016
1 parent bd02e0b commit 67d6780
Show file tree
Hide file tree
Showing 7 changed files with 97 additions and 6 deletions.
10 changes: 10 additions & 0 deletions results/klebsiella/kleb_pneumoniae/assembly/get-max-covg.sh
Original file line number Diff line number Diff line change
@@ -0,0 +1,10 @@
#!/bin/bash

# usage: get-max-covg.sh <genome-size> <in.fa>
# output: <number_of_contigs> <sum_of_lengths> <longest_len> <shortest_len>

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}'
80 changes: 80 additions & 0 deletions results/klebsiella/kleb_pneumoniae/assembly/make-dot-plots.sh
Original file line number Diff line number Diff line change
@@ -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
2 changes: 1 addition & 1 deletion results/klebsiella/kleb_pneumoniae/cortex/analysis.sh
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
2 changes: 1 addition & 1 deletion results/klebsiella/kleb_pneumoniae/freebayes/analysis.sh
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
5 changes: 3 additions & 2 deletions results/klebsiella/kleb_pneumoniae/indels/kleb-plot-hist.R
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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'
Expand Down
2 changes: 1 addition & 1 deletion results/klebsiella/kleb_pneumoniae/mcrun/analysis.sh
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
2 changes: 1 addition & 1 deletion results/klebsiella/kleb_pneumoniae/platypus/analysis.sh
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down

0 comments on commit 67d6780

Please sign in to comment.