#!/bin/bash # Notice - this is a custom Perl script as described in the 2016 Nature Genetics paper # "A high-quality carrot genome assembly provides new insights into carotenoid # accumulation and Asterid genome evolution" # This program was referenced by this line: # ... "Top 1% of FST values were determined and visualized by a custom Perl script." ... # Please be aware that this is not a general purpose program, it is highly specific to # our analysis, and is provided for informative purposes only. Note that the helper # program 34.fstplots.bins.pl is called by this program # Plot Megan's fst values on chromosomes with Jan 29 data pfx=`echo $0 | sed -e "s/^.*\///" -e "s/\..*$//"` binsize=100000 ymax=271 if [ "$1" != "" ]; then binsize="$1"; fi if [ "$2" != "" ]; then ymax="$2"; fi fstdir="/gen/0065/scripts/megan/20150129_sorted_redo_top_1_fst/" outdir="$pfx.outputdir.bis" report="/gen/0065/scripts/notebook/$pfx.combinedreport.bis.htm" htmlheader=' VCRU Bioinformatics - Password Protected Pages - Lab Notebook

This page was last updated on

' htmlfooter=' ' ##### Initialization if [ ! -d "notebook/$outdir" ]; then echo "Creating output directory \"notebook/$outdir\"" mkdir "notebook/$outdir" r=$?; if [ $r != 0 ]; then echo "Error code $r line $LINENO script $0, halting"; exit $r; fi fi ##### Main loop echo "$htmlheader" > $report for i in 1 2 3 4 5 6 7 ; do for j in 2 3 4 5 6 7 ; do if [ $i -lt $j ] ; then echo "Plotting $chr $i vs. $j $(date)" fstfile="$fstdir/sorted_redo_top_1_fst_group${i}_group${j}.weir.fst" if [ ! -s "$fstfile" ]; then echo "Skipping missing i=$i j=$j file \"$fstfile\"" else plot="$outdir/$pfx.all.$i.$j.png" tn="$outdir/$pfx.all.$i.$j.tn.png" tmpdata="$pfx.$i.$j.bis.DELETE.tsv" echo "Converting coordinates to combine chromosomes" ./34.fstplots.bins.pl "$fstfile" "$tmpdata" $binsize "$i vs. $j" r=$?; if [ $r != 0 ]; then echo "Error code $r line $LINENO script $0, halting"; exit $r; fi bb.plot \ --infile="$tmpdata" \ --outfile="notebook/$plot" \ --xaxis=nt \ --yaxis=Number of variants \ --xsize=4096 \ --ymin=0 --ymax=$ymax \ --type=F \ --linewidth=1 \ --xcolumn=1 \ --ycolumn=2 --ylabel="Count" \ --ycolumn=3 --ylabel="Sum" \ --thumbnail \ --xmax=362150989 \ --extra="set xtics 0,10000000" --extra="set mxtics 10" \ --extra="set arrow from 51478430,0 to 51478430,$ymax nohead lw 1 lc rgb 'black'" \ --extra="set arrow from 95400503,0 to 95400503,$ymax nohead lw 1 lc rgb 'black'" \ --extra="set arrow from 145730752,0 to 145730752,$ymax nohead lw 1 lc rgb 'black'" \ --extra="set arrow from 181662705,0 to 181662705,$ymax nohead lw 1 lc rgb 'black'" \ --extra="set arrow from 223632787,0 to 223632787,$ymax nohead lw 1 lc rgb 'black'" \ --extra="set arrow from 260252785,0 to 260252785,$ymax nohead lw 1 lc rgb 'black'" \ --extra="set arrow from 296640290,0 to 296640290,$ymax nohead lw 1 lc rgb 'black'" \ --extra="set arrow from 328396095,0 to 328396095,$ymax nohead lw 1 lc rgb 'black'" r=$?; if [ $r != 0 ]; then echo "Error code $r line $LINENO script $0, halting"; exit $r; fi rm "$tmpdata" echo "
  • All chromosomes $i vs. $j —" >> "$report" echo "
    " >> "$report" echo "
    • " >> "$report" fi # if not missing input file fi # if i -lt j done # for j done # for i echo "" >> "$report" echo "$htmlfooter" >> $report cp -puv $0 notebook/ echo "Done" exit 0 #eof