module load sratoolkit bwa samtools bedtools R python/2.7.9
set -x

# Download and index reference
REF=dmel-all-chromosome-r5.12.fasta
echo "Downloading and indexing reference"
[ -e $REF ] || ( wget -q ftp://ftp.flybase.net/genomes/Drosophila_melanogaster/dmel_r5.12_FB2008_09/fasta/dmel-all-chromosome-r5.12.fasta.gz && gzip -d $REF.gz )
echo "Making samtools index"
[ -e $REF.fai ] || samtools faidx $REF
echo "Making bwa index"
[ -e $REF.bwt ] || bwa index $REF

# Download and rename data
grep "KC167" lubelsky2014.csv | cut -f 5,8 | sed 's/GSM[^:][^:]*: //' | sed 's/Replicate /R/' | sed 's/ ex[^	][^	]*//' | sed "s/ WT//" | sed "s/ /_/g" | awk '{print "[ -e "$1".fastq ] || (prefetch "$2" && fastq-dump --split-3 "$2" && mv "$2".fastq "$1".fastq)"}' | xargs -L 1 -I {} -P 4 bash -c {}

# Combine replicates
ls *_R[12].fastq | sed "s/_R[12].fastq//" | sort -u | awk '{print "[ -e "$1".fastq ] || cat "$1"_R[12].fastq > "$1".fastq"}' | xargs -L 1 -I {} -P 4 bash -c {}

# Align reads and keep unique alignments
for f in *phase.fastq
do
	OUT=${f/fastq/bam}
	[ -e $OUT ] || ( bwa mem -t 16 -M -k 16 $REF $f 2> ${f/fastq/log} | samtools view -uS -F 0x104 - | samtools sort -@ 16 -m 8G -o $OUT -O bam -T ${OUT%%.bam} - )
done

# Convert bams to 10kb bedgraph
BED=${REF}.10000.bed
export LC_ALL=C
[ -e $BED ] || bedtools makewindows -g ${REF}.fai -w 10000 | sort -S 20G -k1,1 -k2,2n > $BED
for f in *phase.bam
do
	BG=${f/bam/bedgraph}
	[ -e $BG ] || echo "bedtools bamtobed -i $f | sort -S 2G -k1,1 -k2,2n | bedtools intersect -a $BED -b - -c -sorted > $BG"
done | xargs -L 1 -P 16 -I {} bash -c {}

# Convert read counts to normalized RPKMs
for f in KC167_*phase.bedgraph
do
	OUT=${f%%-phase.bedgraph}_RPKM.bedgraph
	if [ ! -e $OUT ]
	then
		TOTAL=$(awk '{sum+=$4} END {print sum}' $f)
		awk -v sum=$TOTAL '{printf("%s\t%i\t%i\t%f\n", $1, $2, $3, $4*1000000/sum)}' $f > $OUT
	fi
done

# Make weighted replication signal
[ -e weighted_RPKM.bedgraph ] || \
paste KC167_Early_S_RPKM.bedgraph \
	<( cut -f 4 KC167_Mid-Early_S_RPKM.bedgraph ) \
	<( cut -f 4 KC167_Mid-Late_S_RPKM.bedgraph ) \
	<( cut -f 4 KC167_Late_S_RPKM.bedgraph ) | awk '{
s1=$4+$5+$6+$7;
if ( s1 == 0 )
	out=0;
else
	out=0.125*$7/s1 + 0.375*$6/s1 + 0.625*$5/s1 + 0.875*$4/s1;
printf("%s\t%i\t%i\t%f\n", $1, $2, $3, out)
}' > weighted_RPKM.bedgraph

# Smooth weighted rpkm with loess
Rscript - << 'EOF'
D = read.table('weighted_RPKM.bedgraph')
library(stats)
fit = loess(D$V4 ~ as.numeric(rownames(D)), span=20/length(D$V4))
Y=predict(fit)
write.table(cbind(D[1:3], Y), "continuous.bedgraph", quote=F, row.names=F, col.names=F, sep="\t")
EOF

# Create discrete early calls from smoothed, continuous signal
awk '{
if ($4 > 0.5)
	out=1;
else
	out=0;
printf("%s\t%i\t%i\t%i\n", $1, $2, $3, out)
}' continuous.bedgraph > continuous_segmented.bedgraph

###################################
# Run repliscan
###################################

ml use $WORK/public/apps/modulefiles
ml repliscan

# Make totalS
TS=KC167_Total-S.bam
[ -e $TS ] || samtools merge -@ 16 $TS *phase.bam

# Make input configurations
cat << EOF > input_2S.txt
TS	KC167_Total-S.bam
ES	KC167_Early_S-phase.bam	KC167_Mid-Early_S-phase.bam
LS	KC167_Mid-Late_S-phase.bam	KC167_Late_S-phase.bam
EOF

cat << EOF > input_4S.txt
TS	KC167_Total-S.bam
ES	KC167_Early_S-phase.bam
EMS	KC167_Mid-Early_S-phase.bam
MLS	KC167_Mid-Late_S-phase.bam
LS	KC167_Late_S-phase.bam
EOF

# Run Repliscan
for S in 2S 4S
do
	INPUT="../input_${S}.txt"
	for bin in 10
	do
		BED=${REF}.${bin}000.bed
		[ -e $BED ] || ( bedtools makewindows -g ${REF}.fai -w ${bin}000 | sort -S 5G -k1,1 -k2,2n > $BED )
		fold=repliscan_${bin}kb_$S
		[ -e $fold ] && rm -rf $fold
		mkdir $fold && ln -s $PWD/*bam $fold/
		echo "cd $fold && repliscan.py -R norm -a sum -r ../$REF -w ${bin}000 --plot $INPUT &> log.out || echo $fold 'failed' && cp ratio_segmentation.gff3 ../${fold}.gff3"
	done
done | xargs -n 1 -P 4 -I {} bash -c {}
